diff --git a/.gitignore b/.gitignore index 438ea6c0b497f4939c67e35f4d4046d7346e33e1..344bad90125b6700c395e43f70a6922b4e1248ee 100644 --- a/.gitignore +++ b/.gitignore @@ -26,6 +26,8 @@ doc/html/ doc/latex/ doc/man/ doc/Doxyfile +doc/RTD/source/SubgridModels/*/*.png +doc/RTD/source/RadiativeTransfer/full_dependency_graph_RT.png examples/*/*/*.xmf examples/*/*/*.dat @@ -73,6 +75,14 @@ examples/SmallCosmoVolume/SmallCosmoVolume_cooling/snapshots/ examples/SmallCosmoVolume/SmallCosmoVolume_hydro/snapshots/ examples/SmallCosmoVolume/SmallCosmoVolume_cooling/CloudyData_UVB=HM2012.h5 examples/SmallCosmoVolume/SmallCosmoVolume_cooling/CloudyData_UVB=HM2012_shielded.h5 +examples/GEAR/AgoraDisk/CloudyData_UVB=HM2012.h5 +examples/GEAR/AgoraDisk/CloudyData_UVB=HM2012_shielded.h5 +examples/GEAR/AgoraDisk/chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 +examples/GEAR/AgoraCosmo/CloudyData_UVB=HM2012_shielded.h5 +examples/GEAR/AgoraCosmo/POPIIsw.h5 +examples/GEAR/ZoomIn/CloudyData_UVB=HM2012.h5 +examples/GEAR/ZoomIn/POPIIsw.h5 +examples/GEAR/ZoomIn/snap/ tests/testActivePair tests/testActivePair.sh @@ -198,6 +208,10 @@ theory/SPH/Flavours/sph_flavours.pdf theory/SPH/EoS/eos.pdf theory/SPH/*.pdf theory/paper_pasc/pasc_paper.pdf +theory/Multipoles/alpha_derivatives.pdf +theory/Multipoles/alpha_powers.pdf +theory/Multipoles/chi_derivatives.pdf +theory/Multipoles/sigma_derivatives.pdf theory/Multipoles/fmm.pdf theory/Multipoles/fmm_standalone.pdf theory/Multipoles/potential.pdf diff --git a/configure.ac b/configure.ac index 40bf209aad79e6b2352589c5f74f99635022a660..47aa6459bcb9af418d14a40881d6abf559ac9524 100644 --- a/configure.ac +++ b/configure.ac @@ -2240,7 +2240,7 @@ fi # chemistry function AC_ARG_WITH([chemistry], [AS_HELP_STRING([--with-chemistry=<function>], - [chemistry function @<:@none, GEAR_*, GEAR_DIFFUSION_*, QLA, EAGLE default: none@:>@ + [chemistry function @<:@none, GEAR_*, QLA, EAGLE default: none@:>@ For GEAR, you need to provide the number of elements (e.g. GEAR_10)] )], [with_chemistry="$withval"], @@ -2260,13 +2260,6 @@ case "$with_chemistry" in none) AC_DEFINE([CHEMISTRY_NONE], [1], [No chemistry function]) ;; - GEAR_DIFFUSION_*) - AC_DEFINE([CHEMISTRY_GEAR_DIFFUSION], [1], [Chemistry taken from the GEAR model including the metal diffusion]) - number_element=${with_chemistry##*_} - AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow]) - with_chemistry_name="GEAR (with $number_element elements and diffusion)" - with_chemistry="GEAR_DIFFUSION" - ;; GEAR_*) AC_DEFINE([CHEMISTRY_GEAR], [1], [Chemistry taken from the GEAR model]) number_element=${with_chemistry#*_} diff --git a/examples/Cooling/CoolingRates/cooling_rates.c b/examples/Cooling/CoolingRates/cooling_rates.c index 55528b5abfa9fed828ee5f3c9edddacd0ea3579a..4426b6404c98a7e7df0833eac37e91e01948c206 100644 --- a/examples/Cooling/CoolingRates/cooling_rates.c +++ b/examples/Cooling/CoolingRates/cooling_rates.c @@ -226,6 +226,9 @@ int main(int argc, char **argv) { // Init cosmology cosmology_init(params, &us, &internal_const, &cosmo); + // Init pressure floor + struct pressure_floor_props pressure_floor; + // Set redshift and associated quantities const float scale_factor = 1.0 / (1.0 + redshift); integertime_t ti_current = @@ -237,7 +240,7 @@ int main(int argc, char **argv) { cooling_init(params, &us, &internal_const, &hydro_properties, &cooling); cooling.H_reion_done = 1; cooling_print(&cooling); - cooling_update(&cosmo, &cooling, &s); + cooling_update(&cosmo, &pressure_floor, &cooling, &s); // Copy over the raw metals into the smoothed metals memcpy(&p.chemistry_data.smoothed_metal_mass_fraction, @@ -323,7 +326,9 @@ int main(int argc, char **argv) { unsigned long long cpufreq = 0; clocks_set_cpufreq(cpufreq); - message("This test is only defined for the EAGLE cooling model."); + message( + "This test is only defined for the EAGLE cooling model and with Gadget-2 " + "SPH."); return 0; } #endif diff --git a/examples/GEAR/AgoraCosmo/README b/examples/GEAR/AgoraCosmo/README index d30e9bf8a23204b1bb39a042a1db7ec84ed3675b..2c3415770d3fdd134f0cfcc17d5429933b5b9e01 100644 --- a/examples/GEAR/AgoraCosmo/README +++ b/examples/GEAR/AgoraCosmo/README @@ -7,3 +7,7 @@ We will use SWIFT to cancel the h- and a-factors from the ICs. The IC consists in a halo evolving to a virial mass of ∼ 1e12 Msun at z = 0 with a relatively quiescent merger history between z = 2 and 0. The cosmological parameters are taken from WMAP7/9+SNe+BAO. The target halo contains the highest-resolution dark matter particles particles of masses = 2.8e5 Msun. + +To run this example, SWIFT must be configured with the following options: + +./configure --with-chemistry=GEAR_10 --with-cooling=grackle_0 --with-pressure-floor=GEAR --with-stars=GEAR --with-star-formation=GEAR --with-feedback=GEAR \ No newline at end of file diff --git a/examples/GEAR/AgoraCosmo/agora_cosmo.yml b/examples/GEAR/AgoraCosmo/agora_cosmo.yml index b4c4b6b9dac625f835e2baca017d1e3f07bec5f5..63ec11ab3387915286f77c62a902ba3d36678452 100644 --- a/examples/GEAR/AgoraCosmo/agora_cosmo.yml +++ b/examples/GEAR/AgoraCosmo/agora_cosmo.yml @@ -93,17 +93,19 @@ GEARStarFormation: maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle n_stars_per_particle: 1 min_mass_frac: 0.5 + density_threashold: 1e-30 # Density threashold (in addition to the pressure floor) in g/cm3 GEARPressureFloor: jeans_factor: 8.75 GEARFeedback: - supernovae_energy_erg: 2.3e51 - yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 + supernovae_energy_erg: 1e51 + supernovae_efficiency: 2.3 + yields_table: POPIIsw.h5 discrete_yields: 0 - yields_table_first_stars: chemistry-PopIII.hdf5 # Table containing the yields of the first stars. + yields_table_first_stars: POPIIsw.h5 # Table containing the yields of the first stars. metallicity_max_first_stars: -1 # Maximal metallicity (in mass fraction) for a first star (-1 to deactivate). - elements: [Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). + elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). GEARChemistry: initial_metallicity: 1e-4 diff --git a/examples/GEAR/AgoraCosmo/getChemistryTable.sh b/examples/GEAR/AgoraCosmo/getChemistryTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..b10fd23964158ee7a38d352dbd0ddd9beb7bdd77 --- /dev/null +++ b/examples/GEAR/AgoraCosmo/getChemistryTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5 diff --git a/examples/GEAR/AgoraCosmo/getGrackleCoolingTable.sh b/examples/GEAR/AgoraCosmo/getGrackleCoolingTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..bf2e19000e1e0559b2b8fab6b82732337e467471 --- /dev/null +++ b/examples/GEAR/AgoraCosmo/getGrackleCoolingTable.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/GEAR/CloudyData_UVB=HM2012_shielded.h5 + + diff --git a/examples/GEAR/AgoraCosmo/getIC.sh b/examples/GEAR/AgoraCosmo/getIC.sh new file mode 100755 index 0000000000000000000000000000000000000000..ed1d1269c2b81d46d19ec11652af41c1ec22e916 --- /dev/null +++ b/examples/GEAR/AgoraCosmo/getIC.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraCosmo/agora_swift.hdf5 + + diff --git a/examples/GEAR/AgoraCosmo/run.sh b/examples/GEAR/AgoraCosmo/run.sh old mode 100644 new mode 100755 index 0ff0482c5c9285074e35e4aff93ff85d325506e4..b2d84a1396a3ed239847a0c01eefa42d15d1d3fb --- a/examples/GEAR/AgoraCosmo/run.sh +++ b/examples/GEAR/AgoraCosmo/run.sh @@ -1,29 +1,31 @@ #!/bin/bash -set -e - -# MUSIC binary -music=~/music/MUSIC -if test -f $music; then - echo "Using the following version of MUSIC $music." -else - echo "MUSIC is not found." - exit +# Get the initial conditions +if [ ! -e agora_swift.hdf5 ] +then + echo "Fetching the initial conditions" + ./getIC.sh fi -# Grab the cooling and yield tables if they are not present. + +# Get the Grackle cooling table if [ ! -e CloudyData_UVB=HM2012_shielded.h5 ] then - echo "Fetching tables..." - ../getChemistryTable.sh - ../../Cooling/getGrackleCoolingTable.sh + echo "Fetching the Cloudy tables required by Grackle..." + ./getGrackleCoolingTable.sh +fi + + +if [ ! -e POPIIsw.h5 ] +then + echo "Fetching the chemistry tables..." + ./getChemistryTable.sh fi -echo "Generating the initial conditions" -$music music.conf -echo "Converting the initial conditions into a SWIFT compatible format" -python3 convert_ic.py -echo "Running SWIFT" + ../../../swift --cooling --feedback --cosmology --limiter --sync --self-gravity --hydro --stars --star-formation --threads=24 agora_cosmo.yml 2>&1 | tee output.log + + + diff --git a/examples/GEAR/AgoraDisk/README b/examples/GEAR/AgoraDisk/README index c51d39ac476250eed6d2032053cc037e73eb188e..d230c751e40f23acb1c491fd75e3e6bbe113fa6b 100644 --- a/examples/GEAR/AgoraDisk/README +++ b/examples/GEAR/AgoraDisk/README @@ -1,3 +1,7 @@ +To run this example, SWIFT must be configured with the following options: + +./configure --with-chemistry=GEAR_10 --with-cooling=grackle_0 --with-pressure-floor=GEAR --with-stars=GEAR --with-star-formation=GEAR --with-feedback=GEAR + An analysis script for this test case can be found in the AGORA project repository: https://bitbucket.org/mornkr/agora-analysis-script/ diff --git a/examples/GEAR/AgoraDisk/agora_disk.yml b/examples/GEAR/AgoraDisk/agora_disk.yml index 6c17b9889dea37e3584e7e555f113722664a80ef..8846d17f8ed9e1da46c076d544c61889b40810d6 100644 --- a/examples/GEAR/AgoraDisk/agora_disk.yml +++ b/examples/GEAR/AgoraDisk/agora_disk.yml @@ -22,14 +22,14 @@ Scheduler: # Parameters governing the time integration TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). - time_end: 0.5 # The end time of the simulation (in internal units). + time_end: 0.25 # The end time of the simulation (in internal units). dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). dt_max: 0.1 # The maximal time-step size of the simulation (in internal units). max_dt_RMS_factor: 0.25 # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities. # Parameters governing the snapshots Snapshots: - basename: agora_disk # Common part of the name of output files + basename: snapshot # Common part of the name of output files time_first: 0. # Time of the first output (in internal units) delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) compression: 4 @@ -82,7 +82,7 @@ GrackleCooling: self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding max_steps: 1000 convergence_limit: 1e-2 - thermal_time_myr: 5 + thermal_time_myr: 0 GEARStarFormation: @@ -90,18 +90,20 @@ GEARStarFormation: maximal_temperature: 1e10 # Upper limit to the temperature of a star forming particle n_stars_per_particle: 1 min_mass_frac: 0.5 + density_threashold: 1.67e-23 # Density threashold in g/cm3 + GEARPressureFloor: jeans_factor: 10 GEARFeedback: - supernovae_efficiency: 1.0 - supernovae_energy_erg: 0.1e51 - yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 + supernovae_energy_erg: 1e51 + supernovae_efficiency: 0.1 + yields_table: POPIIsw.h5 discrete_yields: 0 - yields_table_first_stars: chemistry-PopIII.hdf5 # Table containing the yields of the first stars. - metallicity_max_first_stars: -1 # Maximal metallicity (in mass fraction) for a first star (-1 to deactivate). - elements: [Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). + yields_table_first_stars: POPIIsw.h5 # Table containing the yields of the first stars. + metallicity_max_first_stars: -5 # Maximal metallicity (in mass fraction) for a first star (-1 to deactivate). + elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). GEARChemistry: initial_metallicity: 1 diff --git a/examples/GEAR/AgoraDisk/changeType.py b/examples/GEAR/AgoraDisk/changeType.py deleted file mode 100644 index 4d0fbaa1a2c7736f0c9cbd399ce72970d0db5dc5..0000000000000000000000000000000000000000 --- a/examples/GEAR/AgoraDisk/changeType.py +++ /dev/null @@ -1,91 +0,0 @@ -#!/usr/bin/env python3 - -from h5py import File -from sys import argv -import numpy as np - -""" -Change particle types in order to match the implemented types -""" - -# number of particle type -N_type = 6 - -debug = 0 - - -def getOption(): - if len(argv) != 2: - raise IOError("You need to provide a filename") - - # get filename and read it - filename = argv[-1] - - return filename - - -def groupName(part_type): - return "PartType%i" % part_type - - -def changeType(f, old, new): - # check if directory exists - old_group = groupName(old) - if old_group not in f: - raise IOError("Cannot find group '%s'" % old) - old = f[old_group] - - new_group = groupName(new) - if new_group not in f: - f.create_group(new_group) - new = f[new_group] - - for name in old: - if debug: - print("Moving '%s' from '%s' to '%s'" % (name, old_group, new_group)) - - tmp = old[name][:] - del old[name] - if name in new: - new_tmp = new[name][:] - if debug: - print("Found previous data:", tmp.shape, new_tmp.shape) - tmp = np.append(tmp, new_tmp, axis=0) - del new[name] - - if debug: - print("With new shape:", tmp.shape) - new.create_dataset(name, tmp.shape) - new[name][:] = tmp - - del f[old_group] - - -def countPart(f): - npart = [] - for i in range(N_type): - name = groupName(i) - if name in f: - grp = f[groupName(i)] - N = grp["Masses"].shape[0] - else: - N = 0 - npart.append(N) - - f["Header"].attrs["NumPart_ThisFile"] = npart - f["Header"].attrs["NumPart_Total"] = npart - f["Header"].attrs["NumPart_Total_HighWord"] = [0] * N_type - - -if __name__ == "__main__": - filename = getOption() - - f = File(filename, "r+") - - changeType(f, 2, 1) - changeType(f, 3, 1) - changeType(f, 4, 1) - - countPart(f) - - f.close() diff --git a/examples/GEAR/AgoraDisk/checkSolution.py b/examples/GEAR/AgoraDisk/checkSolution.py new file mode 100644 index 0000000000000000000000000000000000000000..97941175caa112c4823f935ad18e2e19cf42e829 --- /dev/null +++ b/examples/GEAR/AgoraDisk/checkSolution.py @@ -0,0 +1,28 @@ +import h5py +import numpy as np + +expected_total_mass = 227413100.0 + +f = h5py.File("snapshot_0025.hdf5", "r") + +# get the total stellar mass +mass = f["StarsParticles/Masses"] +total_mass = sum(mass) * 1e10 + +error = np.fabs((expected_total_mass - total_mass) / expected_total_mass) + +if error < 0.01: + print(48 * "!") + print(r"Well done !") + print(r"The stellar mass is %g Msol," % (total_mass)) + print(r"The expected stellar mass is %g Msol" % (expected_total_mass)) + print(r"This represents an error of %d %% !" % (100 * error)) + print(48 * "!") +else: + print(48 * "!") + print(r"Too bad !") + print(r"The stellar mass is %g Msol," % (total_mass)) + print(r"While the expected stellar mass is %g Msol" % (expected_total_mass)) + print(r"This represents an error of %d %% !" % (100 * error)) + print(r"This is beyond the requirements.") + print(48 * "!") diff --git a/examples/GEAR/AgoraDisk/cleanupSwift.py b/examples/GEAR/AgoraDisk/cleanupSwift.py deleted file mode 100644 index 90537e31147e333366a7879d5bdd76c6e8d91788..0000000000000000000000000000000000000000 --- a/examples/GEAR/AgoraDisk/cleanupSwift.py +++ /dev/null @@ -1,49 +0,0 @@ -#!/usr/bin/env python3 - -# ./translate_particles.py filename output_name -from h5py import File -import sys -from shutil import copyfile - -NPartType = 1 -filename = sys.argv[-2] -out = sys.argv[-1] - -copyfile(filename, out) - -f = File(out) - -for i in range(6): - name = "PartType{}/ElementAbundances".format(i) - if name in f: - del f[name] - -for i in range(NPartType): - name = "PartType%i" % i - if name not in f: - continue - - grp = f[name + "/SmoothingLengths"] - grp[:] *= 1.823 - - # fix issue due to the name of densities - fields = [ - ("Density", "Densities"), - ("Entropies", "Entropies"), - ("InternalEnergy", "InternalEnergies"), - ("SmoothingLength", "SmoothingLengths"), - ] - - for field in fields: - if field[1] in f[name] and field[0] not in f[name]: - f[name + "/" + field[0]] = f[name + "/" + field[1]] - - -cosmo = f["Cosmology"].attrs -head = f["Header"].attrs -head["OmegaLambda"] = cosmo["Omega_lambda"] -head["Omega0"] = cosmo["Omega_b"] -head["HubbleParam"] = cosmo["H0 [internal units]"] -head["Time"] = head["Time"][0] - -f.close() diff --git a/examples/GEAR/AgoraDisk/getChemistryTable.sh b/examples/GEAR/AgoraDisk/getChemistryTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..b10fd23964158ee7a38d352dbd0ddd9beb7bdd77 --- /dev/null +++ b/examples/GEAR/AgoraDisk/getChemistryTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5 diff --git a/examples/GEAR/AgoraDisk/getGrackleCoolingTable.sh b/examples/GEAR/AgoraDisk/getGrackleCoolingTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..e3eb106240709c80151a48625567d2cd78e5f568 --- /dev/null +++ b/examples/GEAR/AgoraDisk/getGrackleCoolingTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5 diff --git a/examples/GEAR/AgoraDisk/getIC.sh b/examples/GEAR/AgoraDisk/getIC.sh index ca7fe3b251d568c58976650d158967214aa3e387..d53df79ffad11c9ded7508eb102d7009802b8072 100755 --- a/examples/GEAR/AgoraDisk/getIC.sh +++ b/examples/GEAR/AgoraDisk/getIC.sh @@ -1,9 +1,11 @@ #!/bin/bash if [ "$#" -ne 1 ]; then - echo "You need to provide the resolution (e.g. ./getIC.sh low)." - echo "The possible options are low, med and high." + echo "You need to provide the resolution (e.g. ./getIC.sh LOW)." + echo "The possible options are LOW and MED." exit fi -wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/$1.hdf5 +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5 +mv snap.$1.hdf5 agora_disk.hdf5 + diff --git a/examples/GEAR/AgoraDisk/run.sh b/examples/GEAR/AgoraDisk/run.sh index 9a32b86a479f2c32108ec973b0f1ccc2dd308b73..67371f8fb837daa0e353b3a1d5e56171c0a23d5d 100755 --- a/examples/GEAR/AgoraDisk/run.sh +++ b/examples/GEAR/AgoraDisk/run.sh @@ -2,16 +2,15 @@ # This example is based on the AGORA disk article (DOI: 10.3847/1538-4357/833/2/202) -# currently only the low resolution is available -sim=low +# set the resolution (LOW or MED) +sim=LOW -rm agora_disk_0*.hdf5 # make run.sh fail if a subcommand fails set -e # Generate the initial conditions if they are not present. -if [ ! -e $sim.hdf5 ] +if [ ! -e agora_disk.hdf5 ] then echo "Fetching initial glass file for the Sedov blast example..." ./getIC.sh $sim @@ -21,31 +20,22 @@ fi if [ ! -e CloudyData_UVB=HM2012.h5 ] then echo "Fetching the Cloudy tables required by Grackle..." - ../../Cooling/getGrackleCoolingTable.sh + ./getGrackleCoolingTable.sh fi -if [ ! -e chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 ] + +if [ ! -e POPIIsw.h5 ] then echo "Fetching the chemistry tables..." - ../getChemistryTable.sh + ./getChemistryTable.sh fi -# copy the initial conditions -cp $sim.hdf5 agora_disk.hdf5 -# Update the particle types -python3 changeType.py agora_disk.hdf5 # Run SWIFT ../../../swift --sync --limiter --cooling --hydro --self-gravity --star-formation --feedback --stars --threads=8 agora_disk.yml 2>&1 | tee output.log -echo "Changing smoothing length to be Gadget compatible" -python3 cleanupSwift.py agora_disk_0000.hdf5 agora_disk_IC.hdf5 -python3 cleanupSwift.py agora_disk_0050.hdf5 agora_disk_500Myr.hdf5 - -echo "Fetching GEAR solution..." -./getSolution.sh -echo "Plotting..." -python3 plotSolution.py +echo "check solution..." +python3 checkSolution.py diff --git a/examples/GEAR/ZoomIn/getChemistryTable.sh b/examples/GEAR/ZoomIn/getChemistryTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..b10fd23964158ee7a38d352dbd0ddd9beb7bdd77 --- /dev/null +++ b/examples/GEAR/ZoomIn/getChemistryTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5 diff --git a/examples/GEAR/ZoomIn/getGrackleCoolingTable.sh b/examples/GEAR/ZoomIn/getGrackleCoolingTable.sh new file mode 100755 index 0000000000000000000000000000000000000000..e3eb106240709c80151a48625567d2cd78e5f568 --- /dev/null +++ b/examples/GEAR/ZoomIn/getGrackleCoolingTable.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5 diff --git a/examples/GEAR/ZoomIn/getIC.sh b/examples/GEAR/ZoomIn/getIC.sh index ed95d78b23d81649a354571f83b7341404f05945..184f3255cddf7da39725c905581c8f5405317611 100755 --- a/examples/GEAR/ZoomIn/getIC.sh +++ b/examples/GEAR/ZoomIn/getIC.sh @@ -1,3 +1,5 @@ #!/bin/bash -# wget https://obswww.unige.ch/~lhausamm/swift/ZoomIC/h177.hdf5 -wget https://obswww.unige.ch/~lhausamm/swift/ZoomIC/h050.hdf5 + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/ZoomIn/h050.hdf5 + + diff --git a/examples/GEAR/ZoomIn/run.sh b/examples/GEAR/ZoomIn/run.sh index 141e6af26f5a1052690c665dab3fdf676a31f34d..a64dd745b4123135825d0435207fc15b3027b923 100755 --- a/examples/GEAR/ZoomIn/run.sh +++ b/examples/GEAR/ZoomIn/run.sh @@ -1,6 +1,30 @@ #!/bin/bash -echo "Fetching initial conditions for the zoom in example..." -./getIC.sh +if [ ! -e h050.hdf5 ] +then + echo "Fetching initial conditions for the zoom in example..." + ./getIC.sh +fi + + +# Get the Grackle cooling table +if [ ! -e CloudyData_UVB=HM2012.h5 ] +then + echo "Fetching the Cloudy tables required by Grackle..." + ./getGrackleCoolingTable.sh +fi + + +if [ ! -e POPIIsw.h5 ] +then + echo "Fetching the chemistry tables..." + ./getChemistryTable.sh +fi + + + + + + ../../../swift --cooling --feedback --cosmology --limiter --sync --self-gravity --hydro --stars --star-formation --threads=8 zoom_in.yml 2>&1 | tee output.log diff --git a/examples/GEAR/ZoomIn/zoom_in.yml b/examples/GEAR/ZoomIn/zoom_in.yml index 9e9135d04deecd15632517b123e24915a5ccfad5..30f4d32aafecc9bfc6da031ece2db9b116e18d12 100644 --- a/examples/GEAR/ZoomIn/zoom_in.yml +++ b/examples/GEAR/ZoomIn/zoom_in.yml @@ -35,7 +35,8 @@ Cosmology: # Parameters governing the snapshots Snapshots: - basename: snap/h050 # Common part of the name of output files + basename: h050 # Common part of the name of output files + subdir: snap time_first: 0. # Time of the first output (in internal units) delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) compression: 4 @@ -95,6 +96,7 @@ GrackleCooling: GEARStarFormation: star_formation_efficiency: 0.01 # star formation efficiency (c_*) maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle + density_threashold: 1.67e-25 # Density threashold in g/cm3 n_stars_per_particle: 4 min_mass_frac: 0.5 @@ -102,12 +104,13 @@ GEARPressureFloor: jeans_factor: 10 GEARFeedback: - supernovae_energy_erg: 0.1e51 - yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 + supernovae_energy_erg: 1e51 + supernovae_efficiency: 0.1 + yields_table: POPIIsw.h5 discrete_yields: 1 - yields_table_first_stars: chemistry-PopIII.hdf5 # Table containing the yields of the first stars. + yields_table_first_stars: POPIIsw.h5 # Table containing the yields of the first stars. metallicity_max_first_stars: -1 # Maximal metallicity (in mass fraction) for a first star (-1 to deactivate). - elements: [Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). + elements: [Fe, Mg, O, C, Al, Ca, Ba, Zn, Eu] # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N). GEARChemistry: initial_metallicity: 0 diff --git a/examples/RadiativeTransferTests/Advection_1D/plotSolution.py b/examples/RadiativeTransferTests/Advection_1D/plotSolution.py index 4feadac2a7e7c5371d89738739315ad1bcbe43b0..47e94a669a6502a16daa0b0345e2944c1f1093b6 100755 --- a/examples/RadiativeTransferTests/Advection_1D/plotSolution.py +++ b/examples/RadiativeTransferTests/Advection_1D/plotSolution.py @@ -116,10 +116,10 @@ def plot_photons(filename, energy_boundaries=None, flux_boundaries=None): ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"]) # Currently, SPHM1RT only works for frequency group = 4 in the code - # However, we only plot 3 frequency groups here, so - # we set ngroups = 3: + # However, we only plot 3 frequency groups here, so + # we set ngroups = 3: if scheme.startswith("SPH M1closure"): - ngroups = 3 + ngroups = 3 for g in range(ngroups): # workaround to access named columns data with swiftsimio visualisaiton diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h index dfd336c88afb0382acb88db3d64002707bb789e1..13ff2f7203422ee79ea138f9eb89c1e377638102 100644 --- a/src/black_holes/EAGLE/black_holes_iact.h +++ b/src/black_holes/EAGLE/black_holes_iact.h @@ -807,7 +807,8 @@ runner_iact_nonsym_bh_gas_feedback( const double u_new = u_init + delta_u; hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new); - hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new); + hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, + u_new); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); diff --git a/src/black_holes/SPIN_JET/black_holes_iact.h b/src/black_holes/SPIN_JET/black_holes_iact.h index e2ae49843d3893886f7e185fa634b2004d946ca1..f252f9d0350e7ee282af4f171547febb295f6938 100644 --- a/src/black_holes/SPIN_JET/black_holes_iact.h +++ b/src/black_holes/SPIN_JET/black_holes_iact.h @@ -930,7 +930,8 @@ runner_iact_nonsym_bh_gas_feedback( const double u_new = u_init + delta_u; hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new); - hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new); + hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, + u_new); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); diff --git a/src/cell_drift.c b/src/cell_drift.c index dd8a0e47065315882b8d24c250aec0c3b99fc6cd..f4a6287581ed44252ae9d04c19e51e8fe170948a 100644 --- a/src/cell_drift.c +++ b/src/cell_drift.c @@ -34,7 +34,6 @@ #include "lightcone/lightcone_array.h" #include "multipole.h" #include "neutrino.h" -#include "pressure_floor.h" #include "rt.h" #include "sink.h" #include "star_formation.h" @@ -344,7 +343,6 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force, mhd_init_part(p); black_holes_init_potential(&p->black_holes_data); chemistry_init_part(p, e->chemistry); - pressure_floor_init_part(p, xp); 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, diff --git a/src/cooling/COLIBRE/cooling.c b/src/cooling/COLIBRE/cooling.c index 5adfe42ef668dc0090429431bbfcea040b119ba8..a1739db9088f3e57a5c75e8d21ca216f8e203736 100644 --- a/src/cooling/COLIBRE/cooling.c +++ b/src/cooling/COLIBRE/cooling.c @@ -69,10 +69,12 @@ static const double bracket_factor = 1.5; * Also calls the additional H reionisation energy injection if need be. * * @param cosmo The current cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The space data, including a pointer to array of particles */ void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s) { /* Extra energy for reionization? */ @@ -85,7 +87,7 @@ void cooling_update(const struct cosmology *cosmo, if (s == NULL) error("Trying to do H reionization on an empty space!"); /* Inject energy to all particles */ - cooling_Hydrogen_reionization(cooling, cosmo, s); + cooling_Hydrogen_reionization(cooling, cosmo, pressure_floor, s); /* Flag that reionization happened */ cooling->H_reion_done = 1; @@ -719,6 +721,7 @@ static INLINE double bisection_iter( * @param cosmo The current cosmological model. * @param hydro_properties the hydro_props struct * @param floor_props Properties of the entropy floor. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -731,6 +734,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *p, struct xpart *xp, const float dt, const float dt_therm, const double time) { @@ -878,7 +882,8 @@ void cooling_cool_part(const struct phys_const *phys_const, /* Update the particle's u and du/dt */ hydro_set_physical_internal_energy(p, xp, cosmo, u_final); - hydro_set_drifted_physical_internal_energy(p, cosmo, u_final); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, + u_final); hydro_set_physical_internal_energy_dt(p, cosmo, 0.); } else { @@ -1382,11 +1387,12 @@ void cooling_split_part(struct part *p, struct xpart *xp, double n) { * * @param cooling The properties of the cooling model. * @param cosmo The cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param s The #space containing the particles. */ -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s) { +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s) { struct part *parts = s->parts; struct xpart *xparts = s->xparts; @@ -1410,7 +1416,8 @@ void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, const float new_u = old_u + extra_heat; hydro_set_physical_internal_energy(p, xp, cosmo, new_u); - hydro_set_drifted_physical_internal_energy(p, cosmo, new_u); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, + new_u); } } } @@ -1576,7 +1583,7 @@ void cooling_restore_tables(struct cooling_function_data *cooling, read_cooling_header(cooling); read_cooling_tables(cooling); - cooling_update(cosmo, cooling, /*space=*/NULL); + cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL); } /** diff --git a/src/cooling/COLIBRE/cooling.h b/src/cooling/COLIBRE/cooling.h index 1cc12cab7204d0059e662ab149612e933077b3e9..e7dba14eff40cb314d15a3d573ba8ac39b61bba9 100644 --- a/src/cooling/COLIBRE/cooling.h +++ b/src/cooling/COLIBRE/cooling.h @@ -33,10 +33,12 @@ struct xpart; struct cosmology; struct hydro_props; struct entropy_floor_properties; +struct pressure_floor_props; struct feedback_props; struct space; void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s); void cooling_cool_part(const struct phys_const *phys_const, @@ -44,6 +46,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *p, struct xpart *xp, const float dt, const float dt_therm, const double time); @@ -164,9 +167,9 @@ float cooling_get_radiated_energy(const struct xpart *xp); void cooling_split_part(struct part *p, struct xpart *xp, double n); -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s); +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s); void cooling_init_backend(struct swift_params *parameter_file, const struct unit_system *us, diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index aaa88c629871acf49da9799d5a4ab0e81d549cf7..4799488ff3318c7cdb3e592c4e27440a0815dad6 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -123,10 +123,12 @@ __attribute__((always_inline)) INLINE void get_redshift_index( * Also calls the additional H reionisation energy injection if need be. * * @param cosmo The current cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The space data, including a pointer to array of particles */ void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s) { /* Current redshift */ @@ -148,7 +150,7 @@ void cooling_update(const struct cosmology *cosmo, if (s == NULL) error("Trying to do H reionization on an empty space!"); /* Inject energy to all particles */ - cooling_Hydrogen_reionization(cooling, cosmo, s); + cooling_Hydrogen_reionization(cooling, cosmo, pressure_floor, s); /* Flag that reionization happened */ cooling->H_reion_done = 1; @@ -370,6 +372,7 @@ INLINE static double bisection_iter( * @param cosmo The current cosmological model. * @param hydro_properties the hydro_props struct * @param floor_props Properties of the entropy floor. + * @param pressure_floor Preopertes of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -383,6 +386,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *restrict p, struct xpart *restrict xp, const float dt, const float dt_therm, @@ -882,11 +886,12 @@ void cooling_split_part(struct part *p, struct xpart *xp, double n) { * * @param cooling The properties of the cooling model. * @param cosmo The cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param s The #space containing the particles. */ -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s) { +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s) { struct part *parts = s->parts; struct xpart *xparts = s->xparts; @@ -909,7 +914,7 @@ void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, const float new_u = old_u + extra_heat; hydro_set_physical_internal_energy(p, xp, cosmo, new_u); - hydro_set_drifted_physical_internal_energy(p, cosmo, new_u); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, new_u); } } @@ -1054,7 +1059,7 @@ void cooling_restore_tables(struct cooling_function_data *cooling, /* Force a re-read of the cooling tables */ cooling->z_index = -10; cooling->previous_z_index = eagle_cooling_N_redshifts - 2; - cooling_update(cosmo, cooling, /*space=*/NULL); + cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL); } /** diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h index 40b1fe7a5163d01e3b237121eed3985f2c1e2f4a..ce590c2cf72399e86dd73dc19d145f2a38f9def4 100644 --- a/src/cooling/EAGLE/cooling.h +++ b/src/cooling/EAGLE/cooling.h @@ -33,10 +33,12 @@ struct xpart; struct cosmology; struct hydro_props; struct entropy_floor_properties; +struct pressure_floor_props; struct space; struct phys_const; void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s); void cooling_cool_part(const struct phys_const *phys_const, @@ -44,17 +46,17 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, - struct part *restrict p, struct xpart *restrict xp, - const float dt, const float dt_therm, const double time); + struct part *p, struct xpart *xp, const float dt, + const float dt_therm, const double time); -float cooling_timestep(const struct cooling_function_data *restrict cooling, - const struct phys_const *restrict phys_const, - const struct cosmology *restrict cosmo, - const struct unit_system *restrict us, +float cooling_timestep(const struct cooling_function_data *cooling, + const struct phys_const *phys_const, + const struct cosmology *cosmo, + const struct unit_system *us, const struct hydro_props *hydro_props, - const struct part *restrict p, - const struct xpart *restrict xp); + const struct part *p, const struct xpart *xp); double cooling_get_electron_pressure( const struct phys_const *phys_const, const struct hydro_props *hydro_props, @@ -131,9 +133,9 @@ float cooling_get_radiated_energy(const struct xpart *restrict xp); void cooling_split_part(struct part *p, struct xpart *xp, double n); -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s); +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s); void cooling_init_backend(struct swift_params *parameter_file, const struct unit_system *us, diff --git a/src/cooling/QLA/cooling.c b/src/cooling/QLA/cooling.c index 60e50ce2f26f3caf1b24d5cb2412d06aa58ab357..dfe6e04edf357f76e72a573ffcd89b73c983a289 100644 --- a/src/cooling/QLA/cooling.c +++ b/src/cooling/QLA/cooling.c @@ -47,6 +47,7 @@ #include "parser.h" #include "part.h" #include "physical_constants.h" +#include "pressure_floor.h" #include "space.h" #include "units.h" @@ -67,10 +68,12 @@ static const double bracket_factor = 1.5; * Also calls the additional H reionisation energy injection if need be. * * @param cosmo The current cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The space data, including a pointer to array of particles */ void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s) { /* Extra energy for reionization? */ @@ -83,7 +86,7 @@ void cooling_update(const struct cosmology *cosmo, if (s == NULL) error("Trying to do H reionization on an empty space!"); /* Inject energy to all particles */ - cooling_Hydrogen_reionization(cooling, cosmo, s); + cooling_Hydrogen_reionization(cooling, cosmo, pressure_floor, s); /* Flag that reionization happened */ cooling->H_reion_done = 1; @@ -363,7 +366,6 @@ static INLINE double bisection_iter( const struct cooling_function_data *cooling, const float abundance_ratio[qla_cooling_N_elementtypes], double dt_cgs, long long ID) { - /* Bracketing */ double u_lower_cgs = max(u_ini_cgs, cooling->umin_cgs); double u_upper_cgs = max(u_ini_cgs, cooling->umin_cgs); @@ -548,6 +550,7 @@ static INLINE double bisection_iter( * @param cosmo The current cosmological model. * @param hydro_properties the hydro_props struct * @param floor_props Properties of the entropy floor. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -560,6 +563,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *p, struct xpart *xp, const float dt, const float dt_therm, const double time) { @@ -702,7 +706,8 @@ void cooling_cool_part(const struct phys_const *phys_const, /* Update the particle's u and du/dt */ hydro_set_physical_internal_energy(p, xp, cosmo, u_final); - hydro_set_drifted_physical_internal_energy(p, cosmo, u_final); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, + u_final); hydro_set_physical_internal_energy_dt(p, cosmo, 0.); } else { @@ -955,11 +960,12 @@ void cooling_split_part(struct part *p, struct xpart *xp, double n) {} * * @param cooling The properties of the cooling model. * @param cosmo The cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param s The #space containing the particles. */ -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s) { +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s) { struct part *parts = s->parts; struct xpart *xparts = s->xparts; @@ -982,7 +988,7 @@ void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, const float new_u = old_u + extra_heat; hydro_set_physical_internal_energy(p, xp, cosmo, new_u); - hydro_set_drifted_physical_internal_energy(p, cosmo, new_u); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, new_u); } } @@ -1131,7 +1137,7 @@ void cooling_restore_tables(struct cooling_function_data *cooling, read_cooling_header(cooling); read_cooling_tables(cooling); - cooling_update(cosmo, cooling, /*space=*/NULL); + cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL); } /** diff --git a/src/cooling/QLA/cooling.h b/src/cooling/QLA/cooling.h index 6e735de2ac24959f73fb7c5770cabb64e00be369..1e3441b90bd9eb6acf7281c50b24cc0ca560b065 100644 --- a/src/cooling/QLA/cooling.h +++ b/src/cooling/QLA/cooling.h @@ -34,10 +34,12 @@ struct xpart; struct cosmology; struct hydro_props; struct entropy_floor_properties; +struct pressure_floor_props; struct phys_const; struct space; void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s); void cooling_cool_part(const struct phys_const *phys_const, @@ -45,6 +47,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *p, struct xpart *xp, const float dt, const float dt_therm, const double time); @@ -68,7 +71,6 @@ float cooling_get_temperature_from_gas( const struct cooling_function_data *cooling, const float rho_phys, const float XH, const float logZZsol, const float u_phys, const int HII_region); - float cooling_get_temperature(const struct phys_const *phys_const, const struct hydro_props *hydro_props, const struct unit_system *us, @@ -165,9 +167,9 @@ float cooling_get_radiated_energy(const struct xpart *xp); void cooling_split_part(struct part *p, struct xpart *xp, double n); -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s); +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s); void cooling_init_backend(struct swift_params *parameter_file, const struct unit_system *us, diff --git a/src/cooling/QLA_EAGLE/cooling.c b/src/cooling/QLA_EAGLE/cooling.c index eacc7d02c5af41de9d84ea45acaf0d52e60c2d75..8e7bb884c912ad24466f8241392d236d2c6740f5 100644 --- a/src/cooling/QLA_EAGLE/cooling.c +++ b/src/cooling/QLA_EAGLE/cooling.c @@ -123,10 +123,12 @@ __attribute__((always_inline)) INLINE void get_redshift_index( * Also calls the additional H reionisation energy injection if need be. * * @param cosmo The current cosmological model. + * @param pressure_floor The properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The space data, including a pointer to array of particles */ void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s) { /* Current redshift */ @@ -148,7 +150,7 @@ void cooling_update(const struct cosmology *cosmo, if (s == NULL) error("Trying to do H reionization on an empty space!"); /* Inject energy to all particles */ - cooling_Hydrogen_reionization(cooling, cosmo, s); + cooling_Hydrogen_reionization(cooling, cosmo, pressure_floor, s); /* Flag that reionization happened */ cooling->H_reion_done = 1; @@ -373,6 +375,7 @@ INLINE static double bisection_iter( * @param cosmo The current cosmological model. * @param hydro_properties the hydro_props struct * @param floor_props Properties of the entropy floor. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -386,6 +389,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *restrict p, struct xpart *restrict xp, const float dt, const float dt_therm, @@ -897,11 +901,12 @@ void cooling_split_part(struct part *p, struct xpart *xp, double n) {} * * @param cooling The properties of the cooling model. * @param cosmo The cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param s The #space containing the particles. */ -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s) { +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s) { struct part *parts = s->parts; struct xpart *xparts = s->xparts; @@ -924,7 +929,7 @@ void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, const float new_u = old_u + extra_heat; hydro_set_physical_internal_energy(p, xp, cosmo, new_u); - hydro_set_drifted_physical_internal_energy(p, cosmo, new_u); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, new_u); } } @@ -1063,7 +1068,7 @@ void cooling_restore_tables(struct cooling_function_data *cooling, /* Force a re-read of the cooling tables */ cooling->z_index = -10; cooling->previous_z_index = qla_eagle_cooling_N_redshifts - 2; - cooling_update(cosmo, cooling, /*space=*/NULL); + cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL); } /** diff --git a/src/cooling/QLA_EAGLE/cooling.h b/src/cooling/QLA_EAGLE/cooling.h index 8fb5b066b0911b4df6f14ee510cd61037109c9d3..0abafc6323be12715c1b5f2de453c99b6dccf58a 100644 --- a/src/cooling/QLA_EAGLE/cooling.h +++ b/src/cooling/QLA_EAGLE/cooling.h @@ -33,10 +33,12 @@ struct xpart; struct cosmology; struct hydro_props; struct entropy_floor_properties; +struct pressure_floor_props; struct space; struct phys_const; void cooling_update(const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct cooling_function_data *cooling, struct space *s); void cooling_cool_part(const struct phys_const *phys_const, @@ -44,6 +46,7 @@ void cooling_cool_part(const struct phys_const *phys_const, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor, const struct cooling_function_data *cooling, struct part *restrict p, struct xpart *restrict xp, const float dt, const float dt_therm, const double time); @@ -139,9 +142,9 @@ float cooling_get_radiated_energy(const struct xpart *restrict xp); void cooling_split_part(struct part *p, struct xpart *xp, double n); -void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, - const struct cosmology *cosmo, - struct space *s); +void cooling_Hydrogen_reionization( + const struct cooling_function_data *cooling, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, struct space *s); void cooling_init_backend(struct swift_params *parameter_file, const struct unit_system *us, diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h index 4076ac07ce6cbc434510da3b47469a05668f2c9e..22c90bc679437cee9e7057263b930bdc4b5a6286 100644 --- a/src/cooling/const_du/cooling.h +++ b/src/cooling/const_du/cooling.h @@ -55,12 +55,14 @@ * given time-step or redshift. * * @param cosmo The current cosmological model. + * @param pressure_floor The properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The #space containing all the particles. */ -INLINE static void cooling_update(const struct cosmology* cosmo, - struct cooling_function_data* cooling, - struct space* s) { +INLINE static void cooling_update( + const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, + struct cooling_function_data* cooling, struct space* s) { // Add content if required. } @@ -75,6 +77,7 @@ INLINE static void cooling_update(const struct cosmology* cosmo, * @param cosmo The current cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props Properties of the entropy floor. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -89,6 +92,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( const struct cosmology* restrict cosmo, const struct hydro_props* hydro_props, const struct entropy_floor_properties* floor_props, + const struct pressure_floor_props* pressure_floor, const struct cooling_function_data* restrict cooling, struct part* restrict p, struct xpart* restrict xp, const float dt, const float dt_therm, const double time) { diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 64a2d35040e25032961d906c38d02e0f9325ef2f..c204996e40d69f8837d7fbeddb4dbc125e98718f 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -51,12 +51,14 @@ * given time-step or redshift. * * @param cosmo The current cosmological model. + * @param pressure_floor The properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The #space containing all the particles. */ -INLINE static void cooling_update(const struct cosmology* cosmo, - struct cooling_function_data* cooling, - struct space* s) { +INLINE static void cooling_update( + const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, + struct cooling_function_data* cooling, struct space* s) { // Add content if required. } @@ -104,6 +106,7 @@ __attribute__((always_inline)) INLINE static double cooling_rate_cgs( * @param cosmo The current cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props Properties of the entropy floor. + * @param pressure_floor The properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the particle' extended data. @@ -113,14 +116,12 @@ __attribute__((always_inline)) INLINE static double cooling_rate_cgs( * units. */ __attribute__((always_inline)) INLINE static void cooling_cool_part( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct hydro_props* hydro_props, + const struct phys_const* phys_const, const struct unit_system* us, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, const struct entropy_floor_properties* floor_props, - const struct cooling_function_data* restrict cooling, - struct part* restrict p, struct xpart* restrict xp, const float dt, - const float dt_therm, const double time) { + const struct pressure_floor_props* pressure_floor, + const struct cooling_function_data* cooling, struct part* p, + struct xpart* xp, const float dt, const float dt_therm, const double time) { /* Nothing to do here? */ if (dt == 0.) return; @@ -182,7 +183,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( const float u_new_com = u_old_com + total_du_dt * dt_therm; const float u_new_phys = u_new_com * cosmo->a_factor_internal_energy; hydro_set_physical_internal_energy(p, xp, cosmo, u_new_phys); - hydro_set_drifted_physical_internal_energy(p, cosmo, u_new_phys); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, + u_new_phys); hydro_set_physical_internal_energy_dt(p, cosmo, 0.); } else { /* Update the internal energy time derivative */ diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c index 48e2072a8f23d7e95cddccd8fa9cfc455479eaea..2c4e093c9b8f986c1a7abec36aa68d17eb8cadc6 100644 --- a/src/cooling/grackle/cooling.c +++ b/src/cooling/grackle/cooling.c @@ -49,15 +49,32 @@ #define GRACKLE_NPART 1 #define GRACKLE_RANK 3 +gr_float cooling_new_energy(const struct phys_const* phys_const, + const struct unit_system* us, + const struct cosmology* cosmo, + const struct hydro_props* hydro_properties, + const struct cooling_function_data* cooling, + const struct part* p, struct xpart* xp, double dt, + double dt_therm); + +gr_float cooling_time(const struct phys_const* phys_const, + const struct unit_system* us, + const struct hydro_props* hydro_properties, + const struct cosmology* cosmo, + const struct cooling_function_data* cooling, + const struct part* p, struct xpart* xp); + /** * @brief Common operations performed on the cooling function at a * given time-step or redshift. * * @param cosmo The current cosmological model. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param s The #space containing all the particles. */ void cooling_update(const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, struct cooling_function_data* cooling, struct space* s) { /* set current time */ if (cooling->redshift == -1) @@ -151,13 +168,12 @@ int cooling_converged(const struct xpart* restrict xp, * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. */ -void cooling_compute_equilibrium( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct hydro_props* hydro_properties, - const struct cosmology* restrict cosmo, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, struct xpart* restrict xp) { +void cooling_compute_equilibrium(const struct phys_const* phys_const, + const struct unit_system* us, + const struct hydro_props* hydro_properties, + const struct cosmology* cosmo, + const struct cooling_function_data* cooling, + const struct part* p, struct xpart* xp) { /* TODO: this can fail spectacularly and needs to be replaced. */ @@ -212,13 +228,12 @@ void cooling_compute_equilibrium( * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. */ -void cooling_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, +void cooling_first_init_part(const struct phys_const* phys_const, + const struct unit_system* us, const struct hydro_props* hydro_props, - const struct cosmology* restrict cosmo, + const struct cosmology* cosmo, const struct cooling_function_data* cooling, - const struct part* restrict p, - struct xpart* restrict xp) { + const struct part* p, struct xpart* xp) { xp->cooling_data.radiated_energy = 0.f; xp->cooling_data.time_last_event = -cooling->thermal_time; @@ -264,7 +279,7 @@ void cooling_first_init_part(const struct phys_const* restrict phys_const, * * @param xp The extended particle data */ -float cooling_get_radiated_energy(const struct xpart* restrict xp) { +float cooling_get_radiated_energy(const struct xpart* xp) { return xp->cooling_data.radiated_energy; } @@ -597,14 +612,13 @@ void cooling_apply_self_shielding( * * @return du / dt */ -gr_float cooling_new_energy( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct hydro_props* hydro_props, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, struct xpart* restrict xp, double dt, - double dt_therm) { +gr_float cooling_new_energy(const struct phys_const* phys_const, + const struct unit_system* us, + const struct cosmology* cosmo, + const struct hydro_props* hydro_props, + const struct cooling_function_data* cooling, + const struct part* p, struct xpart* xp, double dt, + double dt_therm) { /* set current time */ code_units units = cooling->units; @@ -676,13 +690,12 @@ gr_float cooling_new_energy( * * @return cooling time */ -gr_float cooling_time(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, +gr_float cooling_time(const struct phys_const* phys_const, + const struct unit_system* us, const struct hydro_props* hydro_props, - const struct cosmology* restrict cosmo, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, - struct xpart* restrict xp) { + const struct cosmology* cosmo, + const struct cooling_function_data* cooling, + const struct part* p, struct xpart* xp) { /* set current time */ code_units units = cooling->units; @@ -747,6 +760,7 @@ gr_float cooling_time(const struct phys_const* restrict phys_const, * @param cosmo The current cosmological model. * @param hydro_props The #hydro_props. * @param floor_props Properties of the entropy floor. + * @param pressure_floor Properties of the pressure floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the particle' extended data. @@ -755,15 +769,15 @@ gr_float cooling_time(const struct phys_const* restrict phys_const, * @param time The current time (since the Big Bang or start of the run) in * internal units. */ -void cooling_cool_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, +void cooling_cool_part(const struct phys_const* phys_const, + const struct unit_system* us, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, const struct entropy_floor_properties* floor_props, - const struct cooling_function_data* restrict cooling, - struct part* restrict p, struct xpart* restrict xp, - const double dt, const double dt_therm, - const double time) { + const struct pressure_floor_props* pressure_floor, + const struct cooling_function_data* cooling, + struct part* p, struct xpart* xp, const double dt, + const double dt_therm, const double time) { /* Nothing to do here? */ if (dt == 0.) return; @@ -822,13 +836,12 @@ void cooling_cool_part(const struct phys_const* restrict phys_const, * @param p #part data. * @param xp Pointer to the #xpart data. */ -float cooling_get_temperature( - const struct phys_const* restrict phys_const, - const struct hydro_props* restrict hydro_props, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, const struct xpart* restrict xp) { +float cooling_get_temperature(const struct phys_const* phys_const, + const struct hydro_props* hydro_props, + const struct unit_system* us, + const struct cosmology* cosmo, + const struct cooling_function_data* cooling, + const struct part* p, const struct xpart* xp) { // TODO use the grackle library /* Physical constants */ @@ -892,13 +905,12 @@ double Cooling_get_ycompton(const struct phys_const* phys_const, * @param p Pointer to the particle data. * @param xp Pointer to the particle extra data */ -float cooling_timestep(const struct cooling_function_data* restrict cooling, - const struct phys_const* restrict phys_const, - const struct cosmology* restrict cosmo, - const struct unit_system* restrict us, +float cooling_timestep(const struct cooling_function_data* cooling, + const struct phys_const* phys_const, + const struct cosmology* cosmo, + const struct unit_system* us, const struct hydro_props* hydro_props, - const struct part* restrict p, - const struct xpart* restrict xp) { + const struct part* p, const struct xpart* xp) { return FLT_MAX; } diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index 5173e9adb4c96bac2c01b557e772f68b46873c04..3584bd7f188c9a6ea5bf0501635f771753318d10 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -24,49 +24,36 @@ * @brief Cooling using the GRACKLE 3.1.1 library. */ -/* Some standard headers. */ -#include <fenv.h> -#include <float.h> -#include <math.h> - -/* The grackle library itself */ -#include <grackle.h> - -/* Local includes. */ -#include "chemistry.h" -#include "cooling_io.h" +/* Local includes */ #include "cooling_properties.h" -#include "entropy_floor.h" #include "error.h" -#include "hydro.h" -#include "parser.h" -#include "part.h" -#include "physical_constants.h" -#include "units.h" +#include "inline.h" + +struct part; +struct xpart; +struct cosmology; +struct hydro_props; +struct entropy_floor_properties; +struct pressure_floor_props; +struct space; +struct phys_const; +struct unit_system; +struct swift_params; /* need to rework (and check) code if changed */ #define GRACKLE_NPART 1 #define GRACKLE_RANK 3 void cooling_update(const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, struct cooling_function_data* cooling, struct space* s); -void cooling_print_fractions(const struct xpart* restrict xp); -int cooling_converged(const struct xpart* restrict xp, - const struct xpart* restrict old, const float limit); -void cooling_compute_equilibrium( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct hydro_props* hydro_properties, - const struct cosmology* restrict cosmo, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, struct xpart* restrict xp); -void cooling_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, + +void cooling_first_init_part(const struct phys_const* phys_const, + const struct unit_system* us, const struct hydro_props* hydro_properties, - const struct cosmology* restrict cosmo, + const struct cosmology* cosmo, const struct cooling_function_data* cooling, - const struct part* restrict p, - struct xpart* restrict xp); + const struct part* p, struct xpart* xp); /** * @brief Returns the subgrid temperature of a particle. @@ -99,51 +86,12 @@ INLINE static float cooling_get_subgrid_density(const struct part* p, float cooling_get_radiated_energy(const struct xpart* restrict xp); void cooling_print_backend(const struct cooling_function_data* cooling); -void cooling_copy_to_grackle1(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho, - gr_float species_densities[12]); -void cooling_copy_to_grackle2(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho, - gr_float species_densities[12]); -void cooling_copy_to_grackle3(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho, - gr_float species_densities[12]); -void cooling_copy_from_grackle1(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho); -void cooling_copy_from_grackle2(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho); -void cooling_copy_from_grackle3(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho); -void cooling_copy_to_grackle(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho, - gr_float species_densities[12]); -void cooling_copy_from_grackle(grackle_field_data* data, const struct part* p, - struct xpart* xp, gr_float rho); -void cooling_apply_self_shielding( - const struct cooling_function_data* restrict cooling, - chemistry_data* restrict chemistry, const struct part* restrict p, - const struct cosmology* cosmo); -gr_float cooling_new_energy( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct hydro_props* hydro_properties, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, struct xpart* restrict xp, double dt, - double dt_therm); - -gr_float cooling_time(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct hydro_props* hydro_properties, - const struct cosmology* restrict cosmo, - const struct cooling_function_data* restrict cooling, - const struct part* restrict p, struct xpart* restrict xp); - void cooling_cool_part(const struct phys_const* restrict phys_const, const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct hydro_props* hydro_properties, const struct entropy_floor_properties* floor_props, + const struct pressure_floor_props* pressure_floor, const struct cooling_function_data* restrict cooling, struct part* restrict p, struct xpart* restrict xp, const double dt, const double dt_therm, diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h index 69128aa3577c072b1c2633cabef820acdb21d486..404f5ff48497f4461c1a1ab0cab8a5a7b163c9c8 100644 --- a/src/cooling/none/cooling.h +++ b/src/cooling/none/cooling.h @@ -43,11 +43,13 @@ * * @param cosmo The current cosmological model. * @param cooling The #cooling_function_data used in the run. + * @param pressure_floor Properties of the pressure floor. * @param s The #space containing all the particles. */ -INLINE static void cooling_update(const struct cosmology* cosmo, - struct cooling_function_data* cooling, - struct space* s) { +INLINE static void cooling_update( + const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, + struct cooling_function_data* cooling, struct space* s) { // Add content if required. } @@ -60,6 +62,8 @@ INLINE static void cooling_update(const struct cosmology* cosmo, * @param us The internal system of units. * @param cosmo The current cosmological model. * @param hydro_props The properties of the hydro scheme. + * @param hydro_properties the hydro_props struct + * @param floor_props Properties of the entropy floor. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. @@ -69,14 +73,13 @@ INLINE static void cooling_update(const struct cosmology* cosmo, * internal units. */ __attribute__((always_inline)) INLINE static void cooling_cool_part( - const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - const struct hydro_props* hydro_props, + const struct phys_const* phys_const, const struct unit_system* us, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, const struct entropy_floor_properties* floor_props, - const struct cooling_function_data* restrict cooling, - struct part* restrict p, struct xpart* restrict xp, const float dt, - const float dt_therm, const double time) {} + const struct pressure_floor_props* pressure_floor, + const struct cooling_function_data* cooling, struct part* p, + struct xpart* xp, const float dt, const float dt_therm, const double time) { +} /** * @brief Computes the cooling time-step. diff --git a/src/drift.h b/src/drift.h index 4f7d6768b78a7c2547df9598abec9338f1d7017e..1e5764d9004c886af045413e04f692765a80f4a9 100644 --- a/src/drift.h +++ b/src/drift.h @@ -148,7 +148,8 @@ __attribute__((always_inline)) INLINE static void drift_part( const struct cosmology *cosmo = e->cosmology; const struct hydro_props *hydro_props = e->hydro_properties; - const struct entropy_floor_properties *floor = e->entropy_floor; + const struct entropy_floor_properties *entropy_floor = e->entropy_floor; + const struct pressure_floor_props *pressure_floor = e->pressure_floor_props; #ifdef SWIFT_DEBUG_CHECKS if (p->ti_drift != ti_old) @@ -200,8 +201,9 @@ __attribute__((always_inline)) INLINE static void drift_part( /* Predict the values of the extra fields */ hydro_predict_extra(p, xp, dt_drift, dt_therm, dt_kick_grav, cosmo, - hydro_props, floor); - mhd_predict_extra(p, xp, dt_drift, dt_therm, cosmo, hydro_props, floor); + hydro_props, entropy_floor, pressure_floor); + mhd_predict_extra(p, xp, dt_drift, dt_therm, cosmo, hydro_props, + entropy_floor); rt_predict_extra(p, xp, dt_drift); /* Compute offsets since last cell construction */ diff --git a/src/engine.c b/src/engine.c index 9bae8120d17c28b5ee146426e5e18516cf25c597..c18b5eca56c27a69f389b44dcc0c0f86fe60b785 100644 --- a/src/engine.c +++ b/src/engine.c @@ -1950,7 +1950,8 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs, /* Update the cooling function */ if ((e->policy & engine_policy_cooling) || (e->policy & engine_policy_temperature)) - cooling_update(e->cosmology, e->cooling_func, e->s); + cooling_update(e->cosmology, e->pressure_floor_props, e->cooling_func, + e->s); #ifdef WITH_CSDS if (e->policy & engine_policy_csds) { @@ -2370,7 +2371,8 @@ int engine_step(struct engine *e) { /* Update the cooling function */ if ((e->policy & engine_policy_cooling) || (e->policy & engine_policy_temperature)) - cooling_update(e->cosmology, e->cooling_func, e->s); + cooling_update(e->cosmology, e->pressure_floor_props, e->cooling_func, + e->s); /* Update the softening lengths */ if (e->policy & engine_policy_self_gravity) diff --git a/src/feedback/EAGLE_kinetic/feedback_iact.h b/src/feedback/EAGLE_kinetic/feedback_iact.h index 9efe800fbf421637f1f8b302f4b15ad3352270ed..7fef5f8708a71a238fa1f8d0baf8856e398ca449 100644 --- a/src/feedback/EAGLE_kinetic/feedback_iact.h +++ b/src/feedback/EAGLE_kinetic/feedback_iact.h @@ -509,7 +509,8 @@ runner_iact_nonsym_feedback_apply( /* Do the energy injection. */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_enrich); - hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new_enrich); + hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, + u_new_enrich); /* Synchronize the particle on the timeline if we've got some extra * thermal energy from the SNII kicks that have not occured */ diff --git a/src/feedback/EAGLE_thermal/feedback_iact.h b/src/feedback/EAGLE_thermal/feedback_iact.h index 89c91a63984be432776494b35e0794522da1246c..341f8217535a0835b8572bc88e7b1393fdba051f 100644 --- a/src/feedback/EAGLE_thermal/feedback_iact.h +++ b/src/feedback/EAGLE_thermal/feedback_iact.h @@ -381,7 +381,8 @@ runner_iact_nonsym_feedback_apply( /* Do the energy injection. */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new_enrich); - hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new_enrich); + hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, + u_new_enrich); /* Finally, SNII stochastic feedback */ @@ -412,7 +413,8 @@ runner_iact_nonsym_feedback_apply( /* Inject energy into the particle */ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new); - hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new); + hydro_set_drifted_physical_internal_energy(pj, cosmo, /*pfloor=*/NULL, + u_new); /* Impose maximal viscosity */ hydro_diffusive_feedback_reset(pj); diff --git a/src/feedback/GEAR/feedback.c b/src/feedback/GEAR/feedback.c index a97fc832418f0c36fd12574c931696015427e1c4..ba63dbf4d9c83fa6594a47b3449896e88a1d0fa0 100644 --- a/src/feedback/GEAR/feedback.c +++ b/src/feedback/GEAR/feedback.c @@ -46,6 +46,7 @@ void feedback_update_part(struct part* p, struct xpart* xp, if (xp->feedback_data.delta_mass == 0) return; const struct cosmology* cosmo = e->cosmology; + const struct pressure_floor_props* pressure_floor = e->pressure_floor_props; /* Turn off the cooling */ xp->cooling_data.time_last_event = e->time; @@ -71,7 +72,7 @@ void feedback_update_part(struct part* p, struct xpart* xp, const float u_new = u + xp->feedback_data.delta_u; hydro_set_physical_internal_energy(p, xp, cosmo, u_new); - hydro_set_drifted_physical_internal_energy(p, cosmo, u_new); + hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, u_new); xp->feedback_data.delta_u = 0.; diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 1a457a66b76dd601d5ca4ef3ac16324c22144c29..6f3269b623a0e2d442934253edebc7ae25449cce 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -44,6 +44,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include <float.h> @@ -624,7 +625,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { const float fac_B = cosmo->a_factor_Balsara_eps; @@ -782,7 +784,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { /* Here we need to update the artificial viscosity */ @@ -897,7 +900,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -935,7 +939,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Store ratio of new internal energy to old internal energy, as we use this * in the drifting of the pressure. */ @@ -1074,7 +1079,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 827d43c6d88b178d94c35844c5f970837c78cf59..3db1af8771c60a6e3370b83237672c590b67deb4 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -352,12 +352,13 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * * @param p The particle of interest. * @param cosmo Cosmology data structure + * @param pressure_floor The properties of the pressure floor. * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { p->entropy = gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, u); @@ -368,8 +369,8 @@ hydro_set_drifted_physical_internal_energy(struct part *p, /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = - pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); + comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = @@ -639,6 +640,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( * @param xp The extended particle data to act upon * @param cosmo The current cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. * @param dt_alpha The time-step used to evolve non-cosmological quantities such * as the artificial viscosity. * @param dt_therm The time-step used to evolve hydrodynamical quantities. @@ -646,7 +648,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps; @@ -667,8 +670,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = - pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); + comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = @@ -746,10 +749,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param p The particle. * @param xp The extended data of this particle. * @param cosmo The cosmological model. + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -761,8 +766,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = - pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); + comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, comoving_pressure, cosmo); /* Compute the new sound speed */ const float soundspeed = @@ -788,12 +793,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param cosmo The cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props The properties of the entropy floor. + * @param pressure_floor The properties of the pressure 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, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the entropy */ p->entropy += p->entropy_dt * dt_therm; @@ -832,8 +839,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( /* Re-compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = - pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); + comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, comoving_pressure, cosmo); /* Compute the new sound speed */ const float soundspeed = @@ -923,10 +930,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * @param xp The extended data. * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme. + * @param pressure_floor The properties of the pressure floor. */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* We read u in the entropy field. We now get (comoving) A from (physical) u * and (physical) rho. Note that comoving A (A') == physical A */ @@ -948,8 +957,8 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( /* Compute the pressure */ float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy); - comoving_pressure = - pressure_floor_get_comoving_pressure(p, comoving_pressure, cosmo); + comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, comoving_pressure, cosmo); /* Compute the sound speed */ const float soundspeed = diff --git a/src/hydro/Gasoline/hydro.h b/src/hydro/Gasoline/hydro.h index e9cad33611181eeb41aee1ce4d23a67807c23bfb..844b88db8e5de399c93c558f6d54a99e6a5979ab 100644 --- a/src/hydro/Gasoline/hydro.h +++ b/src/hydro/Gasoline/hydro.h @@ -347,12 +347,13 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * * @param p The particle of interest. * @param cosmo Cosmology data structure + * @param pressure_floor The properties of the pressure floor. * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { /* There is no need to use the floor here as this function is called in the * feedback, so the new value of the internal energy should be strictly * higher than the old value. */ @@ -364,7 +365,7 @@ hydro_set_drifted_physical_internal_energy(struct part *p, /* Compute the sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -593,14 +594,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( * @param xp The extended particle data to act upon. * @param cosmo The cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + /* Compute the sound speed */ const float pressure = hydro_get_comoving_pressure(p); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -753,6 +757,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( * @param xp The extended particle data to act upon * @param cosmo The current cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. * @param dt_alpha The time-step used to evolve non-cosmological quantities such * as the artificial viscosity. * @param dt_therm The time-step used to evolve hydrodynamical quantities. @@ -760,7 +765,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { /* Here we need to update the artificial viscosity alpha */ const float d_shock_indicator_dt = dt_alpha == 0.f ? 0.f @@ -839,10 +845,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param p The particle. * @param xp The extended data of this particle. * @param cosmo The cosmological model. + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; p->v[1] = xp->v_full[1]; @@ -854,7 +862,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Compute the sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -882,12 +890,15 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param cosmo The cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props The properties of the entropy floor. + * @param pressure_floor The properties of the pressure 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, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { + /* Predict the internal energy */ p->u += p->u_dt * dt_therm; @@ -926,7 +937,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( /* Compute the new sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -1011,10 +1022,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * @param xp The extended particle to act upon * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme. + * @param pressure_floor The properties of the pressure floor. */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* 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; @@ -1037,7 +1050,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h index 7e3ce039b0d0b5daaa14bc0e3789cdaffb19ec4e..698141905c9dcbe44a53065579308a7f7f0c3ed8 100644 --- a/src/hydro/Gizmo/hydro.h +++ b/src/hydro/Gizmo/hydro.h @@ -38,6 +38,7 @@ #include "hydro_setters.h" #include "hydro_space.h" #include "hydro_velocities.h" +#include "pressure_floor.h" #include <float.h> @@ -421,7 +422,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part* restrict p, struct xpart* restrict xp, - const struct cosmology* cosmo, const struct hydro_props* hydro_props) { + const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct pressure_floor_props* pressure_floor) { /* Initialize time step criterion variables */ p->timestepvars.vmax = 0.; @@ -483,7 +485,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient( __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, const float dt_therm) { + const struct pressure_floor_props* pressure_floor, const float dt_alpha, + const float dt_therm) { hydro_part_reset_gravity_fluxes(p); p->flux.dt = dt_therm; @@ -519,7 +522,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part* restrict p, const struct xpart* restrict xp, - const struct cosmology* cosmo) { + const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor) { // MATTHIEU: Apply the entropy floor here. } @@ -535,7 +539,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( */ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( struct part* p, struct xpart* xp, const struct cosmology* cosmo, - const struct hydro_props* hydro_props) { + const struct hydro_props* hydro_props, + const struct pressure_floor_props* pressure_floor) { p->conserved.energy /= cosmo->a_factor_internal_energy; } @@ -552,7 +557,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part* p, struct xpart* xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology* cosmo, const struct hydro_props* hydro_props, - const struct entropy_floor_properties* floor_props) { + const struct entropy_floor_properties* floor_props, + const struct pressure_floor_props* pressure_floor) { /* skip the drift if we are using Lloyd's algorithm */ hydro_gizmo_lloyd_skip_drift(); diff --git a/src/hydro/Gizmo/hydro_setters.h b/src/hydro/Gizmo/hydro_setters.h index 640912ccbb226a74f629afb0ac630079f6332854..2715f7baff09198891fec29dd42960167a24a916 100644 --- a/src/hydro/Gizmo/hydro_setters.h +++ b/src/hydro/Gizmo/hydro_setters.h @@ -19,6 +19,8 @@ #ifndef SWIFT_GIZMO_HYDRO_SETTERS_H #define SWIFT_GIZMO_HYDRO_SETTERS_H +#include "pressure_floor.h" + /** * @brief Set the primitive variables for the given particle to the given * values. @@ -291,9 +293,9 @@ hydro_set_physical_internal_energy(struct part* p, struct xpart* xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part* p, const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor, const float u) { hydro_set_comoving_internal_energy(p, u / cosmo->a_factor_internal_energy); } diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index 6b67431123a0b87b1189510682eb9587c6cb00a9..f50e76407700af0873db7d892e02ce58acf93ee3 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -44,6 +44,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" /** * @brief Returns the comoving internal energy of a particle at the last @@ -359,9 +360,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor const float u) { p->u = u / cosmo->a_factor_internal_energy; @@ -585,7 +586,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Resets the variables that are required for a gradient calculation. @@ -667,7 +669,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps; @@ -767,7 +770,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -812,7 +816,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the internal energy */ p->u += p->u_dt * dt_therm; @@ -937,7 +942,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ diff --git a/src/hydro/None/hydro.h b/src/hydro/None/hydro.h index 4e7dd9ad6b6ad96e0f5170a93c3733e836fe70c1..e883b321e79852699049f8adffdf23d28c26e15b 100644 --- a/src/hydro/None/hydro.h +++ b/src/hydro/None/hydro.h @@ -35,6 +35,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include <float.h> @@ -382,9 +383,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { error("Empty implementation"); } @@ -521,7 +522,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Resets the variables that are required for a gradient calculation. @@ -582,7 +584,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) {} + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) {} /** * @brief Reset acceleration fields of a particle @@ -605,7 +608,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) {} + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Predict additional particle fields forward in time when drifting @@ -629,7 +633,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) {} + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Finishes the force calculation. @@ -684,7 +689,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Initialises the particles for the first time diff --git a/src/hydro/Phantom/hydro.h b/src/hydro/Phantom/hydro.h index 4259bd8bacaf558ffa4b6fc89aa0ddfeeb228cba..66b7e047b5f15c53c2f0902f1403ce3e9c1919a7 100644 --- a/src/hydro/Phantom/hydro.h +++ b/src/hydro/Phantom/hydro.h @@ -42,6 +42,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include <float.h> @@ -375,9 +376,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { p->u = u / cosmo->a_factor_internal_energy; @@ -599,7 +600,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { const float fac_B = cosmo->a_factor_Balsara_eps; @@ -747,7 +749,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { /* Here we need to update the artificial viscosity */ @@ -835,7 +838,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -877,7 +881,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the internal energy */ p->u += p->u_dt * dt_therm; @@ -1003,7 +1008,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index daa9c39b19c6225ed96bf8daedda9229341a5c73..bb43b34dac6b029848077056784d35de0e465841 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -46,6 +46,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" /* * Note: Define PLANETARY_SPH_NO_BALSARA to disable the Balsara (1995) switch @@ -383,9 +384,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { p->u = u / cosmo->a_factor_internal_energy; @@ -586,7 +587,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Resets the variables that are required for a gradient calculation. @@ -665,7 +667,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps; @@ -774,7 +777,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -820,7 +824,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the internal energy */ p->u += p->u_dt * dt_therm; @@ -932,7 +937,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Compute the pressure */ const float pressure = diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 143f1abc88456857f9ebdaa30a527914476b3399..5dc00b8615f7b49e1ff9eef21be33eac1928d996 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -216,14 +216,16 @@ hydro_get_drifted_physical_entropy(const struct part *restrict p, * * @param p The particle of interest. * @param cosmo The cosmological model. + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void hydro_update_soundspeed( - struct part *restrict p, const struct cosmology *cosmo) { + struct part *restrict p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Compute the sound speed -- see theory section for justification */ /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */ - const float comoving_pressure = - pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); + const float comoving_pressure = pressure_floor_get_comoving_pressure( + p, pressure_floor, p->pressure_bar, cosmo); p->force.soundspeed = gas_soundspeed_from_pressure(p->rho, comoving_pressure); /* Also update the signal velocity; this could be a huge change! */ @@ -399,11 +401,12 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @param p The particle of interest. * @param cosmo Cosmology data structure * @param u The physical internal energy + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void -hydro_set_drifted_physical_internal_energy(struct part *p, - const struct cosmology *cosmo, - const float u) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { /* Store ratio of new internal energy to old internal energy, as we use this * in the drifting of the pressure. */ @@ -419,10 +422,11 @@ hydro_set_drifted_physical_internal_energy(struct part *p, p->pressure_bar *= internal_energy_ratio; /* Update variables. */ - hydro_update_soundspeed(p, cosmo); + hydro_update_soundspeed(p, cosmo, pressure_floor); const float comoving_pressure_with_floor = - pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, p->pressure_bar, + cosmo); p->force.pressure_bar_with_floor = comoving_pressure_with_floor; } @@ -700,6 +704,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( * @param xp The extended particle data to act upon * @param cosmo The current cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. * @param dt_alpha The time-step used to evolve non-cosmological quantities such * as the artificial viscosity. * @param dt_therm The time-step used to evolve hydrodynamical quantities. @@ -707,7 +712,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_B = cosmo->a_factor_Balsara_eps; @@ -720,7 +726,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float abs_div_v = fabsf(p->density.div_v); /* Compute the sound speed -- see theory section for justification */ - hydro_update_soundspeed(p, cosmo); + hydro_update_soundspeed(p, cosmo, pressure_floor); const float soundspeed = hydro_get_comoving_soundspeed(p); /* Compute the Balsara switch */ @@ -760,7 +766,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Get the pressures */ const float comoving_pressure_with_floor = - pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, p->pressure_bar, + cosmo); /* Update variables. */ p->force.f = grad_h_term; @@ -797,10 +804,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param p The particle. * @param xp The extended data of this particle. * @param cosmo The cosmological model. + * @param pressure_floor The properties of the pressure floor. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -811,7 +820,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( p->u = xp->u_full; /* Compute the sound speed */ - hydro_update_soundspeed(p, cosmo); + hydro_update_soundspeed(p, cosmo, pressure_floor); } /** @@ -831,12 +840,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param cosmo The cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props The properties of the entropy floor. + * @param pressure_floor The properties of the pressure 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, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Store ratio of new internal energy to old internal energy, as we use this * in the drifting of the pressure. */ @@ -888,11 +899,12 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->pressure_bar *= internal_energy_ratio; /* Compute the new sound speed */ - hydro_update_soundspeed(p, cosmo); + hydro_update_soundspeed(p, cosmo, pressure_floor); /* update the required variables */ const float comoving_pressure_with_floor = - pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, p->pressure_bar, + cosmo); p->force.pressure_bar_with_floor = comoving_pressure_with_floor; } @@ -972,10 +984,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * @param xp The extended particle to act upon * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme. + * @param pressure_floor The properties of the pressure floor. */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index 2829053b32da9c25d0956d0de55c6b45355ab0a7..69170dc8d793c7e7471c4105d2e1fb5da8646d91 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -48,6 +48,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" #include <float.h> @@ -380,9 +381,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { /* Store ratio of new internal energy to old internal energy, as we use this * in the drifting of the pressure. */ @@ -689,7 +690,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_B = cosmo->a_factor_Balsara_eps; @@ -795,7 +797,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -833,7 +836,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Store ratio of new internal energy to old internal energy, as we use this * in the drifting of the pressure. */ @@ -971,7 +975,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index 3939c59bc8fe2c418e734521ef3bb2e73f862269..312cd353d0afac1c05f00df7e91592314e4f6abc 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -42,6 +42,7 @@ #include "hydro_space.h" #include "kernel_hydro.h" #include "minmax.h" +#include "pressure_floor.h" /** * @brief Returns the comoving internal energy of a particle at the last @@ -176,9 +177,9 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * @param u The physical entropy */ __attribute__((always_inline)) INLINE static void -hydro_set_drifted_physical_internal_energy(struct part *p, - const struct cosmology *cosmo, - const float u) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { error("To be implemented"); } @@ -544,7 +545,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) {} /** * @brief Resets the variables that are required for a gradient calculation. @@ -624,7 +626,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { const float fac_mu = cosmo->a_factor_mu; @@ -716,7 +719,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -758,7 +762,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the entropy */ p->entropy += p->entropy_dt * dt_therm; @@ -881,7 +886,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* We read u in the entropy field. We now get S from u */ xp->entropy_full = diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h index cd65b60499d5f6071b4433d886b37412bf5280c5..f04c95abb3665bc6f20b7fc9d70f96c3ba32ebdd 100644 --- a/src/hydro/SPHENIX/hydro.h +++ b/src/hydro/SPHENIX/hydro.h @@ -369,12 +369,13 @@ hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, * * @param p The particle of interest. * @param cosmo Cosmology data structure + * @param pressure_floor The #pressure_floor_props used. * @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) { +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { /* There is no need to use the floor here as this function is called in the * feedback, so the new value of the internal energy should be strictly @@ -387,7 +388,7 @@ hydro_set_drifted_physical_internal_energy(struct part *p, /* Compute the sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -648,10 +649,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( * @param xp The extended particle data to act upon. * @param cosmo The cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The #pressure_floor_props used. */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part *restrict p, struct xpart *restrict xp, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { const float fac_B = cosmo->a_factor_Balsara_eps; @@ -666,7 +669,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( /* Compute the sound speed */ const float pressure = hydro_get_comoving_pressure(p); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -812,6 +815,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( * @param xp The extended particle data to act upon * @param cosmo The current cosmological model. * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The #pressure_floor_props used. * @param dt_alpha The time-step used to evolve non-cosmological quantities such * as the artificial viscosity. * @param dt_therm The time-step used to evolve hydrodynamical quantities. @@ -819,7 +823,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { /* Here we need to update the artificial viscosity */ @@ -829,7 +834,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed; const float pressure = hydro_get_comoving_pressure(p); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed_physical = gas_soundspeed_from_pressure(p->rho, pressure_including_floor) * cosmo->a_factor_sound_speed; @@ -956,10 +961,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param p The particle. * @param xp The extended data of this particle. * @param cosmo The cosmological model. + * @param pressure_floor The #pressure_floor_props used. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part *restrict p, const struct xpart *restrict xp, - const struct cosmology *cosmo) { + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { /* Re-set the predicted velocities */ p->v[0] = xp->v_full[0]; @@ -972,7 +979,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Compute the sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -1000,12 +1007,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param cosmo The cosmological model. * @param hydro_props The properties of the hydro scheme. * @param floor_props The properties of the entropy floor. + * @param pressure_floor The properties of the pressure 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, float dt_kick_grav, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *floor_props) { + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { /* Predict the internal energy */ p->u += p->u_dt * dt_therm; @@ -1045,7 +1054,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( /* Compute the new sound speed */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); @@ -1132,10 +1141,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * @param xp The extended particle to act upon * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme. + * @param pressure_floor The properties of the pressure floor. */ __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) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { /* Convert the physcial internal energy to the comoving one. */ /* u' = a^(3(g-1)) u */ @@ -1165,7 +1176,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure_including_floor = - pressure_floor_get_comoving_pressure(p, pressure, cosmo); + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure_including_floor); diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h index a4f88d08e97a86461f975ef7c8a9a47db20db9c7..2343d7706103b598a186aa1a092c423d1156ecc5 100644 --- a/src/hydro/Shadowswift/hydro.h +++ b/src/hydro/Shadowswift/hydro.h @@ -27,6 +27,7 @@ #include "hydro_gradients.h" #include "hydro_properties.h" #include "hydro_space.h" +#include "pressure_floor.h" #include "voronoi_algorithm.h" #include <float.h> @@ -307,7 +308,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( __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, const float dt_therm) { + const struct pressure_floor_props* pressure_floor, const float dt_alpha, + const float dt_therm) { /* Initialize time step criterion variables */ p->timestepvars.vmax = 0.0f; @@ -341,7 +343,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( */ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( struct part* restrict p, struct xpart* restrict xp, - const struct cosmology* cosmo, const struct hydro_props* hydro_props) { + const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct pressure_floor_props* pressure_floor) { /* Initialize time step criterion variables */ p->timestepvars.vmax = 0.; @@ -403,7 +406,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( struct part* restrict p, const struct xpart* restrict xp, - const struct cosmology* cosmo) {} + const struct cosmology* cosmo, + const struct pressure_floor_props* pressure_floor) {} /** * @brief Converts the hydrodynamic variables from the initial condition file to @@ -424,7 +428,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( */ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( struct part* p, struct xpart* xp, const struct cosmology* cosmo, - const struct hydro_props* hydro_props) {} + const struct hydro_props* hydro_props, + const struct pressure_floor_props* pressure_floor) {} /** * @brief Extra operations to be done during the drift @@ -437,8 +442,10 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part* p, struct xpart* xp, float dt_drift, float dt_therm, - const struct cosmology* cosmo, const struct hydro_props* hydro_props, - const struct entropy_floor_properties* floor_props) {} + float dt_kick_grav, const struct cosmology* cosmo, + const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props, + const struct pressure_floor_props* pressure_floor) {} /** * @brief Set the particle acceleration after the flux loop. @@ -458,9 +465,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt Physical time step. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( - struct part* p, struct xpart* xp, float dt, float dt_grav, float dt_hydro, - float dt_kick_corr, const struct cosmology* cosmo, - const struct hydro_props* hydro_props, + struct part* p, struct xpart* xp, float dt, float dt_grav, + float dt_grav_mesh, float dt_hydro, float dt_kick_corr, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, const struct entropy_floor_properties* floor_props) { /* Update the conserved variables. We do this here and not in the kick, @@ -624,17 +631,6 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities( v[2] = p->v[2]; } -/** - * @brief Returns the density of a particle - * - * @param p The particle of interest - */ -__attribute__((always_inline)) INLINE static float hydro_get_density( - const struct part* restrict p) { - - return p->primitives.rho; -} - /** * @brief Modifies the thermal state of a particle to the imposed internal * energy diff --git a/src/lightcone/lightcone_particle_io.c b/src/lightcone/lightcone_particle_io.c index 2d98f46d20996ad78dbbd73c882125e08f985869..b58feec5c42569ee082d4b7939bb89c7f2bcc788 100644 --- a/src/lightcone/lightcone_particle_io.c +++ b/src/lightcone/lightcone_particle_io.c @@ -406,7 +406,7 @@ int lightcone_store_gas(const struct engine *e, struct lightcone_props *props, data->mass = hydro_get_mass(p); data->a = a_cross; data->h = p->h; - data->rho = p->rho; + data->rho = hydro_get_comoving_density(p); data->temperature = cooling_get_temperature( e->physical_constants, e->hydro_properties, e->internal_units, e->cosmology, e->cooling_func, p, xp); diff --git a/src/pressure_floor.h b/src/pressure_floor.h index 07f46c6e26434ca0b6acd061ba772c4cc019f2ed..6c4a62310310aa390520cb56bd009a942b31185b 100644 --- a/src/pressure_floor.h +++ b/src/pressure_floor.h @@ -34,13 +34,28 @@ /* Check if pressure floor is implemented in hydro */ #ifndef PRESSURE_FLOOR_NONE -#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH) || defined(SPHENIX_SPH) +#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH) || defined(SPHENIX_SPH) || \ + defined(GASOLINE_SPH) /* Implemented */ #else #error Pressure floor not implemented with this hydro scheme #endif +#endif + +/* Check if the pressure floor is implemented in the stellar feedback */ +#ifdef PRESSURE_FLOOR_GEAR +#if defined(FEEDBACK_EAGLE_THERMAL) || defined(FEEDBACK_EAGLE_KINETIC) +#error Pressure floor not implemented in this stellar feedback scheme +#endif +#endif +/* Check if the pressure floor is implemented in the AGN feedback */ +#ifdef PRESSURE_FLOOR_GEAR +#if defined(BLACK_HOLES_EAGLE) || defined(BLACK_HOLES_SPIN_JET) +#error Pressure floor not implemented in this AGN feedback scheme +#endif #endif + /* Import the right pressure floor definition */ #if defined(PRESSURE_FLOOR_NONE) #include "./pressure_floor/none/pressure_floor.h" diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h index e5e7083139dbe4894a2c40bbd46b9fd9dd0e94b1..362b831b81aa7b6339442f8194ab8c56e938dbc7 100644 --- a/src/pressure_floor/GEAR/pressure_floor.h +++ b/src/pressure_floor/GEAR/pressure_floor.h @@ -21,8 +21,12 @@ /* Forward declaration */ struct cosmology; +struct pressure_floor_props; + __attribute__((always_inline)) static INLINE float -pressure_floor_get_comoving_pressure(const struct part* p, const float pressure, +pressure_floor_get_comoving_pressure(const struct part* p, + const struct pressure_floor_props* floor, + const float pressure_comoving, const struct cosmology* cosmo); #include "adiabatic_index.h" @@ -58,13 +62,15 @@ struct pressure_floor_props { * Note that the particle is not updated!! * * @param p The #part. + * @param pflorr The properties of the pressure floor. * @param pressure_comoving The comoving pressure without any pressure floor. * @param cosmo The #cosmology model. * - * @return The comoving pressure with the floor. + * @return The comoving pressure with the floor applied. */ __attribute__((always_inline)) static INLINE float pressure_floor_get_comoving_pressure(const struct part* p, + const struct pressure_floor_props* pfloor, const float pressure_comoving, const struct cosmology* cosmo) { @@ -73,7 +79,7 @@ pressure_floor_get_comoving_pressure(const struct part* p, /* Compute the pressure floor */ float floor = kernel_gamma * kernel_gamma * p->h * p->h * rho * - pressure_floor_props.constants * cosmo->a_inv; + pfloor->constants * cosmo->a_inv; floor *= a_coef * rho * hydro_one_over_gamma; return fmaxf(pressure_comoving, floor); @@ -133,58 +139,6 @@ __attribute__((always_inline)) INLINE static void pressure_floor_print_snapshot( #endif -/** - * @brief Finishes the density calculation for the pressure floor properties. - * - * @param p The particle to act upon - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void pressure_floor_end_density( - struct part* restrict p, const struct cosmology* cosmo) {} - -/** - * @brief Sets all the pressure floor fields to sensible values when the #part - * has 0 ngbs. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void -pressure_floor_part_has_no_neighbours(struct part* restrict p, - struct xpart* restrict xp, - const struct cosmology* cosmo) {} - -/** - * @brief Sets the pressure_floor properties of the (x-)particles to a valid - * start state. - * - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void pressure_floor_init_part( - struct part* restrict p, struct xpart* restrict xp) {} - -/** - * @brief Sets the pressure_floor properties of the (x-)particles to a valid - * start state. - * - * @param phys_const The physical constant in internal units. - * @param us The unit system. - * @param cosmo The current cosmological model. - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void -pressure_floor_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - struct part* restrict p, - struct xpart* restrict xp) { - - pressure_floor_init_part(p, xp); -} - /** * @brief Write a pressure_floor struct to the given FILE as a stream of bytes. * @@ -218,10 +172,6 @@ __attribute__((always_inline)) INLINE static void pressure_floor_struct_restore( message("restoring pressure_floor..."); pressure_floor_print(pressure_floor); - - /* copy to the global variable */ - pressure_floor_props.n_jeans = pressure_floor->n_jeans; - pressure_floor_props.constants = pressure_floor->constants; } #endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */ diff --git a/src/pressure_floor/none/pressure_floor.h b/src/pressure_floor/none/pressure_floor.h index 148a8303d81cb5be766a34585bd7ebe943e87b46..3d1cb3afb4a5c05fbdbe07bba807532886a78580 100644 --- a/src/pressure_floor/none/pressure_floor.h +++ b/src/pressure_floor/none/pressure_floor.h @@ -38,26 +38,6 @@ struct unit_system; */ struct pressure_floor_props {}; -/** - * @brief Compute the physical pressure floor of a given #part. - * - * Note that the particle is not updated!! - * - * Since this is the 'none' model for floor, there is no floor and - * we just return the physical pressure that was received. - * - * @param p The #part. - * @param physical_pressure The physical pressure without any pressure floor. - * @param cosmo The #cosmology model. - * - * @return The physical pressure with the floor. - */ -static INLINE float pressure_floor_get_physical_pressure( - const struct part* p, const float physical_pressure, - const struct cosmology* cosmo) { - return physical_pressure; -} - /** * @brief Compute the comoving pressure floor of a given #part. * @@ -67,14 +47,15 @@ static INLINE float pressure_floor_get_physical_pressure( * we just return the comoving pressure that was received. * * @param p The #part. + * @param pfloor The pressure floor properties. * @param comoving_pressure The comoving pressure without any pressure floor. * @param cosmo The #cosmology model. * * @return The comoving pressure with the floor. */ static INLINE float pressure_floor_get_comoving_pressure( - const struct part* p, const float comoving_pressure, - const struct cosmology* cosmo) { + const struct part* p, const struct pressure_floor_props* pfloor, + const float comoving_pressure, const struct cosmology* cosmo) { return comoving_pressure; } @@ -119,63 +100,6 @@ INLINE static void pressure_floor_print_snapshot(hid_t h_grp) { } #endif -/** - * @brief Finishes the density calculation for the pressure floor properties. - * - * Nothing to do here. - * - * @param p The particle to act upon - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void pressure_floor_end_density( - struct part* restrict p, const struct cosmology* cosmo) {} - -/** - * @brief Sets all the pressure floor fields to sensible values when the #part - * has 0 ngbs. - * - * Nothing to do here. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - * @param cosmo The current cosmological model. - */ -__attribute__((always_inline)) INLINE static void -pressure_floor_part_has_no_neighbours(struct part* restrict p, - struct xpart* restrict xp, - const struct cosmology* cosmo) {} - -/** - * @brief Sets the pressure_floor properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void pressure_floor_init_part( - struct part* restrict p, struct xpart* restrict xp) {} - -/** - * @brief Sets the pressure_floor properties of the (x-)particles to a valid - * start state. - * - * Nothing to do here. - * - * @param phys_const The physical constant in internal units. - * @param us The unit system. - * @param cosmo The current cosmological model. - * @param p Pointer to the particle data. - * @param xp Pointer to the extended particle data. - */ -__attribute__((always_inline)) INLINE static void -pressure_floor_first_init_part(const struct phys_const* restrict phys_const, - const struct unit_system* restrict us, - const struct cosmology* restrict cosmo, - struct part* restrict p, - struct xpart* restrict xp) {} - /** * @brief Write a pressure_floor struct to the given FILE as a stream of bytes. * diff --git a/src/runner_ghost.c b/src/runner_ghost.c index c5ba6ecc21fa2a1696c80ee38efd505568d58f0a..48c8055b8cfb8287d649b0036da4af67c01eba85 100644 --- a/src/runner_ghost.c +++ b/src/runner_ghost.c @@ -32,8 +32,6 @@ #include "engine.h" #include "feedback.h" #include "mhd.h" -#include "pressure_floor.h" -#include "pressure_floor_iact.h" #include "rt.h" #include "space_getsid.h" #include "star_formation.h" @@ -1004,6 +1002,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { const double time_base = e->time_base; const struct cosmology *cosmo = e->cosmology; const struct hydro_props *hydro_props = e->hydro_properties; + const struct pressure_floor_props *pressure_floor = e->pressure_floor_props; TIMER_TIC; @@ -1051,7 +1050,8 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) { } /* Compute variables required for the force loop */ - hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha, dt_therm); + hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, dt_alpha, + dt_therm); mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha); timestep_limiter_prepare_force(p, xp); rt_prepare_force(p); @@ -1093,6 +1093,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const struct chemistry_global_data *chemistry = e->chemistry; const struct star_formation *star_formation = e->star_formation; const struct hydro_props *hydro_props = e->hydro_properties; + const struct pressure_floor_props *pressure_floor = e->pressure_floor_props; const int with_cosmology = (e->policy & engine_policy_cosmology); const int with_rt = (e->policy & engine_policy_rt); @@ -1200,7 +1201,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { hydro_end_density(p, cosmo); mhd_end_density(p, cosmo); chemistry_end_density(p, chemistry, cosmo); - pressure_floor_end_density(p, cosmo); star_formation_end_density(p, xp, star_formation, cosmo); /* Are we using the alternative definition of the @@ -1252,7 +1252,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* The force variables are set in the extra ghost. */ /* Compute variables required for the gradient loop */ - hydro_prepare_gradient(p, xp, cosmo, hydro_props); + hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor); mhd_prepare_gradient(p, xp, cosmo, hydro_props); /* The particle gradient values are now set. Do _NOT_ @@ -1288,7 +1288,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* As of here, particle force variables will be set. */ /* Compute variables required for the force loop */ - hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha, dt_therm); + hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, + dt_alpha, dt_therm); mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha); timestep_limiter_prepare_force(p, xp); rt_prepare_force(p); @@ -1380,7 +1381,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { hydro_init_part(p, hs); mhd_init_part(p); chemistry_init_part(p, chemistry); - pressure_floor_init_part(p, xp); star_formation_init_part(p, star_formation); tracers_after_init(p, xp, e->internal_units, e->physical_constants, with_cosmology, e->cosmology, @@ -1406,7 +1406,6 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { hydro_part_has_no_neighbours(p, xp, cosmo); mhd_part_has_no_neighbours(p, xp, cosmo); chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo); - pressure_floor_part_has_no_neighbours(p, xp, cosmo); star_formation_part_has_no_neighbours(p, xp, star_formation, cosmo); rt_part_has_no_neighbours(p); @@ -1433,7 +1432,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* The force variables are set in the extra ghost. */ /* Compute variables required for the gradient loop */ - hydro_prepare_gradient(p, xp, cosmo, hydro_props); + hydro_prepare_gradient(p, xp, cosmo, hydro_props, pressure_floor); mhd_prepare_gradient(p, xp, cosmo, hydro_props); /* The particle gradient values are now set. Do _NOT_ @@ -1469,7 +1468,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { /* As of here, particle force variables will be set. */ /* Compute variables required for the force loop */ - hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha, dt_therm); + hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, dt_alpha, + dt_therm); mhd_prepare_force(p, xp, cosmo, hydro_props, dt_alpha); timestep_limiter_prepare_force(p, xp); rt_prepare_force(p); diff --git a/src/runner_others.c b/src/runner_others.c index 35aead47ba64e7f84f6faddc70e9f42987952498..2112a7a2b4717fe3bd4661a8f7b7089c8b3637a7 100644 --- a/src/runner_others.c +++ b/src/runner_others.c @@ -127,6 +127,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { const struct unit_system *us = e->internal_units; const struct hydro_props *hydro_props = e->hydro_properties; const struct entropy_floor_properties *entropy_floor_props = e->entropy_floor; + const struct pressure_floor_props *pressure_floor = e->pressure_floor_props; const double time_base = e->time_base; const integertime_t ti_current = e->ti_current; struct part *restrict parts = c->hydro.parts; @@ -175,8 +176,8 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { /* Let's cool ! */ cooling_cool_part(constants, us, cosmo, hydro_props, - entropy_floor_props, cooling_func, p, xp, dt_cool, - dt_therm, time); + entropy_floor_props, pressure_floor, cooling_func, p, + xp, dt_cool, dt_therm, time); } } } diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c index 7b31072f34fa9407da2b2c282131354b96e9d1df..80ba9e2912339d367127ed450cd2f32442fc18e1 100644 --- a/src/runner_time_integration.c +++ b/src/runner_time_integration.c @@ -361,6 +361,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, const int timer) { const struct engine *e = r->e; const struct cosmology *cosmo = e->cosmology; const struct hydro_props *hydro_props = e->hydro_properties; + const struct pressure_floor_props *pressure_floor = e->pressure_floor_props; const struct entropy_floor_properties *entropy_floor = e->entropy_floor; const int with_cosmology = (e->policy & engine_policy_cosmology); const int periodic = e->s->periodic; @@ -459,7 +460,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, const int timer) { #endif /* Prepare the values to be drifted */ - hydro_reset_predicted_values(p, xp, cosmo); + hydro_reset_predicted_values(p, xp, cosmo, pressure_floor); } } diff --git a/src/space.c b/src/space.c index d955654f9afbe212d532f3076d5a00e5cc8104d2..4c2d21087f3d8f563484ce0071a669142cd4cc32 100644 --- a/src/space.c +++ b/src/space.c @@ -762,6 +762,7 @@ void space_convert_quantities_mapper(void *restrict map_data, int count, struct space *s = (struct space *)extra_data; const struct cosmology *cosmo = s->e->cosmology; const struct hydro_props *hydro_props = s->e->hydro_properties; + const struct pressure_floor_props *floor = s->e->pressure_floor_props; struct part *restrict parts = (struct part *)map_data; const ptrdiff_t index = parts - s->parts; struct xpart *restrict xparts = s->xparts + index; @@ -770,7 +771,8 @@ void space_convert_quantities_mapper(void *restrict map_data, int count, * creation */ for (int k = 0; k < count; k++) { if (parts[k].time_bin <= num_time_bins) { - hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props); + hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props, + floor); mhd_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props); } } diff --git a/src/space_first_init.c b/src/space_first_init.c index 71dbd26b50d9e9c2e6a19603e6d19ea4f22d82b5..ad98957bfefed0ea209296e7b542fc89ac1b5123 100644 --- a/src/space_first_init.c +++ b/src/space_first_init.c @@ -34,7 +34,6 @@ #include "mhd.h" #include "neutrino.h" #include "particle_splitting.h" -#include "pressure_floor.h" #include "rt.h" #include "sink.h" #include "star_formation.h" @@ -123,9 +122,6 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, /* Also initialise the chemistry */ chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k], &xp[k]); - /* Also initialise the pressure floor */ - pressure_floor_first_init_part(phys_const, us, cosmo, &p[k], &xp[k]); - /* Also initialise the star formation */ star_formation_first_init_part(phys_const, us, cosmo, star_formation, &p[k], &xp[k]); diff --git a/src/space_init.c b/src/space_init.c index c59ef1d2406deafa1284e8d090ad659e26a3ed80..900e841ba6e38065f470a12935c730fd2d65c8f8 100644 --- a/src/space_init.c +++ b/src/space_init.c @@ -31,7 +31,6 @@ #include "engine.h" #include "gravity.h" #include "mhd.h" -#include "pressure_floor.h" #include "rt.h" #include "sink.h" #include "star_formation.h" @@ -55,7 +54,6 @@ void space_init_parts_mapper(void *restrict map_data, int count, mhd_init_part(&parts[k]); black_holes_init_potential(&parts[k].black_holes_data); chemistry_init_part(&parts[k], e->chemistry); - pressure_floor_init_part(&parts[k], &xparts[k]); rt_init_part(&parts[k]); rt_reset_part(&parts[k], e->cosmology); star_formation_init_part(&parts[k], e->star_formation); diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index 2bf63d0e732613f5b0ba38680e96275e09e342b0..66337734fac502d86dcdb298360f7eae19f90532 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -411,9 +411,20 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density( /* Copy the velocity divergence */ xp->sf_data.div_v = p->density.div_v; xp->sf_data.div_v += hydro_dimension * cosmo->H; +#elif MINIMAL_SPH + /* Copy the velocity divergence */ + xp->sf_data.div_v = p->density.div_v; + xp->sf_data.div_v += hydro_dimension * cosmo->H; +#elif GASOLINE_SPH + /* Copy the velocity divergence */ + xp->sf_data.div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] + + p->viscosity.velocity_gradient[1][1] + + p->viscosity.velocity_gradient[2][2]); +#elif HOPKINS_PU_SPH + xp->sf_data.div_v = p->density.div_v; #else -#error This scheme is not implemented. Different scheme apply the Hubble flow \ - at different place. Be careful about it. +#error \ + "This scheme is not implemented. Note that Different scheme apply the Hubble flow in different places. Be careful about it." #endif } diff --git a/tests/test125cells.c b/tests/test125cells.c index 5716625aab829b403b90755038267785697dcec7..9d2c3638668d5714478fa3a36d95ed24793a235c 100644 --- a/tests/test125cells.c +++ b/tests/test125cells.c @@ -615,6 +615,9 @@ int main(int argc, char *argv[]) { lightcone_array_properties.nr_lightcones = 0; engine.lightcone_array_properties = &lightcone_array_properties; + struct pressure_floor_props pressure_floor; + engine.pressure_floor_props = &pressure_floor; + /* Construct some cells */ struct cell *cells[125]; struct cell *inner_cells[27]; diff --git a/tests/test27cells.c b/tests/test27cells.c index 24f9143ae4055500bec5c18bac1bc054a49ae8a6..b463a3b328b5733847ccafaab88a58ebbc4673e5 100644 --- a/tests/test27cells.c +++ b/tests/test27cells.c @@ -502,6 +502,9 @@ int main(int argc, char *argv[]) { lightcone_array_properties.nr_lightcones = 0; engine.lightcone_array_properties = &lightcone_array_properties; + struct pressure_floor_props pressure_floor; + engine.pressure_floor_props = &pressure_floor; + /* Construct some cells */ struct cell *cells[27]; struct cell *main_cell; diff --git a/tests/testActivePair.c b/tests/testActivePair.c index 4e27bbdacba0b1112ca9018109b718985cc5d051..ed02a6890d07321f81b2a8f6454b42177ad77f11 100644 --- a/tests/testActivePair.c +++ b/tests/testActivePair.c @@ -34,7 +34,8 @@ /* Typdef function pointer for interaction function. */ typedef void (*interaction_func)(struct runner *, struct cell *, struct cell *); typedef void (*init_func)(struct cell *, const struct cosmology *, - const struct hydro_props *); + const struct hydro_props *, + const struct pressure_floor_props *); typedef void (*finalise_func)(struct cell *, const struct cosmology *); /** @@ -180,8 +181,11 @@ void clean_up(struct cell *ci) { /** * @brief Initializes all particles field to be ready for a density calculation */ -void zero_particle_fields_density(struct cell *c, const struct cosmology *cosmo, - const struct hydro_props *hydro_props) { +void zero_particle_fields_density( + struct cell *c, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + for (int pid = 0; pid < c->hydro.count; pid++) { #if defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) c->hydro.parts[pid].geometry.wcorr = 1.0f; @@ -194,8 +198,11 @@ void zero_particle_fields_density(struct cell *c, const struct cosmology *cosmo, /** * @brief Initializes all particles field to be ready for a force calculation */ -void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo, - const struct hydro_props *hydro_props) { +void zero_particle_fields_force( + struct cell *c, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + for (int pid = 0; pid < c->hydro.count; pid++) { struct part *p = &c->hydro.parts[pid]; struct xpart *xp = &c->hydro.xparts[pid]; @@ -262,7 +269,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo, #endif /* And prepare for a round of force tasks. */ - hydro_prepare_force(p, xp, cosmo, hydro_props, 0., 0.); + hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, 0., 0.); hydro_reset_acceleration(p); } } @@ -345,19 +352,21 @@ void test_pair_interactions(struct runner *runner, struct cell **ci, interaction_func vec_interaction, init_func init, finalise_func finalise) { + const struct engine *e = runner->e; + runner_do_hydro_sort(runner, *ci, 0x1FFF, 0, 0, 0); runner_do_hydro_sort(runner, *cj, 0x1FFF, 0, 0, 0); /* Zero the fields */ - init(*ci, runner->e->cosmology, runner->e->hydro_properties); - init(*cj, runner->e->cosmology, runner->e->hydro_properties); + init(*ci, e->cosmology, e->hydro_properties, e->pressure_floor_props); + init(*cj, e->cosmology, e->hydro_properties, e->pressure_floor_props); /* Run the test */ vec_interaction(runner, *ci, *cj); /* Let's get physical ! */ - finalise(*ci, runner->e->cosmology); - finalise(*cj, runner->e->cosmology); + finalise(*ci, e->cosmology); + finalise(*cj, e->cosmology); /* Dump if necessary */ dump_particle_fields(swiftOutputFileName, *ci, *cj); @@ -365,15 +374,15 @@ void test_pair_interactions(struct runner *runner, struct cell **ci, /* Now perform a brute-force version for accuracy tests */ /* Zero the fields */ - init(*ci, runner->e->cosmology, runner->e->hydro_properties); - init(*cj, runner->e->cosmology, runner->e->hydro_properties); + init(*ci, e->cosmology, e->hydro_properties, e->pressure_floor_props); + init(*cj, e->cosmology, e->hydro_properties, e->pressure_floor_props); /* Run the brute-force test */ serial_interaction(runner, *ci, *cj); /* Let's get physical ! */ - finalise(*ci, runner->e->cosmology); - finalise(*cj, runner->e->cosmology); + finalise(*ci, e->cosmology); + finalise(*cj, e->cosmology); dump_particle_fields(bruteForceOutputFileName, *ci, *cj); } @@ -538,6 +547,7 @@ int main(int argc, char *argv[]) { struct engine engine; struct cosmology cosmo; struct hydro_props hydro_props; + struct pressure_floor_props pressure_floor; struct phys_const prog_const; struct runner *runner; static long long partId = 0; @@ -628,6 +638,7 @@ int main(int argc, char *argv[]) { engine.cosmology = &cosmo; hydro_props_init_no_hydro(&hydro_props); engine.hydro_properties = &hydro_props; + engine.pressure_floor_props = &pressure_floor; if (posix_memalign((void **)&runner, SWIFT_STRUCT_ALIGNMENT, sizeof(struct runner)) != 0) { diff --git a/tests/testComovingCooling.c b/tests/testComovingCooling.c index 9189e718844b89f731614f1e16ce59363c89e16d..56bbbc16df3c40ad456305476ecc2ab6f81b97d9 100644 --- a/tests/testComovingCooling.c +++ b/tests/testComovingCooling.c @@ -23,6 +23,8 @@ #if defined(CHEMISTRY_EAGLE) && defined(COOLING_EAGLE) && defined(GADGET2_SPH) +#include "../src/cooling/EAGLE/cooling_rates.h" + /* * @brief Assign particle density and entropy corresponding to the * hydrogen number density and internal energy specified. @@ -69,6 +71,9 @@ int main(int argc, char **argv) { struct part p; struct xpart xp; struct phys_const phys_const; + struct hydro_props hydro_properties; + struct entropy_floor_properties floor_props; + struct pressure_floor_props pressure_floor; struct cooling_function_data cooling; struct cosmology cosmo; char *parametersFileName = "./testCooling.yml"; @@ -109,22 +114,24 @@ int main(int argc, char **argv) { cosmology_init(params, &us, &phys_const, &cosmo); cosmology_print(&cosmo); + /* Init hydro_props */ + hydro_props_init(&hydro_properties, &phys_const, &us, params); + + /* Init entropy floor */ + entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params); + + /* Init the pressure floor */ + pressure_floor_init(&pressure_floor, &phys_const, &us, &hydro_properties, + params); + /* Set dt */ const int timebin = 38; float dt_cool, dt_therm; - /* Init hydro_props */ - struct hydro_props hydro_properties; - hydro_props_init(&hydro_properties, &phys_const, &us, params); - /* Init cooling */ cooling_init(params, &us, &phys_const, &hydro_properties, &cooling); cooling_print(&cooling); - cooling_update(&cosmo, &cooling, 0); - - /* Init entropy floor */ - struct entropy_floor_properties floor_props; - entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params); + cooling_update(&cosmo, &pressure_floor, &cooling, 0); /* Cooling function needs to know the minimal energy. Set it to the lowest * internal energy in the cooling table. */ @@ -178,11 +185,12 @@ int main(int argc, char **argv) { cosmology_get_therm_kick_factor(&cosmo, ti_begin, ti_begin + ti_step); cooling_init(params, &us, &phys_const, &hydro_properties, &cooling); - cooling_update(&cosmo, &cooling, 0); + cooling_update(&cosmo, &pressure_floor, &cooling, 0); /* compute implicit solution */ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, - &floor_props, &cooling, &p, &xp, dt_cool, dt_therm); + &floor_props, &pressure_floor, &cooling, &p, &xp, + dt_cool, dt_therm, 0); du_dt_check = hydro_get_physical_internal_energy_dt(&p, &cosmo); /* Now we can test the cooling at various redshifts and compare the result @@ -199,11 +207,12 @@ int main(int argc, char **argv) { /* Load the appropriate tables */ cooling_init(params, &us, &phys_const, &hydro_properties, &cooling); - cooling_update(&cosmo, &cooling, 0); + cooling_update(&cosmo, &pressure_floor, &cooling, 0); /* compute implicit solution */ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, - &floor_props, &cooling, &p, &xp, dt_cool, dt_therm); + &floor_props, &pressure_floor, &cooling, &p, &xp, + dt_cool, dt_therm, 0.); du_dt_implicit = hydro_get_physical_internal_energy_dt(&p, &cosmo); /* check if the two solutions are consistent */ diff --git a/tests/testCooling.c b/tests/testCooling.c index 62a63e4ee07c0616f009994203ea94865c380425..2785c9274d267aa6b8258807e3bd7db88afb79af 100644 --- a/tests/testCooling.c +++ b/tests/testCooling.c @@ -23,6 +23,8 @@ #if defined(CHEMISTRY_EAGLE) && defined(COOLING_EAGLE) && defined(GADGET2_SPH) +#include "../src/cooling/EAGLE/cooling_rates.h" + /* * @brief Assign particle density and entropy corresponding to the * hydrogen number density and internal energy specified. @@ -69,6 +71,9 @@ int main(int argc, char **argv) { struct part p; struct xpart xp; struct phys_const phys_const; + struct hydro_props hydro_properties; + struct entropy_floor_properties floor_props; + struct pressure_floor_props pressure_floor; struct cooling_function_data cooling; struct cosmology cosmo; char *parametersFileName = "./testCooling.yml"; @@ -110,22 +115,24 @@ int main(int argc, char **argv) { cosmology_init(params, &us, &phys_const, &cosmo); cosmology_print(&cosmo); + /* Init hydro_props */ + hydro_props_init(&hydro_properties, &phys_const, &us, params); + + /* Init entropy floor */ + entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params); + + /* Init the pressure floor */ + pressure_floor_init(&pressure_floor, &phys_const, &us, &hydro_properties, + params); + /* Set dt */ const int timebin = 38; float dt_cool, dt_therm; - /* Init hydro_props */ - struct hydro_props hydro_properties; - hydro_props_init(&hydro_properties, &phys_const, &us, params); - /* Init cooling */ cooling_init(params, &us, &phys_const, &hydro_properties, &cooling); cooling_print(&cooling); - cooling_update(&cosmo, &cooling, 0); - - /* Init entropy floor */ - struct entropy_floor_properties floor_props; - entropy_floor_init(&floor_props, &phys_const, &us, &hydro_properties, params); + cooling_update(&cosmo, &pressure_floor, &cooling, 0); /* Cooling function needs to know the minimal energy. Set it to the lowest * internal energy in the cooling table. */ @@ -169,7 +176,7 @@ int main(int argc, char **argv) { /* update nh, u, z */ cosmology_update(&cosmo, &phys_const, ti_current); cooling_init(params, &us, &phys_const, &hydro_properties, &cooling); - cooling_update(&cosmo, &cooling, 0); + cooling_update(&cosmo, &pressure_floor, &cooling, 0); set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh_cgs, u_cgs, ti_current); @@ -186,8 +193,8 @@ int main(int argc, char **argv) { for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) { p.entropy_dt = 0; cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, - &floor_props, &cooling, &p, &xp, - dt_cool / n_subcycle, dt_therm / n_subcycle); + &floor_props, &pressure_floor, &cooling, &p, &xp, + dt_cool / n_subcycle, dt_therm / n_subcycle, 0.); xp.entropy_full += p.entropy_dt * dt_therm / n_subcycle; } du_dt_check = hydro_get_physical_internal_energy_dt(&p, &cosmo); @@ -199,7 +206,8 @@ int main(int argc, char **argv) { /* compute implicit solution */ cooling_cool_part(&phys_const, &us, &cosmo, &hydro_properties, - &floor_props, &cooling, &p, &xp, dt_cool, dt_therm); + &floor_props, &pressure_floor, &cooling, &p, &xp, + dt_cool, dt_therm, 0.); du_dt_implicit = hydro_get_physical_internal_energy_dt(&p, &cosmo); /* check if the two solutions are consistent */ diff --git a/tests/testCooling.yml b/tests/testCooling.yml index fab3ae8959112a078e92621a0031ee14dbc2bc0d..fdd82a76d9115f927ede3d2c8bf5118c6aba30c3 100644 --- a/tests/testCooling.yml +++ b/tests/testCooling.yml @@ -11,7 +11,7 @@ Cosmology: h: 0.6777 # Reduced Hubble constant a_begin: 0.5 # Initial scale-factor of the simulation a_end: 1.0 # Final scale factor of the simulation - Omega_m: 0.307 # Matter density parameter + Omega_cdm: 0.2615 # Matter density parameter Omega_lambda: 0.693 # Dark-energy density parameter Omega_b: 0.0455 # Baryon density parameter @@ -89,19 +89,6 @@ EAGLECooling: He_reion_z_sigma: 0.5 He_reion_eV_p_H: 2.0 -#EAGLEChemistry: -# InitMetallicity: 0.014 -# InitAbundance_Hydrogen: 0.70649785 -# InitAbundance_Helium: 0.28055534 -# InitAbundance_Carbon: 2.0665436e-3 -# InitAbundance_Nitrogen: 8.3562563e-4 -# InitAbundance_Oxygen: 5.4926244e-3 -# InitAbundance_Neon: 1.4144605e-3 -# InitAbundance_Magnesium: 5.907064e-4 -# InitAbundance_Silicon: 6.825874e-4 -# InitAbundance_Iron: 1.1032152e-3 -# CalciumOverSilicon: 0.0941736 -# SulphurOverSilicon: 0.6054160 EAGLEChemistry: # Solar abundances init_abundance_metal: 0.014 init_abundance_Hydrogen: 0.70649785 @@ -114,10 +101,12 @@ EAGLEChemistry: # Solar abundances init_abundance_Silicon: 6.825874e-4 init_abundance_Iron: 1.1032152e-3 - GearChemistry: InitialMetallicity: 0.01295 +GEARPressureFloor: + jeans_factor: 10 + # Parameters for the EAGLE "equation of state" EAGLEEntropyFloor: Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3. diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c index ed50bf380a567d79d8b13b248cd36ccd124d788e..13c123f01f07f757a621b6fe1ac262313a2847bd 100644 --- a/tests/testPeriodicBC.c +++ b/tests/testPeriodicBC.c @@ -509,6 +509,9 @@ int main(int argc, char *argv[]) { lightcone_array_properties.nr_lightcones = 0; engine.lightcone_array_properties = &lightcone_array_properties; + struct pressure_floor_props pressure_floor; + engine.pressure_floor_props = &pressure_floor; + /* Construct some cells */ struct cell *cells[dim * dim * dim]; static long long partId = 0;