Merge branch 'AGORA' into 'master'

Agora

See merge request !572

examples/AgoraDisk/agora_disk.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085e21   # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+Scheduler:
+  max_top_level_cells: 8
+  
+# 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).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            agora_disk # 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
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # 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).
+  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
+  
+# 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.
+  minimal_temperature:   10      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./agora_disk.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  shift:    [674.1175, 674.1175, 674.1175]   # (Optional) A shift to apply to all particles read from the ICs (in internal units).
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda_cgs:                  1.0e-22    # Cooling rate (in cgs units)
+  minimum_temperature:         1.0e2      # Minimal temperature (Kelvin)
+  mean_molecular_weight:       0.59       # Mean molecular weight
+  hydrogen_mass_abundance:     0.75       # Hydrogen mass abundance (dimensionless)
+  cooling_tstep_mult:          1.0        # Dimensionless pre-factor for the time-step condition
+
+# 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

examples/AgoraDisk/change_type.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)
+
+    changeType(f, 2, 1)
+    changeType(f, 3, 1)
+    changeType(f, 4, 1)
+
+    countPart(f)
+
+    f.close()

examples/AgoraDisk/cleanup_swift.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)
+
+name = "PartType0/ElementAbundance"
+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 + "/SmoothingLength"]
+    grp[:] *= 1.823
+
+f.close()

examples/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."
+    exit
+fi
+
+wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/$1

examples/AgoraDisk/getSolution.sh

+#!/bin/bash
+
+OPTIND=1
+
+with_cooling=0
+
+function show_help {
+    echo "Valid options are:"
+    echo "\t -h \t Show this help"
+    echo "\t -C \t Download solution with cooling"
+}
+
+while getopts "h?C" opt; do
+    case "$opt" in
+	h|\?)
+	    show_help
+	    exit
+	    ;;
+	C)
+	    with_cooling=1
+	    ;;
+    esac
+done
+
+# cleanup work space
+rm snapshot_0000.hdf5 snapshot_0500.hdf5
+
+if [ $with_cooling -eq 1 ]; then
+    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
+else
+    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
+fi
+
+

examples/AgoraDisk/run.sh

+#!/bin/bash
+
+# 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
+
+# enable cooling or not
+cooling=0
+
+# make run.sh fail if a subcommand fails
+set -e
+
+# define flags
+flag=
+if [ $cooling -eq 1 ]
+then
+    flag=-C
+fi
+
+# Generate the initial conditions if they are not present.
+if [ ! -e $sim.hdf5 ]
+then
+    echo "Fetching initial glass file for the Sedov blast example..."
+    ./getIC.sh $sim.hdf5
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+    echo "Fetching the Cloudy tables required by Grackle..."
+    ../getCoolingTable.sh
+fi
+
+# copy the initial conditions
+cp $sim.hdf5 agora_disk.hdf5
+# Update the particle types
+python3 change_type.py agora_disk.hdf5
+
+# Run SWIFT
+../swift $flag -s -G -t 4 agora_disk.yml 2>&1 | tee output.log
+
+
+echo "Changing smoothing length to be Gadget compatible"
+python3 cleanup_swift.py agora_disk_0000.hdf5 agora_disk_IC.hdf5
+python3 cleanup_swift.py agora_disk_0050.hdf5 agora_disk_500Myr.hdf5
+
+echo "Fetching GEAR solution..."
+./getSolution.sh $flag
+
+echo "Plotting..."
+python3 plot_solution.py $flag

examples/CoolingBox/run.sh

@@ -17,7 +17,7 @@ fi
 if [ ! -e CloudyData_UVB=HM2012.h5 ]
 then
     echo "Fetching the Cloudy tables required by Grackle..."
-    ./getCoolingTable.sh
+    ../getCoolingTable.sh
 fi
 
 # Run SWIFT

examples/KeplerianRing/README.md

@@ -2,7 +2,8 @@ Keplerian Ring Test
 ===================
 
 This test uses the '2D Keplerian Ring' to assess the accuracy of SPH
-calculations. For more information on the test, please see Cullen & Dehnen 2009
+calculations. It is known to fail with all schemes in SWIFT. For more
+information on the test, please see Cullen & Dehnen 2009
 (http://arxiv.org/abs/1006.1524) Section 4.3.
 
 It is key that for this test in particular that the initial conditions are

examples/KeplerianRing/run.sh

@@ -10,3 +10,9 @@ fi
 
 rm -rf keplerian_ring_*.hdf5
 ../swift -g -s -t 1 -v 1 keplerian_ring.yml 2>&1 | tee output.log
+
+echo
+echo
+echo "This test is expected to fail, please read the README.md file"
+echo
+echo