Merge branch 'gear_feedback' into 'master'

Gear

See merge request !987

configure.ac

@@ -1439,13 +1439,14 @@ case "$with_subgrid" in
    ;;
    GEAR)
 	with_subgrid_cooling=grackle_0
-	with_subgrid_chemistry=GEAR
+	with_subgrid_chemistry=GEAR_10
 	with_subgrid_pressure_floor=GEAR
 	with_subgrid_stars=GEAR
 	with_subgrid_star_formation=GEAR
-	with_subgrid_feedback=none
+	with_subgrid_feedback=GEAR
 	with_subgrid_black_holes=none
-	with_subgrid_task_order=GEAR
+  # GEAR's order is not used anymore
+	with_subgrid_task_order=none
 	enable_fof=no
    ;;
    EAGLE)
@@ -1791,7 +1792,8 @@ esac
 #  chemistry function
 AC_ARG_WITH([chemistry],
    [AS_HELP_STRING([--with-chemistry=<function>],
-      [chemistry function @<:@none, GEAR, EAGLE default: none@:>@]
+      [chemistry function @<:@none, GEAR_*, EAGLE default: none@:>@
+      For GEAR, you need to provide the number of elements (e.g. GEAR_10)]
    )],
    [with_chemistry="$withval"],
    [with_chemistry="none"]
@@ -1809,8 +1811,10 @@ case "$with_chemistry" in
    none)
       AC_DEFINE([CHEMISTRY_NONE], [1], [No chemistry function])
    ;;
-   GEAR)
+   GEAR_*)
       AC_DEFINE([CHEMISTRY_GEAR], [1], [Chemistry taken from the GEAR model])
+      number_element=${with_chemistry:5}
+      AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow])
    ;;
    EAGLE)
       AC_DEFINE([CHEMISTRY_EAGLE], [1], [Chemistry taken from the EAGLE model])
@@ -1913,6 +1917,9 @@ case "$with_feedback" in
    EAGLE)
       AC_DEFINE([FEEDBACK_EAGLE], [1], [EAGLE stellar feedback and evolution model])
    ;;
+   GEAR)
+      AC_DEFINE([FEEDBACK_GEAR], [1], [GEAR stellar feedback and evolution model])
+   ;;
    none)
       AC_DEFINE([FEEDBACK_NONE], [1], [No feedback])
    ;;
@@ -2139,6 +2146,12 @@ AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
 # Check if using EAGLE feedback
 AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"])
 
+# check if using grackle cooling
+AM_CONDITIONAL([HAVEGRACKLECOOLING], [test ${with_cooling:0:7} == "grackle"])
+
+# check if using gear feedback
+AM_CONDITIONAL([HAVEGEARFEEDBACK], [test $with_feedback == "GEAR"])
+
 # Handle .in files.
 AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
 AC_CONFIG_FILES([argparse/Makefile tools/Makefile logger/Makefile logger/tests/Makefile])

examples/GEAR/AgoraDisk/agora_disk.yml

@@ -7,15 +7,25 @@ InternalUnitSystem:
   UnitTemp_in_cgs:     1             # Kelvin
 
 Scheduler:
-  max_top_level_cells: 8
+  max_top_level_cells: 16
+  cell_extra_sparts:  5000       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_max_size:             8000   # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
+  cell_sub_size_pair_hydro:  2560 # (Optional) Maximal number of hydro-hydro interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_hydro:  3200     # (Optional) Maximal number of hydro-hydro interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_stars:  2560 # (Optional) Maximal number of hydro-star interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_stars:  3200     # (Optional) Maximal number of hydro-star interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_grav:   2560 # (Optional) Maximal number of interactions per sub-pair gravity task  (this is the default value).
+  cell_sub_size_self_grav:   3200     # (Optional) Maximal number of interactions per sub-self gravity task  (this is the default value).
+  cell_split_size:           100       # (Optional) Maximal number of particles per cell (this is the default value).
   
 # 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).
   dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-4  # The maximal 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
@@ -29,10 +39,12 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  eta:                    0.05    # Constant dimensionless multiplier for time integration.
   theta:                  0.7     # Opening angle (Multipole acceptance criterion)
