Fix GEAR pressure floor

.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

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#*_}

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

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

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

examples/GEAR/AgoraCosmo/getChemistryTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5

examples/GEAR/AgoraCosmo/getGrackleCoolingTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/GEAR/CloudyData_UVB=HM2012_shielded.h5
+
+

examples/GEAR/AgoraCosmo/getIC.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraCosmo/agora_swift.hdf5
+
+

examples/GEAR/AgoraCosmo/run.sh

 #!/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
+
+
+

examples/GEAR/AgoraDisk/README

+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/
 

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

examples/GEAR/AgoraDisk/changeType.py

-#!/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()

examples/GEAR/AgoraDisk/checkSolution.py

+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 * "!")

examples/GEAR/AgoraDisk/cleanupSwift.py

-#!/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()

examples/GEAR/AgoraDisk/getChemistryTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5

examples/GEAR/AgoraDisk/getGrackleCoolingTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/GEAR/AgoraDisk/getIC.sh

 #!/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
+

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

examples/GEAR/ZoomIn/getChemistryTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5

examples/GEAR/ZoomIn/getGrackleCoolingTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5