diff --git a/configure.ac b/configure.ac
index bfda3fa72e9058eccf8834b283d5ec36116bee37..a16131c8bb8de1ab9c71faa40b326ceeb757f9ef 100644
--- a/configure.ac
+++ b/configure.ac
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
-AC_INIT([SWIFT],[0.8.1],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.8.2],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# We want to stop when given unrecognised options. No subdirs so this is safe.
@@ -1371,7 +1371,7 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-du, anarchy-pu default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -1408,6 +1408,9 @@ case "$with_hydro" in
planetary)
AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
;;
+ anarchy-du)
+ AC_DEFINE([ANARCHY_DU_SPH], [1], [ANARCHY (DU) SPH])
+ ;;
anarchy-pu)
AC_DEFINE([ANARCHY_PU_SPH], [1], [ANARCHY (PU) SPH])
;;
diff --git a/doc/RTD/source/conf.py b/doc/RTD/source/conf.py
index 327c5524820772c4c891418cff72a3c22a4fb2b0..468e92d87b05e8d76890c904ec2c9a55b76a6038 100644
--- a/doc/RTD/source/conf.py
+++ b/doc/RTD/source/conf.py
@@ -25,7 +25,7 @@ author = 'SWIFT Team'
# The short X.Y version
version = '0.8'
# The full version, including alpha/beta/rc tags
-release = '0.8.1'
+release = '0.8.2'
# -- General configuration ---------------------------------------------------
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 4dcc1e0381adcb0405dd5b502c3d64896269a10f..8649e34a77f68c2f8b14b4619ce0fd4a616390fc 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: EAGLE-L0012N0188-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -127,8 +131,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index e1d9100a8f87bb8670fcb82ce4451084ecb1375a..c0e6ffe9aabf9d6ea84b034a81543e0b347e3cc5 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: EAGLE-L0025N0376-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -128,8 +132,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 3efe5ba972903dbc3644581d214285603e75ca78..082422ec6acba19acdc02c293a82b0fef3e0e23c 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: EAGLE-L0050N0752-Ref
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
@@ -128,8 +132,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 0067311c413a670ec6696ae5f6f85bc61e2c83b3..d6ecfad02ae93d652cebb7d97933787e24d5ee16 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -135,8 +135,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index c2124f53d9ced16d36b9f74589faae186bdccbd4..a13925e6079ce2aeec3adbdef9a46040fd87677b 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -137,8 +137,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 7805ea93c4fecf8e971d3c6902bb04a0a2e86cb8..06e47d45335380de2abd7edd3017043d735f352a 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -137,8 +137,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index b468e3ad3c7a3cb3497f652debe44677d55d600c..0c993af67652a5c81489e0827ab4e246b586e05b 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -127,6 +127,42 @@ EAGLEEntropyFloor:
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
+# EAGLE feedback model
+EAGLEFeedback:
+ use_SNII_feedback: 1 # Global switch for SNII thermal (stochastic) feedback.
+ use_SNIa_feedback: 1 # Global switch for SNIa thermal (continuous) feedback.
+ use_AGB_enrichment: 1 # Global switch for enrichement from AGB stars.
+ use_SNII_enrichment: 1 # Global switch for enrichement from SNII stars.
+ use_SNIa_enrichment: 1 # Global switch for enrichement from SNIa stars.
+ filename: ./yieldtables/ # Path to the directory containing the EAGLE yield tables.
+ IMF_min_mass_Msun: 0.1 # Minimal stellar mass considered for the Chabrier IMF in solar masses.
+ IMF_max_mass_Msun: 100.0 # Maximal stellar mass considered for the Chabrier IMF in solar masses.
+ SNII_min_mass_Msun: 6.0 # Minimal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
+ SNII_max_mass_Msun: 100.0 # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
+ SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
+ SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
+ SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
+ SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
+ SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
+ SNII_energy_fraction_n_n: 0.8686 # Power-law for the birth density dependance of the SNII energy fraction.
+ SNIa_max_mass_Msun: 8.0 # Maximal mass considered for SNIa feedback and enrichment in solar masses.
+ SNIa_timescale_Gyr: 2.0 # Time-scale of the exponential decay of the SNIa rates in Gyr.
+ SNIa_efficiency_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
+ SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
+ AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
+ SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
+ SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
+ SNII_yield_factor_Nitrogen: 1.0 # (Optional) Correction factor to apply to the Nitrogen yield from the SNII channel.
+ SNII_yield_factor_Oxygen: 1.0 # (Optional) Correction factor to apply to the Oxygen yield from the SNII channel.
+ SNII_yield_factor_Neon: 1.0 # (Optional) Correction factor to apply to the Neon yield from the SNII channel.
+ SNII_yield_factor_Magnesium: 2.0 # (Optional) Correction factor to apply to the Magnesium yield from the SNII channel.
+ SNII_yield_factor_Silicon: 1.0 # (Optional) Correction factor to apply to the Silicon yield from the SNII channel.
+ SNII_yield_factor_Iron: 0.5 # (Optional) Correction factor to apply to the Iron yield from the SNII channel.
+
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
@@ -135,4 +171,4 @@ EAGLEAGN:
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
-
\ No newline at end of file
+
diff --git a/examples/HydroTests/SodShock_1D/plotSolution.py b/examples/HydroTests/SodShock_1D/plotSolution.py
index 12ae8a9cf26fb715281ed00bf565c5c8d9a234fb..a7e6d374bac616440dace666b85c3e7ade479bcd 100644
--- a/examples/HydroTests/SodShock_1D/plotSolution.py
+++ b/examples/HydroTests/SodShock_1D/plotSolution.py
@@ -1,6 +1,7 @@
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 2019 Josh Borrow (josh.borrow@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
@@ -270,7 +271,7 @@ ylim(0.01, 1.1)
# Internal energy profile -------------------------
subplot(234)
-scatter(x, u, marker='.', c=alpha_diff, s=4.0)
+plot(x, u, '.', color='r', ms=4.0)
plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
diff --git a/examples/HydroTests/SodShock_1D/sodShock.yml b/examples/HydroTests/SodShock_1D/sodShock.yml
index b936f50f6c2c7d9078c49fbad868aa5334498957..69554b4db733166fc5dbb6d198966fd8f9b8d49c 100644
--- a/examples/HydroTests/SodShock_1D/sodShock.yml
+++ b/examples/HydroTests/SodShock_1D/sodShock.yml
@@ -27,10 +27,6 @@ Statistics:
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- viscosity_alpha_min: 0.01
- viscosity_alpha: 0.01
- viscosity_alpha_max: 2.0
- viscosity_length: 0.02
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
index e4e084f30020872befa7641f59c876dbcf2ef669..3018f77f3f920e37ed55ee1206eb82e7498f4823 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
@@ -123,8 +123,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/PMillennium/PMillennium-1536/p-mill-1536.yml b/examples/PMillennium/PMillennium-1536/p-mill-1536.yml
index ea44572f90eb69330c973946cb5533a4c58b2c82..770cc4a03e8912eed51d1851b81578496a42070c 100644
--- a/examples/PMillennium/PMillennium-1536/p-mill-1536.yml
+++ b/examples/PMillennium/PMillennium-1536/p-mill-1536.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Planck-Millennium simulation - 1536^3
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
diff --git a/examples/PMillennium/PMillennium-384/p-mill-384.yml b/examples/PMillennium/PMillennium-384/p-mill-384.yml
index 0040f580c8e67182949875d7d85cee2f75851654..428f893300bf3431f20518a6a1259651b3c67226 100644
--- a/examples/PMillennium/PMillennium-384/p-mill-384.yml
+++ b/examples/PMillennium/PMillennium-384/p-mill-384.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Planck-Millennium simulation - 384^3
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
diff --git a/examples/PMillennium/PMillennium-768/p-mill-768.yml b/examples/PMillennium/PMillennium-768/p-mill-768.yml
index 5ba97af72513dca5fdccd216f0ede78d0e279b0a..339fad2c7059d8c8bcab2878fc8958a3bd1977d7 100644
--- a/examples/PMillennium/PMillennium-768/p-mill-768.yml
+++ b/examples/PMillennium/PMillennium-768/p-mill-768.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Planck-Millennium simulation - 768^3
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun
diff --git a/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml b/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml
index 35a6e2762f91707ed43bd5f0107c3dd8a53e12e4..1221cefb92e51e54bfe73a0414b95a89fcda592a 100644
--- a/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml
+++ b/examples/SantaBarbara/SantaBarbara-128/santa_barbara.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Santa-Barbara cluster comparison project 128^3
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
diff --git a/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml b/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml
index 9f6585e0c589d3d016921eb85dd6af13a0643784..04bb18191b750aa996e43b58b24451135bceedca 100644
--- a/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml
+++ b/examples/SantaBarbara/SantaBarbara-256/santa_barbara.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Santa-Barbara cluster comparison project 256^3
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 Msun
diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index a5e2845a6bbe6923c96303254cb6ba8906c5dce8..669d9e3f7c78e0ebf2b68ff0c3a3bcb5369b4a50 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -100,8 +100,8 @@ EAGLEFeedback:
SNII_wind_delay_Gyr: 0.03 # Time in Gyr between a star's birth and the SNII thermal feedback event.
SNII_delta_T_K: 3.16228e7 # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
SNII_energy_erg: 1.0e51 # Energy of one SNII explosion in ergs.
- SNII_energy_fraction_min: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
- SNII_energy_fraction_max: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_min: 0.3 # Minimal fraction of energy applied in a SNII feedback event.
+ SNII_energy_fraction_max: 3.0 # Maximal fraction of energy applied in a SNII feedback event.
SNII_energy_fraction_Z_0: 0.0012663729 # Pivot point for the metallicity dependance of the SNII energy fraction (metal mass fraction).
SNII_energy_fraction_n_0_H_p_cm3: 0.67 # Pivot point for the birth density dependance of the SNII energy fraction in cm^-3.
SNII_energy_fraction_n_Z: 0.8686 # Power-law for the metallicity dependance of the SNII energy fraction.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 375c5abc22d117b7322dd49e37109430b3ef0e32..61e7ab470292d0575881ae0382103988e4748c74 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -1,3 +1,7 @@
+# Define some meta-data about the simulation
+MetaData:
+ run_name: Name of the sim in less than 256 characters. # The name of the simulation. This is written into the snapshot headers.
+
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
diff --git a/src/Makefile.am b/src/Makefile.am
index 3cb445ffe7db5afb1f34dac72cfcbd8b2d1bffec..effb989c3d3b6e204f62dab4d95e66a7289799a8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -95,15 +95,31 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
gravity/Potential/gravity_debug.h gravity/Potential/gravity_part.h \
equation_of_state.h \
equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
- hydro.h hydro_io.h \
+ hydro.h hydro_io.h hydro_parameters.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
+ hydro/Minimal/hydro_parameters.h \
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
+ hydro/Default/hydro_parameters.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
+ hydro/Gadget2/hydro_parameters.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
+ hydro/PressureEntropy/hydro_parameters.h \
+ hydro/PressureEnergy/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+ hydro/PressureEnergy/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+ hydro/PressureEnergy/hydro_parameters.h \
+ hydro/PressureEnergyMorrisMonaghanAV/hydro.h hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h \
+ hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h \
+ hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
+ hydro/AnarchyPU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+ hydro/AnarchyPU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+ hydro/AnarchyPU/hydro_parameters.h \
+ hydro/AnarchyDU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
+ hydro/AnarchyDU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
+ hydro/AnarchyDU/hydro_parameters.h \
hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
hydro/GizmoMFV/hydro_part.h \
@@ -115,6 +131,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFV/hydro_slope_limiters.h \
hydro/GizmoMFV/hydro_unphysical.h \
hydro/GizmoMFV/hydro_velocities.h \
+ hydro/GizmoMFV/hydro_parameters.h \
hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
hydro/GizmoMFM/hydro_part.h \
@@ -125,6 +142,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFM/hydro_slope_limiters_face.h \
hydro/GizmoMFM/hydro_slope_limiters.h \
hydro/GizmoMFM/hydro_unphysical.h \
+ hydro/GizmoMFM/hydro_parameters.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
@@ -141,6 +159,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Shadowswift/voronoi3d_cell.h \
hydro/Shadowswift/voronoi_algorithm.h \
hydro/Shadowswift/voronoi_cell.h \
+ hydro/Shadowswift/hydro_parameters.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \
riemann/riemann_checks.h \
diff --git a/src/cell.c b/src/cell.c
index 86a0bee9ee65499f7aa4fadd471e8f7a87ac0162..d2b102a974ff92d0548fe072b70c44ce417343cf 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -64,6 +64,7 @@
#include "scheduler.h"
#include "space.h"
#include "space_getsid.h"
+#include "star_formation.h"
#include "stars.h"
#include "timers.h"
#include "tools.h"
@@ -4207,6 +4208,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
if (part_is_active(p, e)) {
hydro_init_part(p, &e->s->hs);
chemistry_init_part(p, e->chemistry);
+ star_formation_init_part(p, e->star_formation);
tracers_after_init(p, xp, e->internal_units, e->physical_constants,
with_cosmology, e->cosmology, e->hydro_properties,
e->cooling_func, e->time);
diff --git a/src/const.h b/src/const.h
index 613a48920e6f26c209faf6e354b82c2ed5be0bf1..d7017699cff56b1de545e6e0ae78013f000b76d0 100644
--- a/src/const.h
+++ b/src/const.h
@@ -20,16 +20,6 @@
#ifndef SWIFT_CONST_H
#define SWIFT_CONST_H
-/* SPH Viscosity constants. */
-/* Cosmology default beta=3.0. Planetary default beta=4.0
- * Alpha can be set in the parameter file.
- * Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_beta 3.0f
-
-/* SPH Thermal conductivity constants. */
-#define const_conductivity_alpha \
- 1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
-
/* Time integration constants. */
#define const_max_u_change 0.1f
diff --git a/src/debug.c b/src/debug.c
index eec70c3eede90aa0fdd275f51cf18422f7ba233a..f8509af805fed90cbdbdfd58849cb7b6352062e9 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -60,6 +60,8 @@
#include "./hydro/Shadowswift/hydro_debug.h"
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_debug.h"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro_debug.h"
#elif defined(ANARCHY_PU_SPH)
#include "./hydro/AnarchyPU/hydro_debug.h"
#else
diff --git a/src/engine.c b/src/engine.c
index 3d51a4a7510de4950b53247bf226ba4ba9018269..e760293d954d7ec521286774a3b445c27bc545fb 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5003,6 +5003,15 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
/* Make the space link back to the engine. */
s->e = e;
+ /* Read the run label */
+ memset(e->run_name, 0, PARSER_MAX_LINE_SIZE);
+ parser_get_opt_param_string(params, "MetaData:run_name", e->run_name,
+ "Untitled SWIFT simulation");
+ if (strlen(e->run_name) == 0) {
+ error("The run name in the parameter file cannot be an empty string.");
+ }
+ if (e->nodeID == 0) message("Running simulation '%s'.", e->run_name);
+
/* Setup the timestep if non-cosmological */
if (!(e->policy & engine_policy_cosmology)) {
e->time_begin =
@@ -6055,14 +6064,20 @@ void engine_recompute_displacement_constraint(struct engine *e) {
* @brief Frees up the memory allocated for this #engine
*/
void engine_clean(struct engine *e) {
+ /* Start by telling the runners to stop. */
+ e->step_props = engine_step_prop_done;
+ swift_barrier_wait(&e->run_barrier);
- for (int i = 0; i < e->nr_threads; ++i) {
+ /* Wait for each runner to come home. */
+ for (int k = 0; k < e->nr_threads; k++) {
+ if (pthread_join(e->runners[k].thread, /*retval=*/NULL) != 0)
+ error("Failed to join runner %i.", k);
#ifdef WITH_VECTORIZATION
- cache_clean(&e->runners[i].ci_cache);
- cache_clean(&e->runners[i].cj_cache);
+ cache_clean(&e->runners[k].ci_cache);
+ cache_clean(&e->runners[k].cj_cache);
#endif
- gravity_cache_clean(&e->runners[i].ci_gravity_cache);
- gravity_cache_clean(&e->runners[i].cj_gravity_cache);
+ gravity_cache_clean(&e->runners[k].ci_gravity_cache);
+ gravity_cache_clean(&e->runners[k].cj_gravity_cache);
}
swift_free("runners", e->runners);
free(e->snapshot_units);
diff --git a/src/engine.h b/src/engine.h
index 1c788c3e60a05aa9394fe35f5bbf245d54e8244b..708dd2e3553a4594000d8175be706031e1c71857 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -96,7 +96,8 @@ enum engine_step_properties {
engine_step_prop_snapshot = (1 << 4),
engine_step_prop_restarts = (1 << 5),
engine_step_prop_stf = (1 << 6),
- engine_step_prop_logger_index = (1 << 7)
+ engine_step_prop_logger_index = (1 << 7),
+ engine_step_prop_done = (1 << 8)
};
/* Some constants */
@@ -458,6 +459,9 @@ struct engine {
/* Maximum number of tasks needed for restarting. */
int restart_max_tasks;
+
+ /* Label of the run */
+ char run_name[PARSER_MAX_LINE_SIZE];
};
/* Function prototypes, engine.c. */
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index f7e714b21efdf6987e9a683a0873211eea8e13b7..87f8e5cc48614b44660b87529a31d64aaa354ca3 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -841,6 +841,11 @@ void feedback_props_init(struct feedback_props* fp,
fp->imf_min_mass_msun =
parser_get_param_double(params, "EAGLEFeedback:IMF_min_mass_Msun");
+ /* Check that it makes sense. */
+ if (fp->imf_max_mass_msun < fp->imf_min_mass_msun) {
+ error("Can't have the max IMF mass smaller than the min IMF mass!");
+ }
+
fp->log10_imf_max_mass_msun = log10(fp->imf_max_mass_msun);
fp->log10_imf_min_mass_msun = log10(fp->imf_min_mass_msun);
@@ -870,6 +875,11 @@ void feedback_props_init(struct feedback_props* fp,
const double SNII_max_mass_msun =
parser_get_param_double(params, "EAGLEFeedback:SNII_max_mass_Msun");
+ /* Check that it makes sense. */
+ if (SNII_max_mass_msun < SNII_min_mass_msun) {
+ error("Can't have the max SNII mass smaller than the min SNII mass!");
+ }
+
fp->log10_SNII_min_mass_msun = log10(SNII_min_mass_msun);
fp->log10_SNII_max_mass_msun = log10(SNII_max_mass_msun);
@@ -887,6 +897,11 @@ void feedback_props_init(struct feedback_props* fp,
fp->n_Z =
parser_get_param_double(params, "EAGLEFeedback:SNII_energy_fraction_n_Z");
+ /* Check that it makes sense. */
+ if (fp->f_E_max < fp->f_E_min) {
+ error("Can't have the maximal energy fraction smaller than the minimal!");
+ }
+
/* Properties of the SNII enrichment model -------------------------------- */
/* Set factors for each element adjusting SNII yield */
diff --git a/src/hydro.h b/src/hydro.h
index 3bf7d2228b528796a8717d6a5ab17fda6f569d25..e5aedeb0c00b92e8bdbf3b707516a90132274f65 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -72,6 +72,10 @@
#include "./hydro/Planetary/hydro.h"
#include "./hydro/Planetary/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro.h"
+#include "./hydro/AnarchyDU/hydro_iact.h"
+#define SPH_IMPLEMENTATION "ANARCHY (Density-Energy) SPH (Borrow+ in prep)"
#elif defined(ANARCHY_PU_SPH)
#include "./hydro/AnarchyPU/hydro.h"
#include "./hydro/AnarchyPU/hydro_iact.h"
diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f14e583f6841ffaf73840b8c04bc6bc5ec96ddf
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro.h
@@ -0,0 +1,1050 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_DU_HYDRO_H
+#define SWIFT_ANARCHY_DU_HYDRO_H
+
+/**
+ * @file AnarchyDU/hydro.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (Non-neighbour loop
+ * equations)
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "dimension.h"
+#include "entropy_floor.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+#include "./hydro_parameters.h"
+
+#include <float.h>
+
+/**
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp) {
+
+ return xp->u_full;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ return xp->u_full * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+ const struct part *restrict p) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_pressure * hydro_get_comoving_pressure(p);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+ const struct part *restrict p, const struct xpart *restrict xp) {
+
+ return gas_entropy_from_internal_energy(p->rho, xp->u_full);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, xp->u_full);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
+
+ return gas_soundspeed_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p);
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+ const struct part *restrict p) {
+
+ return p->rho;
+}
+
+/**
+ * @brief Returns the comoving density of a particle.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a3_inv * p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part *restrict p) {
+
+ return p->mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
+/**
+ * @brief Returns the velocities drifted to the current time of a particle.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle.
+ * @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
+ * @param dt_kick_grav The time (for gravity accelerations) since the last kick.
+ * @param v (return) The velocities at the current time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
+ const struct part *restrict p, const struct xpart *xp, float dt_kick_hydro,
+ float dt_kick_grav, float v[3]) {
+
+ v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
+ xp->a_grav[0] * dt_kick_grav;
+ v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
+ xp->a_grav[1] * dt_kick_grav;
+ v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
+ xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
+
+ return p->u_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ * @param cosmo Cosmology data structure
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy_dt(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u_dt * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
+
+ p->u_dt = du_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy_dt(struct part *restrict p,
+ const struct cosmology *cosmo,
+ float du_dt) {
+
+ p->u_dt = du_dt / cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the physical entropy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param entropy The physical entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
+ struct part *p, struct xpart *xp, const struct cosmology *cosmo,
+ const float entropy) {
+
+ /* Note there is no conversion from physical to comoving entropy */
+ const float comoving_entropy = entropy;
+ xp->u_full = gas_internal_energy_from_entropy(p->rho, comoving_entropy);
+}
+
+/**
+ * @brief Sets the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy(struct part *p, struct xpart *xp,
+ const struct cosmology *cosmo,
+ const float u) {
+
+ xp->u_full = u / cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the drifted physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_drifted_physical_internal_energy(struct part *p,
+ const struct cosmology *cosmo,
+ const float u) {
+
+ p->u = u / cosmo->a_factor_internal_energy;
+
+ /* Now recompute the extra quantities */
+
+ /* Compute the sound speed */
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float pressure = hydro_get_comoving_pressure(p);
+
+ /* Update variables. */
+ p->force.soundspeed = soundspeed;
+ p->force.pressure = pressure;
+}
+
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ p->viscosity.alpha = alpha;
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ hydro_set_viscosity_alpha(p,
+ hydro_props_default_viscosity_alpha_feedback_reset);
+}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties,
+ const struct cosmology *restrict cosmo) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+ (cosmo->a_factor_sound_speed * p->viscosity.v_sig);
+
+ return dt_cfl;
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
+/**
+ * @brief Operations performed when a particle gets removed from the
+ * simulation volume.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void hydro_remove_part(
+ const struct part *p, const struct xpart *xp) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part *restrict p, const struct hydro_space *hs) {
+
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->rho = 0.f;
+ p->density.rho_dh = 0.f;
+
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+
+ p->viscosity.div_v = 0.f;
+ p->diffusion.laplace_u = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
+
+ const float rho_inv = 1.f / p->rho;
+ const float a_inv2 = cosmo->a2_inv;
+
+ /* Finish calculation of the velocity curl components */
+ p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+
+ /* Finish calculation of the velocity divergence */
+ p->viscosity.div_v *= h_inv_dim_plus_one * rho_inv * a_inv2;
+ p->viscosity.div_v += cosmo->H * hydro_dimension;
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * This function is called after the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ const float fac_B = cosmo->a_factor_Balsara_eps;
+
+ /* Compute the norm of the curl */
+ const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the norm of div v */
+ const float abs_div_v = fabsf(p->viscosity.div_v);
+
+ /* Compute the sound speed -- see theory section for justification */
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float pressure = hydro_get_comoving_pressure(p);
+
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
+
+ /* Compute the "grad h" term */
+ const float rho_inv = 1.f / p->rho;
+ const float grad_h_term =
+ 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+
+ /* Update variables. */
+ p->force.f = grad_h_term;
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+ p->force.balsara = balsara;
+}
+
+/**
+ * @brief Resets the variables that are required for a gradient calculation.
+ *
+ * This function is called after hydro_prepare_gradient.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_gradient(
+ struct part *restrict p) {
+
+ p->viscosity.v_sig = 2.f * p->force.soundspeed;
+}
+
+/**
+ * @brief Finishes the gradient calculation.
+ *
+ * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
+ * in which case no gradients are used.
+ *
+ * This method also initializes the force loop variables.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part *p) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Include the extra factors in the del^2 u */
+
+ p->diffusion.laplace_u *= 2.f * h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->rho = p->mass * kernel_root * h_inv_dim;
+ p->viscosity.v_sig = 0.f;
+ p->density.wcount = kernel_root * h_inv_dim;
+ p->density.rho_dh = 0.f;
+ p->density.wcount_dh = 0.f;
+
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+
+ /* Probably not shocking, so this is safe to do */
+ p->viscosity.div_v = 0.f;
+ p->diffusion.laplace_u = 0.f;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ * @param dt_alpha The time-step used to evolve non-cosmological quantities such
+ * as the artificial viscosity.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const float dt_alpha) {
+
+ /* Here we need to update the artificial viscosity */
+
+ /* We use in this function that h is the radius of support */
+ const float kernel_support_physical = p->h * cosmo->a * kernel_gamma;
+ const float kernel_support_physical_inv = 1.f / kernel_support_physical;
+ const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed;
+ const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo);
+
+ const float sound_crossing_time_inverse =
+ soundspeed_physical * kernel_support_physical_inv;
+
+ /* Construct time differential of div.v implicitly following the ANARCHY spec
+ */
+
+ const float div_v_dt =
+ dt_alpha == 0.f
+ ? 0.f
+ : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha;
+
+ /* Construct the source term for the AV; if shock detected this is _positive_
+ * as div_v_dt should be _negative_ before the shock hits */
+ const float S = kernel_support_physical * kernel_support_physical *
+ max(0.f, -1.f * div_v_dt);
+ /* 0.25 factor comes from our definition of v_sig (sum of soundspeeds rather
+ * than mean). */
+ /* Note this is v_sig_physical squared, not comoving */
+ const float v_sig_square = 0.25 * v_sig_physical * v_sig_physical;
+
+ /* Calculate the current appropriate value of the AV based on the above */
+ const float alpha_loc =
+ hydro_props->viscosity.alpha_max * S / (v_sig_square + S);
+
+ if (alpha_loc > p->viscosity.alpha) {
+ /* Reset the value of alpha to the appropriate value */
+ p->viscosity.alpha = alpha_loc;
+ } else {
+ /* Integrate the alpha forward in time to decay back to alpha = alpha_loc */
+ p->viscosity.alpha =
+ alpha_loc + (p->viscosity.alpha - alpha_loc) *
+ expf(-dt_alpha * sound_crossing_time_inverse *
+ hydro_props->viscosity.length);
+ }
+
+ /* Check that we did not hit the minimum */
+ p->viscosity.alpha =
+ max(p->viscosity.alpha, hydro_props->viscosity.alpha_min);
+
+ /* Set our old div_v to the one for the next loop */
+ p->viscosity.div_v_previous_step = p->viscosity.div_v;
+
+ /* Now for the diffusive alpha */
+
+ const float diffusion_timescale_physical_inverse =
+ v_sig_physical * kernel_support_physical_inv;
+
+ const float sqrt_u_inv = 1.f / sqrtf(p->u);
+
+ /* Calculate initial value of alpha dt before bounding */
+ /* Evolution term: following Schaller+ 2015. This is made up of several
+ cosmology factors: physical smoothing length, sound speed from laplace(u) /
+ sqrt(u), and the 1 / a^2 coming from the laplace operator. */
+ float alpha_diff_dt = hydro_props->diffusion.beta * kernel_support_physical *
+ p->diffusion.laplace_u * cosmo->a_factor_sound_speed *
+ sqrt_u_inv * cosmo->a2_inv;
+
+ /* Decay term: not documented in Schaller+ 2015 but was present
+ * in the original EAGLE code and in appendix of Schaye+ 2015 */
+ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) *
+ diffusion_timescale_physical_inverse;
+
+ float new_diffusion_alpha = p->diffusion.alpha;
+ new_diffusion_alpha += alpha_diff_dt * dt_alpha;
+
+ /* Consistency checks to ensure min < alpha < max */
+ new_diffusion_alpha =
+ min(new_diffusion_alpha, hydro_props->diffusion.alpha_max);
+ new_diffusion_alpha =
+ max(new_diffusion_alpha, hydro_props->diffusion.alpha_min);
+
+ p->diffusion.alpha = new_diffusion_alpha;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part *restrict p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->u_dt = 0.0f;
+ p->force.h_dt = 0.0f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+ struct part *restrict p, const struct xpart *restrict xp) {
+
+ /* Re-set the predicted velocities */
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
+ /* Re-set the entropy */
+ p->u = xp->u_full;
+
+ /* Compute the sound speed */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
+ * @param cosmo The cosmological model.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param floor_props The properties of the entropy floor.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+ float dt_therm, const struct cosmology *cosmo,
+ const struct hydro_props *hydro_props,
+ const struct entropy_floor_properties *floor_props) {
+
+ /* Predict the internal energy */
+ p->u += p->u_dt * dt_therm;
+
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, floor_u);
+ p->u = max(p->u, min_u);
+
+ const float h_inv = 1.f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density and weighted pressure */
+ const float w2 = -hydro_dimension * w1;
+ if (fabsf(w2) < 0.2f) {
+ const float expf_approx =
+ approx_expf(w2); /* 4th order expansion of exp(w) */
+ p->rho *= expf_approx;
+ } else {
+ const float expf_exact = expf(w2);
+ p->rho *= expf_exact;
+ }
+
+ /* Compute the new sound speed */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param dt_grav The time-step for this kick (for gravity quantities).
+ * @param dt_hydro The time-step for this kick (for hydro quantities).
+ * @param dt_kick_corr The time-step for this kick (for gravity corrections).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ * @param floor_props The properties of the entropy floor.
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ float dt_grav, float dt_hydro, float dt_kick_corr,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const struct entropy_floor_properties *floor_props) {
+
+ /* Integrate the internal energy forward in time */
+ const float delta_u = p->u_dt * dt_therm;
+
+ /* Do not decrease the energy by more than a factor of 2*/
+ xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full);
+
+ /* Check against entropy floor */
+ const float floor_A = entropy_floor(p, cosmo, floor_props);
+ const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ /* Take highest of both limits */
+ const float energy_min = max(min_u, floor_u);
+
+ if (xp->u_full < energy_min) {
+ xp->u_full = energy_min;
+ p->u_dt = 0.f;
+ }
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+ /* Convert the physcial internal energy to the comoving one. */
+ /* u' = a^(3(g-1)) u */
+ const float factor = 1.f / cosmo->a_factor_internal_energy;
+ p->u *= factor;
+ xp->u_full = p->u;
+
+ /* Apply the minimal energy limit */
+ const float min_comoving_energy =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+ if (xp->u_full < min_comoving_energy) {
+ xp->u_full = min_comoving_energy;
+ p->u = min_comoving_energy;
+ p->u_dt = 0.f;
+ }
+
+ /* Set the initial value of the artificial viscosity based on the non-variable
+ schemes for safety */
+
+ p->viscosity.alpha = hydro_props->viscosity.alpha;
+ /* Initialise this here to keep all the AV variables together */
+ p->viscosity.div_v_previous_step = 0.f;
+
+ /* Set the initial values for the thermal diffusion */
+ p->diffusion.alpha = hydro_props->diffusion.alpha;
+
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->time_bin = 0;
+ p->wakeup = time_bin_not_awake;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+ xp->a_grav[0] = 0.f;
+ xp->a_grav[1] = 0.f;
+ xp->a_grav[2] = 0.f;
+ xp->u_full = p->u;
+
+ hydro_reset_acceleration(p);
+ hydro_init_part(p, NULL);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_H */
diff --git a/src/hydro/AnarchyDU/hydro_debug.h b/src/hydro/AnarchyDU/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..bc72bec62f8ddff6086516a32a4b86651e6aa66a
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro_debug.h
@@ -0,0 +1,47 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_ANARCHY_DU_HYDRO_DEBUG_H
+#define SWIFT_ANARCHY_DU_HYDRO_DEBUG_H
+
+/**
+ * @file AnarchyDU/hydro_debug.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (Debugging routines)
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
+ "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+ "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+ "alpha=%.3e \n"
+ "time_bin=%d\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+ xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+ p->u, p->u_dt, p->viscosity.v_sig, hydro_get_comoving_pressure(p), p->h,
+ p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+ p->viscosity.alpha, p->time_bin);
+}
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_DEBUG_H */
diff --git a/src/hydro/AnarchyDU/hydro_iact.h b/src/hydro/AnarchyDU/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..d1d3b1f6c0e4941539044d14bda1dd9f7a5ca615
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro_iact.h
@@ -0,0 +1,583 @@
+/*******************************************************************************
+ * This file is part* of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_DU_HYDRO_IACT_H
+#define SWIFT_ANARCHY_DU_HYDRO_IACT_H
+
+/**
+ * @file AnarchyDU/hydro_iact.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (interaction routines)
+ */
+
+#include "adiabatic_index.h"
+#include "minmax.h"
+
+#include "./hydro_parameters.h"
+
+/**
+ * @brief Density interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, const float* dx, float hi, float hj, struct part* pi,
+ struct part* pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+ float dv[3], curlvr[3];
+
+ const float r = sqrtf(r2);
+
+ /* Get the masses. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ pj->rho += mi * wj;
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+ /* Now we need to compute the div terms */
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+ const float facj = mi * wj_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->viscosity.div_v -= faci * dvdr;
+ pj->viscosity.div_v -= facj * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+
+ /* Negative because of the change in sign of dx & dv. */
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
+}
+
+/**
+ * @brief Density interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, const float* dx, float hi, float hj, struct part* pi,
+ const struct part* pj, float a, float H) {
+
+ float wi, wi_dx;
+ float dv[3], curlvr[3];
+
+ /* Get the masses. */
+ const float mj = pj->mass;
+
+ /* Get r and r inverse. */
+ const float r = sqrtf(r2);
+
+ const float h_inv = 1.f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ const float r_inv = 1.f / r;
+ const float faci = mj * wi_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->viscosity.div_v -= faci * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+ float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+ struct part* restrict pj, float a, float H) {
+
+ /* We need to construct the maximal signal velocity between our particle
+ * and all of it's neighbours */
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.f / r;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Cosmology terms for the signal velocity */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij;
+
+ /* Update if we need to */
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
+ pj->viscosity.v_sig = max(pj->viscosity.v_sig, new_v_sig);
+
+ /* Calculate Del^2 u for the thermal diffusion coefficient. */
+ /* Need to get some kernel values F_ij = wi_dx */
+ float wi, wi_dx, wj, wj_dx;
+
+ const float ui = r / hi;
+ const float uj = r / hj;
+
+ kernel_deval(ui, &wi, &wi_dx);
+ kernel_deval(uj, &wj, &wj_dx);
+
+ const float delta_u_factor = (pi->u - pj->u) * r_inv;
+ pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
+ pj->diffusion.laplace_u -= pi->mass * delta_u_factor * wj_dx / pi->rho;
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+ float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+ struct part* restrict pj, float a, float H) {
+
+ /* We need to construct the maximal signal velocity between our particle
+ * and all of it's neighbours */
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.f / r;
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Cosmology terms for the signal velocity */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Add Hubble flow */
+
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float new_v_sig = ci + cj - const_viscosity_beta * mu_ij;
+
+ /* Update if we need to */
+ pi->viscosity.v_sig = max(pi->viscosity.v_sig, new_v_sig);
+
+ /* Calculate Del^2 u for the thermal diffusion coefficient. */
+ /* Need to get some kernel values F_ij = wi_dx */
+ float wi, wi_dx;
+
+ const float ui = r / hi;
+
+ kernel_deval(ui, &wi, &wi_dx);
+
+ const float delta_u_factor = (pi->u - pj->u) * r_inv;
+ pi->diffusion.laplace_u += pj->mass * delta_u_factor * wi_dx / pj->rho;
+}
+
+/**
+ * @brief Force interaction between two part*icles.
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, const float* dx, float hi, float hj, struct part* pi,
+ struct part* pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ const float mj = pj->mass;
+ const float mi = pi->mass;
+
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Includes the hubble flow term; not used for du/dt */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
+ const_viscosity_beta * mu_ij;
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = rhoi + rhoj;
+ const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
+ const float visc =
+ -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* Compute gradient terms */
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+ const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+ /* Diffusion term */
+ const float v_diff =
+ max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+ const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+ /* wi_dx + wj_dx / 2 is F_ij */
+ const float diff_du_term =
+ alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
+ const float du_dt_j = sph_du_term_j + visc_du_term - diff_du_term;
+
+ /* Internal energy time derivative */
+ pi->u_dt += du_dt_i * mj;
+ pj->u_dt += du_dt_j * mi;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+}
+
+/**
+ * @brief Force interaction between two part*icles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, const float* dx, float hi, float hj, struct part* pi,
+ const struct part* pj, float a, float H) {
+
+ /* Cosmological factors entering the EoMs */
+ const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+ const float a2_Hubble = a * a * H;
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Recover some data */
+ // const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Includes the hubble flow term; not used for du/dt */
+ const float dvdr_Hubble = dvdr + a2_Hubble * r2;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = min(dvdr_Hubble, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
+ const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
+ const_viscosity_beta * mu_ij;
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = rhoi + rhoj;
+ const float alpha = pi->viscosity.alpha + pj->viscosity.alpha;
+ const float visc =
+ -0.25f * alpha * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
+
+ /* Compute gradient terms */
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
+
+ /* SPH acceleration term */
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ /* Get the time derivative for u. */
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
+
+ /* Diffusion term */
+ const float v_diff =
+ max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+ const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+ /* wi_dx + wj_dx / 2 is F_ij */
+ const float diff_du_term =
+ alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
+
+ /* Internal energy time derivative */
+ pi->u_dt += du_dt_i * mj;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+}
+
+/**
+ * @brief Timestep limiter loop
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_limiter(
+ float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+ struct part* restrict pj, float a, float H) {
+
+ /* Nothing to do here if both particles are active */
+}
+
+/**
+ * @brief Timestep limiter loop (non-symmetric version)
+ *
+ * @param r2 Comoving square distance between the two part*icles.
+ * @param dx Comoving vector separating both part*icles (pi - pj).
+ * @param hi Comoving smoothing-length of part*icle i.
+ * @param hj Comoving smoothing-length of part*icle j.
+ * @param pi First part*icle.
+ * @param pj Second part*icle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_limiter(
+ float r2, const float* dx, float hi, float hj, struct part* restrict pi,
+ struct part* restrict pj, float a, float H) {
+
+ /* Wake up the neighbour? */
+ if (pi->viscosity.v_sig >
+ const_limiter_max_v_sig_ratio * pj->viscosity.v_sig) {
+
+ pj->wakeup = time_bin_awake;
+ }
+}
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_IACT_H */
diff --git a/src/hydro/AnarchyDU/hydro_io.h b/src/hydro/AnarchyDU/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..fcec103b1e634a793268851935392e51a1611a97
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro_io.h
@@ -0,0 +1,223 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR DURPOSE. 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_ANARCHY_DU_HYDRO_IO_H
+#define SWIFT_ANARCHY_DU_HYDRO_IO_H
+
+/**
+ * @file AnarchyDU/hydro_io.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (i/o routines)
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+#include "./hydro_parameters.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void hydro_read_particles(struct part* parts,
+ struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, rho);
+}
+
+INLINE static void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_entropy(p, xp);
+}
+
+INLINE static void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_pressure(p);
+}
+
+INLINE static void convert_part_pos(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = p->x[0];
+ ret[1] = p->x[1];
+ ret[2] = p->x[2];
+ }
+}
+
+INLINE static void convert_part_vel(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+INLINE static void convert_part_potential(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
+INLINE static void convert_viscosity(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+ ret[0] = p->viscosity.alpha;
+}
+
+INLINE static void convert_diffusion(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+ ret[0] = p->diffusion.alpha;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void hydro_write_particles(const struct part* parts,
+ const struct xpart* xparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 12;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[7] = io_make_output_field_convert_part(
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
+ list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, xparts, convert_S);
+ list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
+ list[10] = io_make_output_field_convert_part("Viscosity", FLOAT, 1,
+ UNIT_CONV_NO_UNITS, parts,
+ xparts, convert_viscosity);
+ list[11] = io_make_output_field_convert_part("Diffusion", FLOAT, 1,
+ UNIT_CONV_NO_UNITS, parts,
+ xparts, convert_diffusion);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+INLINE static void hydro_write_flavour(hid_t h_grpsph) {
+
+ /* Viscosity and thermal conduction */
+ /* Nothing in this minimal model... */
+ io_write_attribute_s(h_grpsph, "Thermal Conductivity Model",
+ "Simple treatment as in Price (2008)");
+ io_write_attribute_s(h_grpsph, "Viscosity Model",
+ "Simplified version of Cullen & Denhen (2011)");
+
+ /* Time integration properties */
+ io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+INLINE static int writeEntropyFlag(void) { return 0; }
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_IO_H */
diff --git a/src/hydro/AnarchyDU/hydro_parameters.h b/src/hydro/AnarchyDU/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..86935344d940c90d542fedf03d2e9347e31c3222
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro_parameters.h
@@ -0,0 +1,298 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_ANARCHY_DU_HYDRO_PARAMETERS_H
+#define SWIFT_ANARCHY_DU_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file AnarchyDU/hydro_parameters.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (default compile-time
+ * parameters).
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/*! Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/*! Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/*! The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 2.0f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/*! The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.1f
+
+/*! Minimal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_min 0.0f
+
+/*! Maximal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_max 2.0f
+
+/*! Decay length for the viscosity scheme. This is scheme dependent. In
+ * non-variable schemes this must be defined but is not used. This also
+ * sets the decay length for the diffusion. */
+#define hydro_props_default_viscosity_length 0.25f
+
+/* Diffusion parameters -- Defaults; can be changed at run-time */
+
+/*! The "initial" diffusion, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.0. */
+#define hydro_props_default_diffusion_alpha 0.0f
+
+/*! Beta coefficient for the diffusion. This controls how fast the
+ * diffusion coefficient peaks, and how high it can get. Chosen to be
+ * very small in schemes where little diffusion is needed, 0.2-1.0 in
+ * schemes (e.g. density-energy) where diffusion is needed to solve
+ * the contact discontinuity problem. */
+#define hydro_props_default_diffusion_beta 0.25f
+
+/*! Maximal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_max 1.0f
+
+/*! Minimal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_min 0.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+
+ /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+ float alpha_max;
+
+ /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+ float alpha_min;
+
+ /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+ float length;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {
+
+ /*! Initialisation value, or the case for constant thermal diffusion coeffs
+ */
+ float alpha;
+
+ /*! Tuning parameter for speed of ramp up/down */
+ float beta;
+
+ /*! Maximal value for alpha_diff */
+ float alpha_max;
+
+ /*! Minimal value for alpha_diff */
+ float alpha_min;
+};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+
+ viscosity->alpha_max =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
+ hydro_props_default_viscosity_alpha_max);
+
+ viscosity->alpha_min =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
+ hydro_props_default_viscosity_alpha_min);
+
+ viscosity->length = parser_get_opt_param_float(
+ params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+ viscosity->alpha_max = hydro_props_default_viscosity_alpha_max;
+ viscosity->alpha_min = hydro_props_default_viscosity_alpha_min;
+ viscosity->length = hydro_props_default_viscosity_length;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message(
+ "Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
+ "min: %.3f, length: %.3f.",
+ viscosity->alpha, viscosity->alpha_max, viscosity->alpha_min,
+ viscosity->length);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (max)", viscosity->alpha_max);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (min)", viscosity->alpha_min);
+ io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
+ viscosity->length);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {
+
+ diffusion->alpha = parser_get_opt_param_float(
+ params, "SPH:diffusion_alpha", hydro_props_default_diffusion_alpha);
+
+ diffusion->beta = parser_get_opt_param_float(
+ params, "SPH:diffusion_beta", hydro_props_default_diffusion_beta);
+
+ diffusion->alpha_max =
+ parser_get_opt_param_float(params, "SPH:diffusion_alpha_max",
+ hydro_props_default_diffusion_alpha_max);
+
+ diffusion->alpha_min =
+ parser_get_opt_param_float(params, "SPH:diffusion_alpha_min",
+ hydro_props_default_diffusion_alpha_min);
+}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {
+ diffusion->alpha = hydro_props_default_diffusion_alpha;
+ diffusion->alpha_max = hydro_props_default_diffusion_alpha_max;
+ diffusion->alpha_min = hydro_props_default_diffusion_alpha_min;
+ diffusion->beta = hydro_props_default_diffusion_beta;
+}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {
+ message(
+ "Artificial diffusion parameters set to alpha: %.3f, max: %.3f, "
+ "min: %.3f, beta: %.3f.",
+ diffusion->alpha, diffusion->alpha_max, diffusion->alpha_min,
+ diffusion->beta);
+}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {
+ io_write_attribute_f(h_grpsph, "Diffusion alpha", diffusion->alpha);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (max)", diffusion->alpha_max);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (min)", diffusion->alpha_min);
+ io_write_attribute_f(h_grpsph, "Diffusion beta", diffusion->beta);
+}
+#endif
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/AnarchyDU/hydro_part.h b/src/hydro/AnarchyDU/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..98bb058af9881fd050e6b9f2101766870ed26a46
--- /dev/null
+++ b/src/hydro/AnarchyDU/hydro_part.h
@@ -0,0 +1,211 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) &
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ANARCHY_DU_HYDRO_PART_H
+#define SWIFT_ANARCHY_DU_HYDRO_PART_H
+
+/**
+ * @file AnarchyDU/hydro_part.h
+ * @brief Density-Energy conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (particle definition)
+ */
+
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+#include "star_formation_struct.h"
+#include "tracers_struct.h"
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Offset between the current position and position at the last sort. */
+ float x_diff_sort[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
+
+ /*! Gravitational acceleration at the last full step. */
+ float a_grav[3];
+
+ /*! Internal energy at the last full step. */
+ float u_full;
+
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+ /* Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
+ /* Additional data used by the tracers */
+ struct star_formation_xpart_data sf_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+ /*! Particle unique ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle predicted velocity. */
+ float v[3];
+
+ /*! Particle acceleration. */
+ float a_hydro[3];
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle smoothing length. */
+ float h;
+
+ /*! Particle internal energy. */
+ float u;
+
+ /*! Time derivative of the internal energy. */
+ float u_dt;
+
+ /*! Particle density. */
+ float rho;
+
+ /* Store viscosity information in a separate struct. */
+ struct {
+
+ /*! Particle velocity divergence */
+ float div_v;
+
+ /*! Particle velocity divergence from previous step */
+ float div_v_previous_step;
+
+ /*! Artificial viscosity parameter */
+ float alpha;
+
+ /*! Signal velocity */
+ float v_sig;
+
+ } viscosity;
+
+ /* Store thermal diffusion information in a separate struct. */
+ struct {
+
+ /*! del^2 u, a smoothed quantity */
+ float laplace_u;
+
+ /*! Thermal diffusion coefficient */
+ float alpha;
+
+ } diffusion;
+
+ /* Store density/force specific stuff. */
+ union {
+
+ /**
+ * @brief Structure for the variables only used in the density loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the density
+ * loop over neighbours and the ghost task.
+ */
+ struct {
+
+ /*! Neighbour number count. */
+ float wcount;
+
+ /*! Derivative of the neighbour number with respect to h. */
+ float wcount_dh;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ /*! Particle velocity curl. */
+ float rot_v[3];
+
+ } density;
+
+ /**
+ * @brief Structure for the variables only used in the force loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the force
+ * loop over neighbours and the ghost, drift and kick tasks.
+ */
+ struct {
+
+ /*! "Grad h" term -- only partial in P-U */
+ float f;
+
+ /*! Particle pressure. */
+ float pressure;
+
+ /*! Particle soundspeed. */
+ float soundspeed;
+
+ /*! Time derivative of smoothing length */
+ float h_dt;
+
+ /*! Balsara switch */
+ float balsara;
+
+ } force;
+ };
+
+ /* Chemistry information */
+ struct chemistry_part_data chemistry_data;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+ /* Need waking up ? */
+ timebin_t wakeup;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_ANARCHY_DU_HYDRO_PART_H */
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 28445ad3ed4745e143919eb800966a7a06697b27..460461d06af954f9058e254fcdf2e9a8ced21d14 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -21,14 +21,11 @@
#define SWIFT_ANARCHY_PU_HYDRO_H
/**
- * @file PressureEnergy/hydro.h
- * @brief P-U conservative implementation of SPH (Non-neighbour loop
- * equations)
- *
- * The thermal variable is the internal energy (u). A simple constant
- * viscosity term with a Balsara switch is implemented.
- *
- * No thermal conduction term is implemented.
+ * @file AnarchyPU/hydro.h
+ * @brief P-U conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (Non-neighbour loop
+ * equations)
*
* This implementation corresponds to the one presented in the SWIFT
* documentation and in Hopkins, "A general class of Lagrangian smoothed
@@ -47,6 +44,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
#include <float.h>
/**
@@ -420,6 +419,29 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ p->viscosity.alpha = alpha;
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ hydro_set_viscosity_alpha(p,
+ hydro_props_default_viscosity_alpha_feedback_reset);
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
diff --git a/src/hydro/AnarchyPU/hydro_debug.h b/src/hydro/AnarchyPU/hydro_debug.h
index 79ab5b96653a3f503c1baf255f4296f0ccc4aca9..f52a63cd31c608b9f1f0e998a56224e9e133f157 100644
--- a/src/hydro/AnarchyPU/hydro_debug.h
+++ b/src/hydro/AnarchyPU/hydro_debug.h
@@ -20,8 +20,10 @@
#ifndef SWIFT_ANARCHY_PU_HYDRO_DEBUG_H
#define SWIFT_ANARCHY_PU_HYDRO_DEBUG_H
/**
- * @file PressureEnergy/hydro_debug.h
- * @brief P-U conservative implementation of SPH (Debugging routines)
+ * @file AnarchyPU/hydro_debug.h
+ * @brief P-U conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (Debugging routines).
*/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h
index f3a8734fbf037b43d509750e018d51c9b01933b4..bf683fb447a3c6aa10ca750e5c967356f9dabbfd 100644
--- a/src/hydro/AnarchyPU/hydro_iact.h
+++ b/src/hydro/AnarchyPU/hydro_iact.h
@@ -21,20 +21,19 @@
#define SWIFT_ANARCHY_PU_HYDRO_IACT_H
/**
- * @file PressureEnergy/hydro_iact.h
- * @brief P-U implementation of SPH (Neighbour loop equations)
+ * @file AnarchyPU/hydro_iact.h
+ * @brief P-U conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion) (Neighbour loop equations).
*
- * The thermal variable is the internal energy (u). A simple constant
- * viscosity term with a Balsara switch is implemented.
- *
- * No thermal conduction term is implemented.
- *
- * See PressureEnergy/hydro.h for references.
+ * See AnarchyPU/hydro.h for references.
*/
#include "adiabatic_index.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two part*icles.
*
diff --git a/src/hydro/AnarchyPU/hydro_io.h b/src/hydro/AnarchyPU/hydro_io.h
index b2a411cf408a2067ea5490869be004f9fd26954d..e1525cc99db8074a99c8d28a73918adb7b7b5319 100644
--- a/src/hydro/AnarchyPU/hydro_io.h
+++ b/src/hydro/AnarchyPU/hydro_io.h
@@ -20,15 +20,12 @@
#ifndef SWIFT_ANARCHY_PU_HYDRO_IO_H
#define SWIFT_ANARCHY_PU_HYDRO_IO_H
/**
- * @file PressureEnergy/hydro_io.h
- * @brief P-U implementation of SPH (i/o routines)
+ * @file AnarchyPU/hydro_io.h
+ * @brief P-U conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion) (i/o routines)
*
- * The thermal variable is the internal energy (u). A simple constant
- * viscosity term with a Balsara switch is implemented.
- *
- * No thermal conduction term is implemented.
- *
- * See PressureEnergy/hydro.h for references.
+ * See AnarchyPU/hydro.h for references.
*/
#include "adiabatic_index.h"
@@ -36,6 +33,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
@@ -205,7 +204,8 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
/* Nothing in this minimal model... */
- io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+ io_write_attribute_s(h_grpsph, "Thermal Conductivity Model",
+ "Simple treatment as in Price (2008)");
io_write_attribute_s(h_grpsph, "Viscosity Model",
"Simplified version of Cullen & Denhen (2011)");
diff --git a/src/hydro/AnarchyPU/hydro_parameters.h b/src/hydro/AnarchyPU/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..6b1114a2677091b435577bf06d78380983c70651
--- /dev/null
+++ b/src/hydro/AnarchyPU/hydro_parameters.h
@@ -0,0 +1,298 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_ANARCHY_PU_HYDRO_PARAMETERS_H
+#define SWIFT_ANARCHY_PU_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file AnarchyPU/hydro_parameters.h
+ * @brief P-PU conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion (default compile-time
+ * parameters).
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 2.0f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.1f
+
+/* Minimal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_min 0.0f
+
+/* Maximal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_max 2.0f
+
+/* Decay length for the viscosity scheme. This is scheme dependent. In
+ * non-variable schemes this must be defined but is not used. This also
+ * sets the decay length for the diffusion. */
+#define hydro_props_default_viscosity_length 0.25f
+
+/* Diffusion parameters -- Defaults; can be changed at run-time */
+
+/* The "initial" diffusion, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.0. */
+#define hydro_props_default_diffusion_alpha 0.0f
+
+/* Beta coefficient for the diffusion. This controls how fast the
+ * diffusion coefficient peaks, and how high it can get. Chosen to be
+ * very small in schemes where little diffusion is needed, 0.2-1.0 in
+ * schemes (e.g. density-energy) where diffusion is needed to solve
+ * the contact discontinuity problem. */
+#define hydro_props_default_diffusion_beta 0.01f
+
+/* Maximal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_max 1.0f
+
+/* Minimal value for the diffusion alpha in variable schemes. */
+#define hydro_props_default_diffusion_alpha_min 0.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+
+ /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+ float alpha_max;
+
+ /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+ float alpha_min;
+
+ /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+ float length;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {
+
+ /*! Initialisation value, or the case for constant thermal diffusion coeffs
+ */
+ float alpha;
+
+ /*! Tuning parameter for speed of ramp up/down */
+ float beta;
+
+ /*! Maximal value for alpha_diff */
+ float alpha_max;
+
+ /*! Minimal value for alpha_diff */
+ float alpha_min;
+};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+
+ viscosity->alpha_max =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
+ hydro_props_default_viscosity_alpha_max);
+
+ viscosity->alpha_min =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
+ hydro_props_default_viscosity_alpha_min);
+
+ viscosity->length = parser_get_opt_param_float(
+ params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+ viscosity->alpha_max = hydro_props_default_viscosity_alpha_max;
+ viscosity->alpha_min = hydro_props_default_viscosity_alpha_min;
+ viscosity->length = hydro_props_default_viscosity_length;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message(
+ "Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
+ "min: %.3f, length: %.3f.",
+ viscosity->alpha, viscosity->alpha_max, viscosity->alpha_min,
+ viscosity->length);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (max)", viscosity->alpha_max);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (min)", viscosity->alpha_min);
+ io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
+ viscosity->length);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {
+
+ diffusion->alpha = parser_get_opt_param_float(
+ params, "SPH:diffusion_alpha", hydro_props_default_diffusion_alpha);
+
+ diffusion->beta = parser_get_opt_param_float(
+ params, "SPH:diffusion_beta", hydro_props_default_diffusion_beta);
+
+ diffusion->alpha_max =
+ parser_get_opt_param_float(params, "SPH:diffusion_alpha_max",
+ hydro_props_default_diffusion_alpha_max);
+
+ diffusion->alpha_min =
+ parser_get_opt_param_float(params, "SPH:diffusion_alpha_min",
+ hydro_props_default_diffusion_alpha_min);
+}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {
+ diffusion->alpha = hydro_props_default_diffusion_alpha;
+ diffusion->alpha_max = hydro_props_default_diffusion_alpha_max;
+ diffusion->alpha_min = hydro_props_default_diffusion_alpha_min;
+ diffusion->beta = hydro_props_default_diffusion_beta;
+}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {
+ message(
+ "Artificial diffusion parameters set to alpha: %.3f, max: %.3f, "
+ "min: %.3f, beta: %.3f.",
+ diffusion->alpha, diffusion->alpha_max, diffusion->alpha_min,
+ diffusion->beta);
+}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {
+ io_write_attribute_f(h_grpsph, "Diffusion alpha", diffusion->alpha);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (max)", diffusion->alpha_max);
+ io_write_attribute_f(h_grpsph, "Diffusion alpha (min)", diffusion->alpha_min);
+ io_write_attribute_f(h_grpsph, "Diffusion beta", diffusion->beta);
+}
+#endif
+
+#endif /* SWIFT_ANARCHY_PU_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/AnarchyPU/hydro_part.h b/src/hydro/AnarchyPU/hydro_part.h
index 2c5022c262587c50f577970e1d1891a42b70491b..3c59f2d7c1b8610eeef816c674038dcab5554130 100644
--- a/src/hydro/AnarchyPU/hydro_part.h
+++ b/src/hydro/AnarchyPU/hydro_part.h
@@ -20,15 +20,12 @@
#ifndef SWIFT_ANARCHY_PU_HYDRO_PART_H
#define SWIFT_ANARCHY_PU_HYDRO_PART_H
/**
- * @file PressureEnergy/hydro_part.h
- * @brief P-U implementation of SPH (Particle definition)
+ * @file AnarchyPU/hydro_part.h
+ * @brief P-U conservative implementation of SPH,
+ * with added ANARCHY physics (Cullen & Denhen 2011 AV,
+ * Price 2008 thermal diffusion) (Particle definition)
*
- * The thermal variable is the internal energy (u). A simple constant
- * viscosity term with a Balsara switch is implemented.
- *
- * No thermal conduction term is implemented.
- *
- * See PressureEnergy/hydro.h for references.
+ * See AnarchyPU/hydro.h for references.
*/
#include "chemistry_struct.h"
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 635c5a7df1cef2cf21d61dc4022eddbd0e4b63a4..4bcc0bcd4fa4ebaf8f43d738230c597f19873ce5 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -30,6 +30,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
#include <float.h>
/**
@@ -375,6 +377,29 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ p->alpha = alpha;
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ hydro_set_viscosity_alpha(p,
+ hydro_props_default_viscosity_alpha_feedback_reset);
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -554,6 +579,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Set the AV property */
p->alpha = hydro_props->viscosity.alpha;
+ /* Set the diffusion parameter */
+ p->alpha_diff = hydro_props->diffusion.alpha;
+
/* Viscosity parameter decay time */
/* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
*/
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 85c586a4e921e38296453b71a2a2b9637971c28c..8ce77b9e353bf38b31dbd7b7a5444c10f3f6b5bc 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -22,6 +22,8 @@
#include "adiabatic_index.h"
+#include "./hydro_parameters.h"
+
/**
* @brief SPH interaction functions following the Gadget-2 version of SPH.
*
@@ -241,7 +243,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Thermal conductivity */
v_sig_u = sqrtf(2.f * hydro_gamma_minus_one *
fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
- tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
+ tc = pi->alpha_diff * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr);
/* Get the common factor out. */
@@ -351,7 +353,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Thermal conductivity */
v_sig_u = sqrtf(2.f * hydro_gamma_minus_one *
fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
- tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
+ tc = pi->alpha_diff * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr);
/* Get the common factor out. */
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 69919c202223fdecc197a87178e59767c02ee16e..858e7e02d8b2ab2a9458f74bb949a5a14faca2cf 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -24,6 +24,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
@@ -177,8 +179,7 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
io_write_attribute_s(h_grpsph, "Thermal Conductivity Model",
"Price (2008) without switch");
- io_write_attribute_f(h_grpsph, "Thermal Conductivity alpha",
- const_conductivity_alpha);
+
io_write_attribute_s(
h_grpsph, "Viscosity Model",
"Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
diff --git a/src/hydro/Default/hydro_parameters.h b/src/hydro/Default/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..f11a2466cd911dbde960aa91e3bb50bc3d51f1f7
--- /dev/null
+++ b/src/hydro/Default/hydro_parameters.h
@@ -0,0 +1,212 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_DEFAULT_HYDRO_PARAMETERS_H
+#define SWIFT_DEFAULT_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file Default/hydro_parameters.h
+ * @brief Default implementation of SPH (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Diffusion parameters -- Defaults; can be changed at run-time */
+
+/* The "initial" diffusion, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.0. */
+#define hydro_props_default_diffusion_alpha 1.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {
+
+ /*! Initialisation value, or the case for constant thermal diffusion coeffs
+ */
+ float alpha;
+};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {
+
+ diffusion->alpha = parser_get_opt_param_float(
+ params, "SPH:diffusion_alpha", hydro_props_default_diffusion_alpha);
+}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {
+ diffusion->alpha = hydro_props_default_diffusion_alpha;
+}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {
+ message("Artificial diffusion parameters set to alpha: %.3f",
+ diffusion->alpha);
+}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {
+ io_write_attribute_f(h_grpsph, "Diffusion alpha", diffusion->alpha);
+}
+#endif
+
+#endif /* SWIFT_DEFAULT_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 2ff41da0a7825d2ddba276c3d20a9d2d4aad3224..c262a2db71048858231b9886977b2b8a8f98cf00 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -89,6 +89,9 @@ struct part {
/* Particle viscosity parameter */
float alpha;
+ /* Particle diffusion parameter */
+ float alpha_diff;
+
/* Store density/force specific stuff. */
union {
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 20c65576cf33488db2a942e26fa6be94da5fb350..c6367ec7b46076d97f5637d360d946dc056d7bd4 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -42,6 +42,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Returns the comoving internal energy of a particle at the last
* time the particle was kicked.
@@ -399,6 +401,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 819ba3889b008e859939ff1d0c2b6ac62a1be30a..118a29d771facbafbb15ecef90d3d3c726e1b92c 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -35,6 +35,8 @@
#include "cache.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two particles.
*
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index db6a899f05b3bf5dd8cbbe1db01ca2c3ae3d14ca..5d68dc5d9f0f4ac0d12f7829f2c2a12d9964fe4a 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -24,6 +24,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/Gadget2/hydro_parameters.h b/src/hydro/Gadget2/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f00554bdfa78fe93e82d6476bd0a2d9d9983614
--- /dev/null
+++ b/src/hydro/Gadget2/hydro_parameters.h
@@ -0,0 +1,189 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GADGET2_HYDRO_PARAMETERS_H
+#define SWIFT_GADGET2_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file Gadget2/hydro_parameters.h
+ * @brief Follows the Gadget-2 version of SPH. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param us: pointer to the internal unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param us: pointer to the internal unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_GADGET2_HYDRO_PARAMETERS_H */
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index c744a189a4a41720d7f0ee1f856fafd43e469033..8c47d529dd4696331b9afbb3e402639116bd06f1 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -652,7 +652,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Apply the minimal energy limit */
const float min_energy =
- hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy * p->conserved.mass) {
p->conserved.energy = min_energy * p->conserved.mass;
p->flux.energy = 0.0f;
@@ -1030,6 +1030,28 @@ hydro_set_drifted_physical_internal_energy(struct part* p,
error("Need implementing");
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part* restrict p, float alpha) {
+ /* Purposefully left empty */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part* restrict p) {
+ /* Purposefully left empty */
+}
+
/**
* @brief Returns the comoving density of a particle
*
@@ -1117,15 +1139,10 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
__attribute__((always_inline)) INLINE static void
hydro_set_init_internal_energy(struct part* p, float u_init) {
- p->conserved.energy = u_init * p->conserved.mass;
-#ifdef GIZMO_TOTAL_ENERGY
- /* add the kinetic energy */
- p->conserved.energy +=
- 0.5f * p->conserved.mass *
- (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
- p->conserved.momentum[2] * p->v[2]);
-#endif
- p->P = hydro_gamma_minus_one * p->rho * u_init;
+ /* We store the initial energy per unit mass in the energy
+ * variable as the conversion to energy will be done later,
+ * in hydro_first_init_part(). */
+ p->conserved.energy = u_init;
}
/**
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index 09d4c7c70ee2bae8a31d10cb4a568c4627c7b3cd..08e3de538a5e2cd775a69abab9bcb11fb09e7a9b 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -25,6 +25,8 @@
#include "hydro_gradients.h"
#include "riemann.h"
+#include "./hydro_parameters.h"
+
#define GIZMO_VOLUME_CORRECTION
/**
diff --git a/src/hydro/GizmoMFM/hydro_parameters.h b/src/hydro/GizmoMFM/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..5c345d6055be330ff7c8145c88edecc609bb9009
--- /dev/null
+++ b/src/hydro/GizmoMFM/hydro_parameters.h
@@ -0,0 +1,162 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_GIZMOMFM_HYDRO_PARAMETERS_H
+#define SWIFT_GIZMOMFM_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file GizmoMFM/hydro_parameters.h
+ * @brief Gizmo-MFM scheme. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ *
+ * These aren't really used in GIZMO but must be defined.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_GIZMOMFM_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index 1f26f29fa5f3f9d496c73116c7608f788f93d605..dc11bcece94bbd0622f4eb48ca2b14d1fcd9d235 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -744,7 +744,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Apply the minimal energy limit */
const float min_energy =
- hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy * p->conserved.mass) {
p->conserved.energy = min_energy * p->conserved.mass;
p->conserved.flux.energy = 0.f;
@@ -1101,6 +1101,28 @@ hydro_set_drifted_physical_internal_energy(struct part* p,
error("Need implementing");
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part* restrict p, float alpha) {
+ /* Purposefully left empty */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part* restrict p) {
+ /* Purposefully left empty */
+}
+
/**
* @brief Returns the comoving density of a particle
*
@@ -1188,15 +1210,10 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
__attribute__((always_inline)) INLINE static void
hydro_set_init_internal_energy(struct part* p, float u_init) {
- p->conserved.energy = u_init * p->conserved.mass;
-#ifdef GIZMO_TOTAL_ENERGY
- /* add the kinetic energy */
- p->conserved.energy += 0.5f * p->conserved.mass *
- (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
-#endif
- p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+ /* We store the initial energy per unit mass in the energy
+ * variable as the conversion to energy will be done later,
+ * in hydro_first_init_part(). */
+ p->conserved.energy = u_init;
}
/**
diff --git a/src/hydro/GizmoMFV/hydro_iact.h b/src/hydro/GizmoMFV/hydro_iact.h
index d882549f8c55018419a2e1730d2ac099bbe1f5ee..2ac5474f0e767879656129f486fb276608d13622 100644
--- a/src/hydro/GizmoMFV/hydro_iact.h
+++ b/src/hydro/GizmoMFV/hydro_iact.h
@@ -25,6 +25,8 @@
#include "hydro_gradients.h"
#include "riemann.h"
+#include "./hydro_parameters.h"
+
#define GIZMO_VOLUME_CORRECTION
/**
diff --git a/src/hydro/GizmoMFV/hydro_parameters.h b/src/hydro/GizmoMFV/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..29b618444e4fa48dcddfa3d688a5dffc07e0433d
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_parameters.h
@@ -0,0 +1,162 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ *****************************************************************************/
+
+#ifndef SWIFT_GIZMOMFV_HYDRO_PARAMETERS_H
+#define SWIFT_GIZMOMFV_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file GizmoMFV/hydro_parameters.h
+ * @brief Gizmo-MFV scheme. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ *
+ * These aren't really used in GIZMO but must be defined.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_GIZMOMFV_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 6eccaf26999d5f2a6dbe29414136919d7a5b7e14..33e01a783c4bf854707d7cf72a67c437e4f3ba0e 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -44,6 +44,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Returns the comoving internal energy of a particle at the last
* time the particle was kicked.
@@ -394,6 +396,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 7fc7a3c67f6c832d70109319ad964e25df30ff4e..8929ede0ac8ef736776e745323d33bb0b2265021 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -35,6 +35,8 @@
#include "adiabatic_index.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two particles.
*
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 1146aa9347d443833cd481103da6f6c57d21fcbf..c6e36e32d6176c6968e70c7ba689b6651a2d1c18 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -38,6 +38,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/Minimal/hydro_parameters.h b/src/hydro/Minimal/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..6aa483847aa6071c3ed2928b488511004c30b597
--- /dev/null
+++ b/src/hydro/Minimal/hydro_parameters.h
@@ -0,0 +1,190 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_MINIMAL_HYDRO_PARAMETERS_H
+#define SWIFT_MINIMAL_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file Minimal/hydro_parameters.h
+ * @brief Minimal conservative implementation of SPH . (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param us: pointer to the internal unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param us: pointer to the internal unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_MINIMAL_HYDRO_PARAMETERS_H */
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index ba9228289fd43e6cf130b58e968de63e4fee3147..7a6d0861cbe0b22f8e687f9781db0a4585c8b569 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -46,6 +46,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/*
* Note: Define PLANETARY_SPH_NO_BALSARA to disable the Balsara (1995) switch
* for the artificial viscosity and use the vanilla Monaghan (1992) instead.
@@ -418,6 +420,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index afebb6a406bd310f38d51dcb32fc25da6b2674b5..3a8d1f4e88fa9df0a1dcea500baf56e46bae7cd1 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -37,6 +37,8 @@
#include "const.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two particles.
*
diff --git a/src/hydro/Planetary/hydro_io.h b/src/hydro/Planetary/hydro_io.h
index 1b84f8d6db295694846ffd26a422ce158aad0c60..64f229be3087fcb225b2b934b613f2eda0d84eba 100644
--- a/src/hydro/Planetary/hydro_io.h
+++ b/src/hydro/Planetary/hydro_io.h
@@ -39,6 +39,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/Planetary/hydro_parameters.h b/src/hydro/Planetary/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..3e15b05b9418bc1ebd558a03e6c87887e5b52c7c
--- /dev/null
+++ b/src/hydro/Planetary/hydro_parameters.h
@@ -0,0 +1,190 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_PLANETARY_HYDRO_PARAMETERS_H
+#define SWIFT_PLANETARY_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file Planetary/hydro_parameters.h
+ * @brief Minimal conservative implementation of SPH . (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Planetary default beta=4.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 4.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 1.5f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 1.5f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_PLANETARY_HYDRO_PARAMETERS_H */
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 2ad1d70e8016e3cf1837acace3c6df20efd152f2..1691781a12e3c4c2299a08a4aac08bc83c1d123d 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -47,6 +47,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
#include <float.h>
/**
@@ -419,6 +421,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -886,22 +910,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
__attribute__((always_inline)) INLINE static void hydro_remove_part(
const struct part *p, const struct xpart *xp) {}
-/**
- * @brief Inputs extra heat to a particle at redshift of reionization
- *
- * We assume a constant density.
- * @param p The particle of interest
- * @param xp The extended particle data
- * @param cosmo Cosmology data structure
- * @param extra_heat The extra internal energy given to the particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_reion_heating(
- struct part *p, struct xpart *xp, const struct cosmology *cosmo,
- float extra_heat) {
-
- const float old_u = xp->u_full * cosmo->a_factor_internal_energy;
- const float new_u = old_u + extra_heat;
- xp->u_full = new_u / cosmo->a_factor_internal_energy;
-}
-
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_H */
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index ae154ea549a52cb24ed7c69453533b7d59b39a85..e2fc6e7e92423397278f95bd30fb1b3d5d6c93a4 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -35,6 +35,8 @@
#include "adiabatic_index.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two part*icles.
*
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 701c12283bf77acef4af77598f57705a2b364fa1..3e645803ebc3cbdbc361e73d776824e2ea7906fa 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -36,6 +36,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/PressureEnergy/hydro_parameters.h b/src/hydro/PressureEnergy/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fccc58a17545e862cc618ce79d37e4addc3a6de
--- /dev/null
+++ b/src/hydro/PressureEnergy/hydro_parameters.h
@@ -0,0 +1,190 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_PARAMETERS_H
+#define SWIFT_PRESSURE_ENERGY_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file PressureEnergy/hydro_parameters.h
+ * @brief P-U implementation of SPH. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index 1a93cb98ebc806131434b442b8e82f42a3ab1578..e477c9f96b2daaf1d786dfc58070da8fd89a7f88 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -21,7 +21,7 @@
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_H
/**
- * @file PressureEnergy/hydro.h
+ * @file PressureEnergyMorrisMonaghanAV/hydro.h
* @brief P-U conservative implementation of SPH (Non-neighbour loop
* equations)
*
@@ -48,6 +48,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
#include <float.h>
/**
@@ -416,6 +418,29 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
p->force.soundspeed = soundspeed;
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ p->alpha = alpha;
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ hydro_set_viscosity_alpha(p,
+ hydro_props_default_viscosity_alpha_feedback_reset);
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
index d0cd5367f94cd90f36cc2b738a63c7963adbd445..7f728e9a5d390a8508730725ff2f8eb82126f616 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
@@ -20,7 +20,7 @@
#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
/**
- * @file PressureEnergy/hydro_debug.h
+ * @file PressureEnergyMorrisMonaghanAV/hydro_debug.h
* @brief P-U conservative implementation of SPH (Debugging routines)
*
* The thermal variable is the internal energy (u). A simple variable
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
index 69da511c7544a71ef381a0889c8b56c80d5211f1..d7521d3387afc97666e88780aa7894e3a1459b4a 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
@@ -21,7 +21,7 @@
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H
/**
- * @file PressureEnergy/hydro_iact.h
+ * @file PressureEnergyMorrisMonaghanAV/hydro_iact.h
* @brief P-U implementation of SPH (Neighbour loop equations)
*
* The thermal variable is the internal energy (u). A simple variable
@@ -36,6 +36,8 @@
#include "adiabatic_index.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Density interaction between two part*icles.
*
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
index 71662f14c61c92d65bcf493b6f5a43b8172e3697..fe44864b2c360373a3eff676775de6e9fac96266 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
@@ -20,7 +20,7 @@
#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IO_H
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IO_H
/**
- * @file PressureEnergy/hydro_io.h
+ * @file PressureEnergyMorrisMonaghanAV/hydro_io.h
* @brief P-U implementation of SPH (i/o routines)
*
* The thermal variable is the internal energy (u). A simple variable
@@ -37,6 +37,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..0c0ae33c29e2f13331a002df3aa2d0fa1162b49d
--- /dev/null
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h
@@ -0,0 +1,231 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PARAMETERS_H
+#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file PressureEnergyMorrisMonaghanAV/hydro_parameters.h
+ * @brief P-U implementation of SPH. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 2.0f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Minimal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_min 0.1f
+
+/* Maximal value for the viscosity alpha in variable schemes. */
+#define hydro_props_default_viscosity_alpha_max 2.0f
+
+/* Decay length for the viscosity scheme. This is scheme dependent. In
+ * non-variable schemes this must be defined but is not used. This also
+ * sets the decay length for the diffusion. */
+#define hydro_props_default_viscosity_length 0.1f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+
+ /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+ float alpha_max;
+
+ /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+ float alpha_min;
+
+ /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+ float length;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+
+ viscosity->alpha_max =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
+ hydro_props_default_viscosity_alpha_max);
+
+ viscosity->alpha_min =
+ parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
+ hydro_props_default_viscosity_alpha_min);
+
+ viscosity->length = parser_get_opt_param_float(
+ params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+ viscosity->alpha_max = hydro_props_default_viscosity_alpha_max;
+ viscosity->alpha_min = hydro_props_default_viscosity_alpha_min;
+ viscosity->length = hydro_props_default_viscosity_length;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message(
+ "Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
+ "min: %.3f, length: %.3f.",
+ viscosity->alpha, viscosity->alpha_max, viscosity->alpha_min,
+ viscosity->length);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (max)", viscosity->alpha_max);
+ io_write_attribute_f(h_grpsph, "Alpha viscosity (min)", viscosity->alpha_min);
+ io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
+ viscosity->length);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
index ecd20938456b04004ed2299fbe1de0c1b8bb50d6..5a520f522a85ac87376d7db8cc78f58c66e3550b 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
@@ -20,7 +20,7 @@
#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PART_H
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_PART_H
/**
- * @file PressureEnergy/hydro_part.h
+ * @file PressureEnergyMorrisMonaghanAV/hydro_part.h
* @brief P-U implementation of SPH (Particle definition)
*
* The thermal variable is the internal energy (u). A simple variable
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 8fdd201cf29ce43e5c2207dd0b219d08e9e34515..a381011069c9a5f0717f1a26ffef1c4bd7eeb817 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -42,6 +42,8 @@
#include "kernel_hydro.h"
#include "minmax.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Returns the comoving internal energy of a particle at the last
* time the particle was kicked.
@@ -181,6 +183,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
error("To be implemented");
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
/**
* @brief Returns the comoving entropy of a particle drifted to the
* current time.
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index 19279adec1f37117cf985e63a18a681ceee4f973..d5045d77f2850408b33a70e9b9d6d880577bea0f 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -19,6 +19,8 @@
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H
+#include "./hydro_parameters.h"
+
/**
* @file PressureEntropy/hydro_iact.h
* @brief Pressure-Entropy implementation of SPH (Particle interactions)
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index e9397bf6108b8bc16658157e424055274f05f23c..5c0bb71d3dcfe77619d18115e27fbafbac3facec 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -35,6 +35,8 @@
#include "io_properties.h"
#include "kernel_hydro.h"
+#include "./hydro_parameters.h"
+
/**
* @brief Specifies which particle fields to read from a dataset
*
diff --git a/src/hydro/PressureEntropy/hydro_parameters.h b/src/hydro/PressureEntropy/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..40fa6efd74648a328f6304aaf28340d17242d0c1
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro_parameters.h
@@ -0,0 +1,190 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_PARAMETERS_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file PressureEnergy/hydro_parameters.h
+ * @brief P-A implementation of SPH. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * hydro_props_default_viscosity_alpha in fixed schemes, and likely
+ * to hydro_props_default_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 0.8f
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/* The "initial" hydro viscosity, or the fixed value for non-variable
+ * schemes. This usually takes the value 0.8. */
+#define hydro_props_default_viscosity_alpha 0.8f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {
+ /*! For the fixed, simple case. Also used to set the initial AV
+ coefficient for variable schemes. */
+ float alpha;
+};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {
+
+ /* Read the artificial viscosity parameters from the file, if they exist,
+ * otherwise set them to the defaults defined above. */
+
+ viscosity->alpha = parser_get_opt_param_float(
+ params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
+}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {
+ viscosity->alpha = hydro_props_default_viscosity_alpha;
+}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {
+ message("Artificial viscosity parameters set to alpha: %.3f",
+ viscosity->alpha);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Alpha viscosity", viscosity->alpha);
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index ed1b000199f3641ff458c05e6e5268a8b183c468..a6b3b41c9c08e6f83ed4d989bd57176507221797 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -874,6 +874,28 @@ hydro_set_drifted_physical_internal_energy(struct part* p,
error("Need implementing");
}
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part* restrict p, float alpha) {
+ /* Purposefully left empty */
+}
+
+/**
+ * @brief Update the value of the viscosity alpha to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_viscosity_alpha_max_feedback(struct part* restrict p) {
+ /* Purposefully left empty */
+}
+
/**
* @brief Returns the comoving pressure of a particle
*
diff --git a/src/hydro/Shadowswift/hydro_parameters.h b/src/hydro/Shadowswift/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..93cb3bb7367b7c0f58f9c130e0e159dd174f0e6a
--- /dev/null
+++ b/src/hydro/Shadowswift/hydro_parameters.h
@@ -0,0 +1,160 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_SHADOWSWIFT_HYDRO_PARAMETERS_H
+#define SWIFT_SHADOWSWIFT_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "../../../config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file Shadowswift/hydro_parameters.h
+ * @brief Shadowswift implementation. (default parameters)
+ *
+ * This file defines a number of things that are used in
+ * hydro_properties.c as defaults for run-time parameters
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity parameters -- FIXED -- MUST BE DEFINED AT COMPILE-TIME */
+
+/* Cosmology default beta=3.0.
+ * Alpha can be set in the parameter file.
+ * Beta is defined as in e.g. Price (2010) Eqn (103) */
+#define const_viscosity_beta 3.0f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {};
+
+/*! Thermal diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+ io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {}
+#endif
+
+#endif /* SWIFT_SHADOWSWIFT_HYDRO_PARAMETERS_H */
\ No newline at end of file
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 60e5593cc0630ef4bc33ab407f6a669b7de1def1..ccfb57371ac28ba4931fc14d547b0b674c77d48d 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -43,6 +43,8 @@
#include "./hydro/Shadowswift/hydro_io.h"
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_io.h"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro_io.h"
#elif defined(ANARCHY_PU_SPH)
#include "./hydro/AnarchyPU/hydro_io.h"
#else
diff --git a/src/hydro_parameters.h b/src/hydro_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..322bea5496893965fc727f957f54a19e57ee3e44
--- /dev/null
+++ b/src/hydro_parameters.h
@@ -0,0 +1,60 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_HYDRO_PARAMETERS_H
+#define SWIFT_HYDRO_PARAMETERS_H
+
+/**
+ * @file src/hydro_parameters.h
+ * @brief Contains all the parameters of the hydro schemes, included from
+ * their own local header.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Import the right hydro header */
+#if defined(MINIMAL_SPH)
+#include "./hydro/Minimal/hydro_parameters.h"
+#elif defined(GADGET2_SPH)
+#include "./hydro/Gadget2/hydro_parameters.h"
+#elif defined(HOPKINS_PE_SPH)
+#include "./hydro/PressureEntropy/hydro_parameters.h"
+#elif defined(HOPKINS_PU_SPH)
+#include "./hydro/PressureEnergy/hydro_parameters.h"
+#elif defined(HOPKINS_PU_SPH_MONAGHAN)
+#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h"
+#elif defined(DEFAULT_SPH)
+#include "./hydro/Default/hydro_parameters.h"
+#elif defined(GIZMO_MFV_SPH)
+#include "./hydro/GizmoMFV/hydro_parameters.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro_parameters.h"
+#elif defined(SHADOWFAX_SPH)
+#include "./hydro/Shadowswift/hydro_parameters.h"
+#elif defined(PLANETARY_SPH)
+#include "./hydro/Planetary/hydro_parameters.h"
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro_parameters.h"
+#elif defined(ANARCHY_PU_SPH)
+#include "./hydro/AnarchyPU/hydro_parameters.h"
+#else
+#error "Invalid choice of SPH variant"
+#endif
+
+#endif /* SWIFT_HYDRO_PARAMETERS_H */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index e5126fe67584a4401061cc795d105d9a0d1d7b07..604a61f1a79d2cd83755b9ba13d73efd4dba8318 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -44,34 +44,6 @@
#define hydro_props_default_init_temp 0.f
#define hydro_props_default_min_temp 0.f
#define hydro_props_default_H_ionization_temperature 1e4
-#define hydro_props_default_viscosity_alpha 0.8f
-
-#ifdef ANARCHY_PU_SPH
-/* This nasty #ifdef is only temporary until we separate the viscosity
- * and hydro components. If it is not removed by July 2019, shout at JB. */
-#undef hydro_props_default_viscosity_alpha
-#define hydro_props_default_viscosity_alpha \
- 0.1f /* Use a very low initial AV paramater for hydrodynamics tests */
-#define hydro_props_default_viscosity_alpha_min \
- 0.0f /* values NOT the same as Schaller+ 2015 */
-#define hydro_props_default_viscosity_alpha_max \
- 2.0f /* values taken from Schaller+ 2015 */
-#define hydro_props_default_viscosity_length \
- 0.25f /* values taken from Schaller+ 2015 */
-#else
-#define hydro_props_default_viscosity_alpha_min \
- 0.1f /* values taken from (price,2004), not used in legacy gadget mode */
-#define hydro_props_default_viscosity_alpha_max \
- 2.0f /* values taken from (price,2004), not used in legacy gadget mode */
-#define hydro_props_default_viscosity_length \
- 0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
-#endif /* ANARCHY_PU_SPH */
-
-/* Following values taken directly from the ANARCHY paper (Schaller+ 2015) */
-#define hydro_props_default_diffusion_alpha 0.0f
-#define hydro_props_default_diffusion_beta 0.01f
-#define hydro_props_default_diffusion_alpha_max 1.0f
-#define hydro_props_default_diffusion_alpha_min 0.0f
/**
* @brief Initialize the global properties of the hydro scheme.
@@ -166,37 +138,15 @@ void hydro_props_init(struct hydro_props *p,
/* Mean molecular mass for fully ionised gas */
p->mu_ionised = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
- /* Read the artificial viscosity parameters from the file, if they exist */
- p->viscosity.alpha = parser_get_opt_param_float(
- params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha);
-
- p->viscosity.alpha_max =
- parser_get_opt_param_float(params, "SPH:viscosity_alpha_max",
- hydro_props_default_viscosity_alpha_max);
-
- p->viscosity.alpha_min =
- parser_get_opt_param_float(params, "SPH:viscosity_alpha_min",
- hydro_props_default_viscosity_alpha_min);
-
- p->viscosity.length = parser_get_opt_param_float(
- params, "SPH:viscosity_length", hydro_props_default_viscosity_length);
+ /* Initialise the implementation-dependent viscosity parameters
+ * (see hydro/SCHEME/hydro_parameters.h for this implementation) */
+ viscosity_init(params, us, phys_const, &(p->viscosity));
/* Same for the thermal diffusion parameters */
- p->diffusion.alpha = parser_get_opt_param_float(
- params, "SPH:diffusion_alpha", hydro_props_default_diffusion_alpha);
-
- p->diffusion.beta = parser_get_opt_param_float(
- params, "SPH:diffusion_beta", hydro_props_default_diffusion_beta);
-
- p->diffusion.alpha_max =
- parser_get_opt_param_float(params, "SPH:diffusion_alpha_max",
- hydro_props_default_diffusion_alpha_max);
+ diffusion_init(params, us, phys_const, &(p->diffusion));
- p->diffusion.alpha_min =
- parser_get_opt_param_float(params, "SPH:diffusion_alpha_min",
- hydro_props_default_diffusion_alpha_min);
-
- /* Compute the initial energy (Note the temp. read is in internal units) */
+ /* Compute the initial energy (Note the temp. read is in internal units)
+ */
/* u_init = k_B T_init / (mu m_p (gamma - 1)) */
double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
u_init *= p->initial_temperature;
@@ -248,11 +198,12 @@ void hydro_props_print(const struct hydro_props *p) {
message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
- message(
- "Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
- "min: %.3f, length: %.3f.",
- p->viscosity.alpha, p->viscosity.alpha_max, p->viscosity.alpha_min,
- p->viscosity.length);
+ /* Print out the implementation-dependent viscosity parameters
+ * (see hydro/SCHEME/hydro_parameters.h for this implementation) */
+ viscosity_print(&(p->viscosity));
+
+ /* Same for the diffusion */
+ diffusion_print(&(p->diffusion));
message(
"Hydrodynamic integration: Max change of volume: %.2f "
@@ -316,22 +267,15 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
p->hydrogen_mass_fraction);
io_write_attribute_f(h_grpsph, "Hydrogen ionization transition temperature",
p->hydrogen_ionization_temperature);
- io_write_attribute_f(h_grpsph, "Alpha viscosity", p->viscosity.alpha);
- io_write_attribute_f(h_grpsph, "Alpha viscosity (max)",
- p->viscosity.alpha_max);
- io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
- p->viscosity.alpha_min);
- io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
- p->viscosity.length);
- io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
io_write_attribute_f(h_grpsph, "Max v_sig ratio (limiter)",
const_limiter_max_v_sig_ratio);
- io_write_attribute_f(h_grpsph, "Diffusion alpha", p->diffusion.alpha);
- io_write_attribute_f(h_grpsph, "Diffusion alpha (max)",
- p->diffusion.alpha_max);
- io_write_attribute_f(h_grpsph, "Diffusion alpha (min)",
- p->diffusion.alpha_min);
- io_write_attribute_f(h_grpsph, "Diffusion beta", p->diffusion.beta);
+
+ /* Write out the implementation-dependent viscosity parameters
+ * (see hydro/SCHEME/hydro_parameters.h for this implementation) */
+ viscosity_print_snapshot(h_grpsph, &(p->viscosity));
+
+ /* Same for the diffusion */
+ diffusion_print_snapshot(h_grpsph, &(p->diffusion));
}
#endif
@@ -370,15 +314,12 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
p->hydrogen_ionization_temperature =
hydro_props_default_H_ionization_temperature;
- p->viscosity.alpha = hydro_props_default_viscosity_alpha;
- p->viscosity.alpha_max = hydro_props_default_viscosity_alpha_max;
- p->viscosity.alpha_min = hydro_props_default_viscosity_alpha_min;
- p->viscosity.length = hydro_props_default_viscosity_length;
+ /* Initialise the implementation-dependent viscosity parameters
+ * (see hydro/SCHEME/hydro_parameters.h for this implementation) */
+ viscosity_init_no_hydro(&(p->viscosity));
- p->diffusion.alpha = hydro_props_default_diffusion_alpha;
- p->diffusion.beta = hydro_props_default_diffusion_beta;
- p->diffusion.alpha_max = hydro_props_default_diffusion_alpha_max;
- p->diffusion.alpha_min = hydro_props_default_diffusion_alpha_min;
+ /* Same for the diffusion */
+ diffusion_init_no_hydro(&(p->diffusion));
}
/**
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index afc8a4b87aeb39f6bff1e61ae39edf391c856b1b..40bf9652b560ab977f4dfaa29745d3e139454e68 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -32,6 +32,7 @@
#endif
/* Local includes. */
+#include "hydro_parameters.h"
#include "restart.h"
/* Forward declarations */
@@ -101,38 +102,10 @@ struct hydro_props {
float mu_ionised;
/*! Artificial viscosity parameters */
- struct {
- /*! For the fixed, simple case. Also used to set the initial AV
- coefficient for variable schemes. */
- float alpha;
-
- /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
- float alpha_max;
-
- /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
- float alpha_min;
-
- /*! The decay length of the artificial viscosity (used in M&M, etc.) */
- float length;
- } viscosity;
+ struct viscosity_global_data viscosity;
/*! Thermal diffusion parameters */
- struct {
-
- /*! Initialisation value, or the case for constant thermal diffusion coeffs
- */
- float alpha;
-
- /*! Tuning parameter for speed of ramp up/down */
- float beta;
-
- /*! Maximal value for alpha_diff */
- float alpha_max;
-
- /*! Minimal value for alpha_diff */
- float alpha_min;
-
- } diffusion;
+ struct diffusion_global_data diffusion;
};
void hydro_props_print(const struct hydro_props *p);
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 8d10c12a7bb8f39e7bbbb797b64750050b2866e0..e11bc753949b5fda8e1fd74755c36981e734f0f8 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1051,6 +1051,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
io_write_attribute_s(h_grp, "Code", "SWIFT");
time_t tm = time(NULL);
io_write_attribute_s(h_grp, "Snapshot date", ctime(&tm));
+ io_write_attribute_s(h_grp, "RunName", e->run_name);
/* GADGET-2 legacy values */
/* Number of particles of each type */
diff --git a/src/part.h b/src/part.h
index ac86c2ab533f08e5ce2cb7e0756a50dd92aec68b..253cec7b4edad08234a751ee20e6c05fcff7e6ec 100644
--- a/src/part.h
+++ b/src/part.h
@@ -77,6 +77,10 @@
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_part.h"
#define hydro_need_extra_init_loop 0
+#elif defined(ANARCHY_DU_SPH)
+#include "./hydro/AnarchyDU/hydro_part.h"
+#define hydro_need_extra_init_loop 0
+#define EXTRA_HYDRO_LOOP
#elif defined(ANARCHY_PU_SPH)
#include "./hydro/AnarchyPU/hydro_part.h"
#define hydro_need_extra_init_loop 0
diff --git a/src/runner.c b/src/runner.c
index 4f79264f94bf9af24fad28a6e234048153cac841..2b509b3557ebc78c84c763e73dc8f85428a11690 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1926,6 +1926,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const struct cosmology *cosmo = e->cosmology;
const struct chemistry_global_data *chemistry = e->chemistry;
const struct star_formation *star_formation = e->star_formation;
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+
const float hydro_h_max = e->hydro_properties->h_max;
const float hydro_h_min = e->hydro_properties->h_min;
const float eps = e->hydro_properties->h_tolerance;
@@ -2071,7 +2074,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Calculate the time-step for passing to hydro_prepare_force, used
* for the evolution of alpha factors (i.e. those involved in the
* artificial viscosity and thermal conduction terms) */
- const int with_cosmology = (e->policy & engine_policy_cosmology);
const double time_base = e->time_base;
const integertime_t ti_current = e->ti_current;
double dt_alpha;
@@ -2168,6 +2170,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
hydro_init_part(p, hs);
chemistry_init_part(p, chemistry);
star_formation_init_part(p, star_formation);
+ tracers_after_init(p, xp, e->internal_units, e->physical_constants,
+ with_cosmology, e->cosmology,
+ e->hydro_properties, e->cooling_func, e->time);
/* Off we go ! */
continue;
@@ -2222,7 +2227,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Calculate the time-step for passing to hydro_prepare_force, used for
* the evolution of alpha factors (i.e. those involved in the artificial
* viscosity and thermal conduction terms) */
- const int with_cosmology = (e->policy & engine_policy_cosmology);
const double time_base = e->time_base;
const integertime_t ti_current = e->ti_current;
double dt_alpha;
@@ -3990,6 +3994,9 @@ void *runner_main(void *data) {
/* Wait at the barrier. */
engine_barrier(e);
+ /* Can we go home yet? */
+ if (e->step_props & engine_step_prop_done) break;
+
/* Re-set the pointer to the previous task, as there is none. */
struct task *t = NULL;
struct task *prev = NULL;
diff --git a/src/serial_io.c b/src/serial_io.c
index 6ffac174347b6bfddc3c85f96d6ce54fdb6c0f58..629b4955066a47484f50b544be6caba55132cf27 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -918,6 +918,7 @@ void write_output_serial(struct engine* e, const char* baseName,
io_write_attribute_s(h_grp, "Code", "SWIFT");
time_t tm = time(NULL);
io_write_attribute_s(h_grp, "Snapshot date", ctime(&tm));
+ io_write_attribute_s(h_grp, "RunName", e->run_name);
/* GADGET-2 legacy values */
/* Number of particles of each type */
diff --git a/src/single_io.c b/src/single_io.c
index 57a83b4196542965b32ab099075ab44e1082ee6e..4dda6ad724cbcd111faa61a1d9d35a0098d85d95 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -768,6 +768,7 @@ void write_output_single(struct engine* e, const char* baseName,
io_write_attribute_s(h_grp, "Code", "SWIFT");
time_t tm = time(NULL);
io_write_attribute_s(h_grp, "Snapshot date", ctime(&tm));
+ io_write_attribute_s(h_grp, "RunName", e->run_name);
/* GADGET-2 legacy values */
/* Number of particles of each type */
diff --git a/src/space.c b/src/space.c
index 1c5023c3b79aab472b40756072c01870eb50dcd4..3ff86f439d20e0e63d2ab61314e5599aedb75dfc 100644
--- a/src/space.c
+++ b/src/space.c
@@ -4051,6 +4051,13 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
#endif
}
+ /* Overwrite the internal energy? */
+ if (u_init > 0.f) {
+ for (int k = 0; k < count; k++) {
+ hydro_set_init_internal_energy(&p[k], u_init);
+ }
+ }
+
/* Initialise the rest */
for (int k = 0; k < count; k++) {
@@ -4059,9 +4066,6 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
logger_part_data_init(&xp[k].logger_data);
#endif
- /* Overwrite the internal energy? */
- if (u_init > 0.f) hydro_set_init_internal_energy(&p[k], u_init);
-
/* Also initialise the chemistry */
chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k], &xp[k]);
@@ -4333,8 +4337,21 @@ void space_init_parts_mapper(void *restrict map_data, int count,
void *restrict extra_data) {
struct part *restrict parts = (struct part *)map_data;
- const struct hydro_space *restrict hs = (struct hydro_space *)extra_data;
- for (int k = 0; k < count; k++) hydro_init_part(&parts[k], hs);
+ const struct engine *restrict e = (struct engine *)extra_data;
+ const struct hydro_space *restrict hs = &e->s->hs;
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+
+ size_t ind = parts - e->s->parts;
+ struct xpart *restrict xparts = e->s->xparts + ind;
+
+ for (int k = 0; k < count; k++) {
+ hydro_init_part(&parts[k], hs);
+ chemistry_init_part(&parts[k], e->chemistry);
+ star_formation_init_part(&parts[k], e->star_formation);
+ tracers_after_init(&parts[k], &xparts[k], e->internal_units,
+ e->physical_constants, with_cosmology, e->cosmology,
+ e->hydro_properties, e->cooling_func, e->time);
+ }
}
/**
@@ -4349,7 +4366,7 @@ void space_init_parts(struct space *s, int verbose) {
if (s->nr_parts > 0)
threadpool_map(&s->e->threadpool, space_init_parts_mapper, s->parts,
- s->nr_parts, sizeof(struct part), 0, &s->hs);
+ s->nr_parts, sizeof(struct part), 0, s->e);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 7fc743725a58ed82c815dc97b3a7b89351c3f2c0..f05d3f83c41b5f26d7828f09aecd76685a5c55b6 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -111,8 +111,9 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->entropy = pressure / pow_gamma(density);
#elif defined(DEFAULT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
- defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
+ defined(ANARCHY_DU_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(PLANETARY_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
@@ -404,9 +405,9 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \
defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
- defined(ANARCHY_PU_SPH)
+ defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
0.f,
-#elif defined(ANARCHY_PU_SPH)
+#elif defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
main_cell->hydro.parts[pid].viscosity.div_v,
#else
main_cell->hydro.parts[pid].density.div_v,
@@ -430,7 +431,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
defined(HOPKINS_PU_SPH_MONAGHAN)
main_cell->hydro.parts[pid].force.v_sig, 0.f,
main_cell->hydro.parts[pid].u_dt
-#elif defined(ANARCHY_PU_SPH)
+#elif defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
main_cell->hydro.parts[pid].viscosity.v_sig, 0.f,
main_cell->hydro.parts[pid].u_dt
#else
diff --git a/tests/test27cells.c b/tests/test27cells.c
index cc34f503304feb56799a2d31baa3416b940202d3..e7c87cfe76382ed33a36744087e40a9161b60bfa 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -289,7 +289,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->hydro.parts[pid].density.rot_v[0],
main_cell->hydro.parts[pid].density.rot_v[1],
main_cell->hydro.parts[pid].density.rot_v[2]
-#elif defined(ANARCHY_PU_SPH)
+#elif defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
/* this is required because of the variable AV scheme */
main_cell->hydro.parts[pid].viscosity.div_v,
main_cell->hydro.parts[pid].density.rot_v[0],
@@ -334,7 +334,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
cj->hydro.parts[pjd].density.rot_v[0],
cj->hydro.parts[pjd].density.rot_v[1],
cj->hydro.parts[pjd].density.rot_v[2]
-#elif defined(ANARCHY_PU_SPH)
+#elif defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
/* this is required because of the variable AV scheme */
cj->hydro.parts[pjd].viscosity.div_v,
cj->hydro.parts[pjd].density.rot_v[0],
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index cbafd325f22b4e47c9c3ffcfaa5086c18196f39a..36bc04b6ffc794bcac556f969c0f710e67c4b24b 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -111,8 +111,9 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
/* Set the thermodynamic variable */
#if defined(GADGET2_SPH)
part->entropy = 1.f;
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
- defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
+ defined(ANARCHY_DU_SPH)
part->u = 1.f;
#elif defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
@@ -196,7 +197,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
p->density.rot_v[2] = 0.f;
p->density.div_v = 0.f;
#endif /* GADGET-2 */
-#ifdef MINIMAL_SPH
+#if defined(MINIMAL_SPH) || defined(ANARCHY_DU_SPH)
p->rho = 1.f;
p->density.rho_dh = 0.f;
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
@@ -219,7 +220,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
p->density.wcount_dh = 0.f;
#endif /* PRESSURE-ENERGY */
-#if defined(ANARCHY_PU_SPH)
+#if defined(ANARCHY_PU_SPH) || defined(ANARCHY_DU_SPH)
/* Initialise viscosity variables */
p->viscosity.alpha = 0.8;
p->viscosity.div_v = 0.f;
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index e14fddd640764c7e22a217fb483791494ba4fae0..0b1fb081b89b82223c86836bef6b0f68714eae2a 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -114,7 +114,7 @@ void prepare_force(struct part *parts, size_t count) {
#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
!defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) && \
!defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
- !defined(ANARCHY_PU_SPH)
+ !defined(ANARCHY_PU_SPH) && !defined(ANARCHY_DU_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -154,8 +154,9 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
p->force.v_sig, p->entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
p->force.v_sig, 0.f, p->force.u_dt
-#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
- defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH)
+#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
+ defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
+ defined(ANARCHY_DU_SPH)
p->force.v_sig, 0.f, p->u_dt
#else
0.f, 0.f, 0.f
@@ -560,7 +561,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
rhoiq[i] = pi_vec.rho;
grad_hiq[i] = pi_vec.force.f;
#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
- !defined(ANARCHY_PU_SPH)
+ !defined(ANARCHY_PU_SPH) && !defined(ANARCHY_DU_SPH)
pOrhoi2q[i] = pi_vec.force.P_over_rho2;
#endif
balsaraiq[i] = pi_vec.force.balsara;
@@ -574,7 +575,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
rhojq[i] = pj_vec[i].rho;
grad_hjq[i] = pj_vec[i].force.f;
#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
- !defined(ANARCHY_PU_SPH)
+ !defined(ANARCHY_PU_SPH) && !defined(ANARCHY_DU_SPH)
pOrhoj2q[i] = pj_vec[i].force.P_over_rho2;
#endif
balsarajq[i] = pj_vec[i].force.balsara;
@@ -657,7 +658,7 @@ void test_force_interactions(struct part test_part, struct part *parts,
VEC_HADD(h_dtSum, piq[0]->force.h_dt);
VEC_HMAX(v_sigSum, piq[0]->force.v_sig);
#if !defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
- !defined(ANARCHY_PU_SPH)
+ !defined(ANARCHY_PU_SPH) && !defined(ANARCHY_DU_SPH)
VEC_HADD(entropy_dtSum, piq[0]->entropy_dt);
#endif