-  comoving_softening:     0.08 # Comoving softening length (in internal units).
-  max_physical_softening: 0.08    # Physical softening length (in internal units).
+  comoving_DM_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_DM_softening: 0.08    # Physical softening length (in internal units).
+  comoving_baryon_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_baryon_softening: 0.08    # Physical softening length (in internal units).
   mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
   
 # Parameters for the hydrodynamics scheme
@@ -40,6 +52,8 @@ SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
   minimal_temperature:   10.      # Kelvin
+  h_min_ratio:           0.1095       # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+  h_max:                     10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
 
 # Parameters related to the initial conditions
 InitialConditions:
@@ -55,13 +69,34 @@ LambdaCooling:
 
 # Cooling with Grackle 2.0
 GrackleCooling:
-  CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
-  WithUVbackground: 1 # Enable or not the UV background
-  Redshift: 0 # Redshift to use (-1 means time based redshift)
-  WithMetalCooling: 1 # Enable or not the metal cooling
-  ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
-  ProvideSpecificHeatingRates: 0 # User provide specific heating rates
-  SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
-  OutputMode: 1 # Write in output corresponding primordial chemistry mode
-  MaxSteps: 1000
-  ConvergenceLimit: 1e-2
+  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1 # Enable or not the UV background
+  redshift: 0 # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1 # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
+  provide_specific_heating_rates: 0 # User provide specific heating rates
+  self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
+  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
+
+
+GEARStarFormation:
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  1e10         # Upper limit to the temperature of a star forming particle
+
+GEARPressureFloor:
+  jeans_factor: 10
+
+GEARFeedback:
+  supernovae_energy_erg: 0.1e51
+  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
+  discrete_yields: 0
+
+GEARChemistry:
+  initial_metallicity: 1
+  scale_initial_metallicity: 1
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.

examples/GEAR/AgoraDisk/cleanupSwift.py

@@ -14,7 +14,7 @@ copyfile(filename, out)
 f = File(out)
 
 for i in range(6):
-    name = "PartType{}/ElementAbundance".format(i)
+    name = "PartType{}/ElementAbundances".format(i)
     if name in f:
         del f[name]
 
@@ -23,13 +23,27 @@ for i in range(NPartType):
     if name not in f:
         continue
 
-    grp = f[name + "/SmoothingLength"]
+    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/getIC.sh

@@ -6,4 +6,4 @@ if [ "$#" -ne 1 ]; then
     exit
 fi
 
-wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/$1.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/$1.hdf5

examples/GEAR/AgoraDisk/getSolution.sh

 #!/bin/bash
 
-# cleanup work space
-rm snapshot_0000.hdf5 snapshot_0500.hdf5
-
-wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
-wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
-
-# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
-# wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
-
-
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0000.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0050.hdf5

examples/GEAR/AgoraDisk/plotSolution.py

@@ -34,7 +34,7 @@ from yt.analysis_modules.halo_analysis.api import *
 from yt.data_objects.particle_filters import add_particle_filter
 from scipy.stats import kde
 from subprocess import call
-#mylog.setLevel(1)
+# mylog.setLevel(1)
 
 from yt.utilities.math_utils import get_cyl_r, get_cyl_theta, get_cyl_z
 
@@ -53,7 +53,7 @@ draw_PDF                         = 2###2         # 0/1/2/3     = OFF/ON/ON with
 draw_pos_vel_PDF                 = 4##4          # 0/1/2/3/4   = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
 draw_star_pos_vel_PDF            = 4##4          # 0/1/2/3/4   = OFF/ON/ON with 1D profile/ON also with 1D dispersion profile/ON also with separate 1D vertical dispersion profiles
 draw_rad_height_PDF              = 2##2          # 0/1/2/3     = OFF/ON/ON with 1D profile/ON with analytic ftn subtracted
-draw_metal_PDF                   = 0##0          # 0/1         = OFF/ON
+draw_metal_PDF                   = 1##1          # 0/1         = OFF/ON
 draw_density_DF                  = 2#2           # 0/1/2       = OFF/ON/ON with difference plot between 1st and 2nd datasets (when 2, dataset_num should be set to 2)
 draw_radius_DF                   = 1#1           # 0/1         = OFF/ON
 draw_star_radius_DF              = 2#2           # 0/1/2       = OFF/ON/ON with SFR profile and K-S plot (when 2, this automatically turns on draw_radius_DF)
@@ -110,9 +110,9 @@ marker_names             = ['s', 'o', 'p', 'v', '^', '<', '>', 'h', '*']
 #             [file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
 codes = ['SWIFT', 'GEAR']
 filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
-              ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]],
+              ["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"]],
              [["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-SFF (with star formation and feedback)
-              ["./snapshot_0000.hdf5", "./snapshot_0500.hdf5"]]] # I did not check the order, they can be switched
+              ["./snapshot_0000.hdf5", "./snapshot_0050.hdf5"]]] # I did not check the order, they can be switched
 
 # codes = ["SWIFT"]
 # filenames = [[["./agora_disk_0000.hdf5", "./agora_disk_0050.hdf5"]],
@@ -139,8 +139,8 @@ filenames = [[["./agora_disk_IC.hdf5", "./agora_disk_500Myr.hdf5"], # Sim-noSFF
 # filenames = [[[file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_000.hdf5', file_location[0]+'GADGET-3/AGORA_ISO_LOW_DRY/snap_iso_dry_010.hdf5']],
 #            [[file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_000.hdf5', file_location[1]+'GADGET-3/AGORA_ISO_LOW_SF_SNII_Thermal_Chevalier_SFT10/snap_iso_sf_010.hdf5']]]
 # codes = ['GEAR']
-# filenames = [[['snapshot_0000', 'snapshot_0500']],
-#            [['snapshot_0000', 'snapshot_0500']]]
+# filenames = [[['snapshot_0000.hdf5', 'snapshot_0500.hdf5']],
+#            [['snapshot_0000.hdf5', 'snapshot_0500']]]
 # codes = ['GIZMO']
 # filenames = [[[file_location[0]+'GIZMO/snapshot_temp_000', file_location[0]+'GIZMO/snapshot_temp_100']],
 #            [[file_location[1]+'GIZMO/snapshot_temp_000', file_location[1]+'GIZMO/snapshot_temp_100']]]
@@ -439,7 +439,7 @@ for time in range(len(times)):
                                         yield line
         if draw_PDF >= 1:
                 fig_PDF              += [plt.figure(figsize=(50, 80))]
-                grid_PDF             += [AxesGrid(fig_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+                grid_PDF             += [AxesGrid(fig_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, int(math.ceil(len(codes)/2.0))), axes_pad = 0.05, add_all = True, share_all = True,
                                                   label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
         if draw_pos_vel_PDF >= 1:
                 fig_pos_vel_PDF      += [plt.figure(figsize=(50, 80))]
@@ -469,7 +469,7 @@ for time in range(len(times)):
                 rad_height_profiles.append([])
         if draw_metal_PDF == 1:
                 fig_metal_PDF        += [plt.figure(figsize=(50, 80))]
-                grid_metal_PDF       += [AxesGrid(fig_metal_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (3, int(math.ceil(len(codes)/3.0))), axes_pad = 0.05, add_all = True, share_all = True,
+                grid_metal_PDF       += [AxesGrid(fig_metal_PDF[time], (0.01,0.01,0.99,0.99), nrows_ncols = (2, int(math.ceil(len(codes)/2.0))), axes_pad = 0.05, add_all = True, share_all = True,
                                                   label_mode = "1", cbar_mode = "single", cbar_location = "right", cbar_size = "2%", cbar_pad = 0.05, aspect = False)]
         if draw_density_DF >= 1:
                 density_DF_xs.append([])
@@ -576,9 +576,9 @@ for time in range(len(times)):
                         PartType_StarBeforeFiltered_to_use = "PartType4"
                         if time != 0:
                             def _FormationTime(field, data):
-                                return pf.arr(data["PartType4", "BirthTime"].d, 'code_time')
+                                return pf.arr(data["PartType4", "BirthTimes"].d, 'code_time')
                             pf.add_field(("PartType4", FormationTimeType_to_use), function=_FormationTime, particle_type=True, take_log=False, units="code_time")
-                        
+
                         MassType_to_use = "Masses"
                 elif codes[code] == "GEAR":
                         PartType_Gas_to_use = "PartType0"
@@ -715,6 +715,20 @@ for time in range(len(times)):
                                 def _metallicity_3(field, data):
                                         return data["deposit", PartType_Gas_to_use+"_smoothed_"+MetallicityType_to_use]
                                 pf.add_field(("gas", "metallicity"), function=_metallicity_3, force_override=True, display_name="Metallicity", particle_type=False, take_log=True, units="")
+                        elif codes[code] == 'SWIFT': # "Metals" in SWIFT is 10-species field ([:,9] is the total metal fraction), so requires a change in _vector_fields in frontends/gadget/io.py: added ("Metals", 10)
+                                def _metallicity_2(field, data):
+                                        if len(data[PartType_Gas_to_use, "SmoothedElementAbundances"].shape) == 1:
+                                                return data[PartType_Gas_to_use, "SmoothedElementAbundances"]
+                                        else:
+                                                return data[PartType_Gas_to_use, "SmoothedElementAbundances"][:,9].in_units("") # in_units("") turned out to be crucial!; otherwise code_metallicity will be used and it will mess things up
+                                # We are creating ("Gas", "Metallicity") here, different from ("Gas", "metallicity") which is auto-generated by yt but doesn't work properly
+                                pf.add_field((PartType_Gas_to_use, MetallicityType_to_use), function=_metallicity_2, display_name="Metallicity", particle_type=True, take_log=True, units="")
+                                # Also creating smoothed field following an example in yt-project.org/docs/dev/cookbook/calculating_information.html; use hardcoded num_neighbors as in frontends/gadget/fields.py
+                                fn = add_volume_weighted_smoothed_field(PartType_Gas_to_use, "Coordinates", MassType_to_use, "SmoothingLength", "Density", MetallicityType_to_use, pf.field_info, nneighbors=64)
+                                # Alias doesn't work -- e.g. pf.field_info.alias(("gas", "metallicity"), fn[0]) -- probably because pf=GadgetDataset()?, not load()?; so I add and replace existing ("gas", "metallicity")
+                                def _metallicity_3(field, data):
+                                        return data["deposit", PartType_Gas_to_use+"_smoothed_"+MetallicityType_to_use]
+                                pf.add_field(("gas", "metallicity"), function=_metallicity_3, force_override=True, display_name="Metallicity", particle_type=False, take_log=True, units="")
                         if draw_metal_map >= 2:
                                 def _metal_mass(field, data):
                                         return data["gas", "cell_mass"] * data["gas", "metallicity"]
@@ -745,7 +759,7 @@ for time in range(len(times)):
                         if draw_metal_map >= 1 or draw_metal_PDF == 1:
                                 def _Metallicity_2(field, data):
                                         return data[(PartType_Gas_to_use, MetallicityType_to_use)]
-                                pf.add_field((PartType_Gas_to_use, "Metallicity_2"), function=_Metallicity_2, take_log=True, particle_type=False, display_name="Metallicity", units="")
+                                pf.add_field((PartType_Gas_to_use, "Metallicity_2"), function=_Metallicity_2, take_log=True, particle_type=True, display_name="Metallicity", units="")
                         def _Density_2_minus_analytic(field, data):
                                 return data[(PartType_Gas_to_use, "density")] - rho_agora_disk(data[(PartType_Gas_to_use, "radius")], data[(PartType_Gas_to_use, "particle_position_cylindrical_z_abs")])
                         pf.add_field((PartType_Gas_to_use, "Density_2_minus_analytic"), function=_Density_2_minus_analytic, take_log=False, particle_type=False, display_name="Density Residual", units="g/cm**3")
@@ -1565,13 +1579,18 @@ for time in range(len(times)):
                         else:
                                 # Because ParticlePhasePlot doesn't yet work for a log-log PDF for some reason, I will do the following trick.
                                 pf.field_info[(PartType_Gas_to_use, "Metallicity_2")].take_log = False
-                                p3 = PhasePlot(sp, (PartType_Gas_to_use, "density"), (PartType_Gas_to_use, "temperature"), (PartType_Gas_to_use, "Metallicity_2"), weight_field=(PartType_Gas_to_use, MassType_to_use), fontsize=12, x_bins=300, y_bins=300)
+                                p3 = ParticlePhasePlot(sp, (PartType_Gas_to_use, "density"), (PartType_Gas_to_use, "temperature"), (PartType_Gas_to_use, "Metallicity_2"),
+                                                       weight_field=(PartType_Gas_to_use, MassType_to_use), fontsize=12, x_bins=300, y_bins=300)
                                 p3.set_zlim((PartType_Gas_to_use, "Metallicity_2"), 0.01, 0.04)
                                 p3.set_cmap((PartType_Gas_to_use, "Metallicity_2"), my_cmap2)
                                 plot3 = p3.plots[(PartType_Gas_to_use, "Metallicity_2")]
 
+                        p3.set_unit("density", "g/cm**3")
                         p3.set_xlim(1e-29, 1e-21)
+                        p3.set_unit("temperature", "K")
                         p3.set_ylim(10, 1e7)
+                        p3.set_log("density", True)
+                        p3.set_log("temperature", True)
 
                         plot3.figure = fig_metal_PDF[time]
                         plot3.axes = grid_metal_PDF[time][code].axes
@@ -2271,7 +2290,7 @@ for time in range(len(times)):
         if draw_metal_PDF == 1:
                 fig_metal_PDF[time].savefig("metal_PDF_%dMyr" % times[time], bbox_inches='tight', pad_inches=0.03, dpi=300)
         if draw_density_DF >= 1:
-                plt.clf()
+#                 plt.clf()
 #		plt.subplot(111, aspect=1)
                 fig = plt.figure(figsize=(8, 8))
                 gridspec.GridSpec(4, 1)

examples/GEAR/AgoraDisk/run.sh

@@ -5,6 +5,8 @@
 # currently only the low resolution is available
 sim=low
 
+rm agora_disk_0*.hdf5
+
 # make run.sh fail if a subcommand fails
 set -e
 
@@ -19,7 +21,13 @@ fi
 if [ ! -e CloudyData_UVB=HM2012.h5 ]
 then
     echo "Fetching the Cloudy tables required by Grackle..."
-    ../../Cooling/getGrackleCoolingTable.sh 
+    ../../Cooling/getGrackleCoolingTable.sh
+fi
+
+if [ ! -e chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 ]
+then
+    echo "Fetching the chemistry tables..."
+    ../getChemistryTable.sh
 fi
 
 # copy the initial conditions
@@ -28,7 +36,7 @@ cp $sim.hdf5 agora_disk.hdf5
 python3 changeType.py agora_disk.hdf5
 
 # Run SWIFT
-../../swift --cooling --hydro --self-gravity --threads=4 agora_disk.yml 2>&1 | tee output.log
+../../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"

examples/GEAR/DwarfGalaxy/dwarf_galaxy.yml

@@ -68,3 +68,31 @@ InitialConditions:
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
 
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1                   # Enable or not the UV background
+  redshift: -1                           # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1                   # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0      # (optional) User provide volumetric heating rates
+  provide_specific_heating_rates: 0        # (optional) User provide specific heating rates
+  self_shielding_method: 0                # (optional) Grackle (<= 3) or Gear self shielding method
+  max_steps: 10000                       # (optional) Max number of step when computing the initial composition
+  convergence_limit: 1e-2                # (optional) Convergence threshold (relative) for initial composition
+  thermal_time_Myr:  5
+
+GearChemistry:
+  initial_metallicity: 0.01295
+
+GEARFeedback:
+  supernovae_energy_erg: 1e50
+  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
+
+GearPressureFloor:
+  jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+
+GEARStarFormation:
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle

examples/GEAR/DwarfGalaxy/run.sh

@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../../swift --feedback --self-gravity --hydro --stars --threads=8 $@ dwarf_galaxy.yml 2>&1 | tee output.log
+../../swift --feedback --limiter --sync --self-gravity --hydro --stars --cooling --star-formation --threads=8 $@ dwarf_galaxy.yml 2>&1 | tee output.log
 

examples/GEAR/ZoomIn/run.sh

@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../../swift --feedback --cosmology --self-gravity --hydro --stars --threads=8 zoom_in.yml 2>&1 | tee output.log
+../../swift --feedback --cosmology  --limiter --sync --self-gravity --hydro --stars --star-formation --threads=8 zoom_in.yml 2>&1 | tee output.log
 

examples/GEAR/ZoomIn/zoom_in.yml

@@ -60,3 +60,30 @@ InitialConditions:
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
 
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1                   # Enable or not the UV background
+  redshift: -1                           # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1                   # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0      # (optional) User provide volumetric heating rates
+  provide_specific_heating_rates: 0        # (optional) User provide specific heating rates
+  self_shielding_method: 0                # (optional) Grackle (<= 3) or Gear self shielding method
+  max_steps: 10000                       # (optional) Max number of step when computing the initial composition
+  convergence_limit: 1e-2                # (optional) Convergence threshold (relative) for initial composition
+  thermal_time_Myr:  5
+
+GearChemistry:
+  initial_metallicity: 0.01295
+
+GEARFeedback:
+  supernovae_energy_erg: 1e50
+  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
+
+GEARStarFormation:
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
+
+GEARPressureFloor:
+  jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.

examples/GEAR/getChemistryTable.sh

+#!/bin/bash
+wget https://obswww.unige.ch/lastro/projects/Clastro/PySWIFTsim/chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5

examples/SubgridTests/SupernovaeFeedback/SN_feedback.yml

-# Define the system of units to use internally. 
-InternalUnitSystem:
-  UnitMass_in_cgs:     1   # Grams
-  UnitLength_in_cgs:   1   # Centimeters
-  UnitVelocity_in_cgs: 1   # Centimeters per second
-  UnitCurrent_in_cgs:  1   # Amperes
-  UnitTemp_in_cgs:     1   # Kelvin
-
-# Values of some physical constants
-PhysicalConstants:
-  G:            0 # (Optional) Overwrite the value of Newton's constant used internally by the code.
-
-# Parameters governing the time integration
-TimeIntegration:
-  time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   5e-2  # The end time of the simulation (in internal units).
-  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
-  basename:            SN_feedback # 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:         1
- 
-# Parameters governing the conserved quantities statistics
-Statistics:
-  delta_time:          1e-3 # Time between statistics output
-
-# Parameters for the hydrodynamics scheme
-SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters related to the initial conditions
-InitialConditions:
-  file_name:                    ./SN_feedback.hdf5          
-  smoothing_length_scaling:     1.
-  periodic:                    1    # Are we running with periodic ICs?
- 
-# Parameters for the stellar models
-Stars:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).

examples/SubgridTests/SupernovaeFeedback/getGlass.sh

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

examples/SubgridTests/SupernovaeFeedback/makeIC.py

-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
-
-import h5py
-from numpy import *
-
-# Generates a swift IC file for the Sedov blast test in a periodic cubic box
-
-# Parameters
-gamma = 5./3.      # Gas adiabatic index
-rho0 = 1.          # Background density
-P0 = 1.e-6         # Background pressure
-E0= 1.             # Energy of the explosion
-N_inject = 15      # Number of particles in which to inject energy
-fileName = "SN_feedback.hdf5" 
-
-#---------------------------------------------------
-glass = h5py.File("glassCube_64.hdf5", "r")
-
-# Read particle positions and h from the glass
-pos = glass["/PartType0/Coordinates"][:,:]
-eps = 1e-6
-pos = (pos - pos.min()) / (pos.max() - pos.min() + eps)
-h = glass["/PartType0/SmoothingLength"][:] * 0.3 * 3.3
-
-numPart = size(h)
-vol = 1.
-Boxsize = 1.
-
-# Generate extra arrays
-v = zeros((numPart, 3))
-ids = linspace(1, numPart, numPart)
-m = zeros(numPart)
-u = zeros(numPart)
-r = zeros(numPart)
-
-r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2 + (pos[:,2] - 0.5)**2)
-m[:] = rho0 * vol / numPart    
-u[:] = P0 / (rho0 * (gamma - 1))
-
-#--------------------------------------------------
-
-star_pos = zeros((1, 3))
-star_pos[:,:] = 0.5 * Boxsize
-
-star_v = zeros((1, 3))
-star_v[:,:] = 0.
-
-# increase mass to keep it at center
-star_m = 1e3 * array([rho0 * vol / numPart])
-star_ids = array([numPart + 1])
-star_h = array([h.max()])
-
-#--------------------------------------------------
-
-#File
-file = h5py.File(fileName, 'w')
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [Boxsize]*3
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 1, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 1, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = 0
-grp.attrs["Dimension"] = 3
-
-#Runtime parameters
-grp = file.create_group("/RuntimePars")