examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.988409870698051e+33 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.0856775814913673e+24 # Mpc
+  UnitVelocity_in_cgs: 1e5       # 1 km / s in cm/s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            2.0e-3     # The end time of the simulation (in internal units).
+  dt_min:              1e-13    # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e-3    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_2/output  # Common part of the name of output files
+  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          1e-5    # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.98848e+43
+  UnitLength_in_cgs:   3.086e+21
+  UnitVelocity_in_cgs: 1e5
+  UnitCurrent_in_cgs:  1
+  UnitTemp_in_cgs:     1
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.0    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          circular_orbits_MW.hdf5 # The file to read
+  periodic:           1
+
+# NFW_MN_PSC potential parameters
+MWPotential2014Potential:
+  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]    #Centre of the potential with respect to centre of the box
+  timestep_mult:   0.0005     # Dimensionless pre-factor for the time-step condition
+  epsilon:         0.001e-3      # Softening size (internal units)
+  #The following parameters are the default ones.
+  # concentration:    9.823403437774843
+  # M_200_Msun:       147.41031542774076e10
+  # H:                127.78254614201471e-2
+  # Mdisk_Msun:       6.8e10
+  # Rdisk_kpc:        3.0
+  # Zdisk_kpc:        0.280
+  # amplitude_Msun_per_kpc3: 1e10
+  # r_1_kpc:          1.0
+  # alpha:            1.8
+  # r_c_kpc:          1.9
+  # potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262]

examples/GravityTests/MWPotential2014_circularorbit/run.sh

+#!/bin/bash
+
+#Creates a directory for the outputs
+DIR=output_1 #First test of units conversion
+if [ -d "$DIR" ];
+then
+    echo "$DIR directory exists. Its content will be removed."
+    rm $DIR/output_*
+else
+    echo "$DIR directory does not exists. It will be created."
+    mkdir $DIR
+fi
+
+DIR=output_2 #Second test of units conversion
+if [ -d "$DIR" ];
+then
+    echo "$DIR directory exists. Its content will be removed."
+    rm $DIR/output_*
+else
+    echo "$DIR directory does not exists. It will be created."
+    mkdir $DIR
+fi
+
+#Clears the previous figures
+echo "Clearing existing figures."
+if [ -f "circular_orbits_simulation_kpc.png" ];
+then
+    rm circular_orbits_simulation_kpc.png
+fi
+
+if [ -f "circular_orbits_simulation_Mpc.png" ];
+then
+    rm circular_orbits_simulation_Mpc.png
+fi
+
+if [ -f "deviation_simulation_kpc.png" ];
+then
+    rm deviation_simulation_kpc.png
+fi
+if [ -f "deviation_simulation_Mpc.png" ];
+then
+    rm deviation_simulation_Mpc.png
+fi
+
+if [ -f "deviation_from_original_data_simulation_kpc.png" ];
+then
+    rm deviation_from_original_data_simulation_kpc.png
+fi
+
+if [ -f "deviation_from_original_data_simulation_Mpc.png" ];
+then
+    rm deviation_from_original_data_simulation_Mpc.png
+fi
+
+#Clears the IC file
+if [ -f "circular_orbits_MW.hdf5" ];
+then
+    rm circular_orbits_MW.hdf5
+fi
+
+
+#Generates the initial conditions
+echo "Generate initial conditions for circular orbits"
+if command -v python3 &>/dev/null; then
+    python3 makeIC.py
+else
+    python3 makeIC.py
+fi
+
+#Runs the simulation
+# self gravity G, external potential g, hydro s, threads t and high verbosity v
+echo "Starting the simulation with the first type of units (kpc)... Look at output_1.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_1.yml 2>&1 > output_1.log
+echo "Simulation ended."
+
+echo "Starting the simulation with the second type of units (Mpc)... Look at output_2.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_2.yml 2>&1 > output_2.log
+echo "Simulation ended."
+
+#Saves the plots
+echo "Save plots of the circular orbits and of the errors"
+if command -v python3 &>/dev/null; then
+    python3 makePlots.py
+else 
+    python3 makePlots.py
+fi

examples/GravityTests/MWPotential2014_df/README

+# Context
+
+This example tests the dynamical friction implemented in the MWPotential2014.
+It compares the orbit predicted by SWIFT with an orbit computed with pNbody, using:
+
+`orbits_integration_MW --t_forward 10 --position 0.1 0.1 100 --velocity 80 0 0 --dynamical_friction -o orbit.csv`
+
+# How to run this example
+
+In a terminal at the root of the "swiftsim" directory, type:
+
+`./autogen.sh`
+
+Then, configure swift to take into account external potentials. Type:
+
+`./configure --with-ext-potential=MWPotential2014`
+
+Feel free to adapt other configurations parameters depending on your system. Then, build the code by typing:
+
+`make`
+
+Finally, to run this example, open a terminal in the present directory and type:
+
+`./run.sh`
+
+In this last command run Swift run on two different parameter files, params_unit_1.yml and params_unit_2.yml.
+Those two sets of parameters only differ by the choice of the internal units. This allow to test the correct 
+units conversion of all parameters involved.

examples/GravityTests/MWPotential2014_df/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Darwin Roduit (darwin.roduit@alumni.epfl.ch)
+#
+# 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 numpy as np
+import h5py as h5
+
+box_size = 1000.0
+N_PARTICLES = 1
+print("Initial conditions written to 'IC.hdf5'")
+
+pos = np.zeros((1, 3))
+pos[0, 0] = 0.1
+pos[0, 1] = 0.1
+pos[0, 2] = 100
+
+
+# pos += np.array(
+#    [box_size / 2, box_size / 2, box_size / 2]
+# )  # shifts the particles to the center of the box
+
+vel = np.zeros((1, 3))
+vel[0, 0] = 80
+vel[0, 1] = 0
+vel[0, 2] = 0
+
+
+ids = np.array([1.0])
+mass = np.array([1.0]) * 1e-10
+
+# File
+file = h5.File("IC.hdf5", "w")
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 3.086e21
+grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e43
+grp.attrs["Unit time in cgs (U_t)"] = 3.086e16
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = box_size
+grp.attrs["NumPart_Total"] = [0, N_PARTICLES, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, N_PARTICLES, 0, 0, 0, 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, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+# Particle group
+grp1 = file.create_group("/PartType1")
+ds = grp1.create_dataset("Velocities", (N_PARTICLES, 3), "f", data=vel)
+ds = grp1.create_dataset("Masses", (N_PARTICLES,), "f", data=mass)
+ds = grp1.create_dataset("ParticleIDs", (N_PARTICLES,), "L", data=ids)
+ds = grp1.create_dataset("Coordinates", (N_PARTICLES, 3), "d", data=pos)
+
+file.close()

examples/GravityTests/MWPotential2014_df/makePlots.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Darwin Roduit (yves.revaz@.epfl.ch)
+#
+# 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 numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+
+#%%Functions
+
+
+def get_positions_and_time(N_snapshots, N_part, output_dir, boxsize):
+    xx = np.zeros((N_part, N_snapshots))
+    yy = np.zeros((N_part, N_snapshots))
+    zz = np.zeros((N_part, N_snapshots))
+    time = np.zeros(N_snapshots)
+    pot = np.zeros((N_part, N_snapshots))
+    for i in range(0, N_snapshots):
+        Data = h5py.File(output_dir + "/output_%04d.hdf5" % i, "r")
+        header = Data["Header"]
+        time[i] = header.attrs["Time"][0]
+        particles = Data["PartType1"]
+        positions = particles["Coordinates"]
+        xx[:, i] = positions[:, 0] - boxsize / 2.0
+        yy[:, i] = positions[:, 1] - boxsize / 2.0
+        zz[:, i] = positions[:, 2] - boxsize / 2.0
+        pot[:, i] = particles["Potentials"][:]
+    return xx, yy, zz, time, pot
+
+
+def plot_orbits(x, y, z, t, color, save_fig_name_suffix):
+    # Plots the orbits
+    fig, ax = plt.subplots(nrows=1, ncols=4, num=1, figsize=(12, 4.1))
+    fig.suptitle("Orbits", fontsize=15)
+    ax[0].clear()
+    ax[1].clear()
+
+    for i in range(0, N_part):
+        ax[0].plot(x[i, :], y[i, :], color[i])
+
+    ax[0].set_aspect("equal", "box")
+    ax[0].set_xlim([-300, 300])
+    ax[0].set_ylim([-300, 300])
+    ax[0].set_ylabel("y (kpc)")
+    ax[0].set_xlabel("x (kpc)")
+
+    for i in range(0, N_part):
+        ax[1].plot(x[i, :], z[i, :], col[i], label="SWIFT solution")
+
+    ax[1].set_aspect("equal", "box")
+    ax[1].set_xlim([-100, 100])
+    ax[1].set_ylim([-100, 100])
+    ax[1].set_ylabel("z (kpc)")
+    ax[1].set_xlabel("x (kpc)")
+
+    for i in range(0, N_part):
+        ax[2].plot(y[i, :], z[i, :], col[i])
+
+    ax[2].set_aspect("equal", "box")
+    ax[2].set_xlim([-100, 100])
+    ax[2].set_ylim([-100, 100])
+    ax[2].set_ylabel("z (kpc)")
+    ax[2].set_xlabel("y (kpc)")
+    plt.tight_layout()
+
+    for i in range(0, N_part):
+        ax[3].plot(t, np.sqrt(x[i, :] ** 2 + y[i, :] ** 2 + z[i, :] ** 2), col[i])
+
+    ax[3].set_aspect("auto", "box")
+    ax[3].set_ylim([0, 100])
+    ax[3].set_ylabel("r (kpc)")
+    ax[3].set_xlabel("t (kpc)")
+    plt.tight_layout()
+
+    # add the reference orbit
+    data = np.genfromtxt("orbit.csv", delimiter=",", skip_header=1)
+    t = data[:, 0]
+    x = data[:, 1]
+    y = data[:, 2]
+    z = data[:, 3]
+
+    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
+
+    ax[0].plot(x, y, "grey", alpha=0.5, lw=5)
+    ax[1].plot(x, z, "grey", alpha=0.5, lw=5, label="pNbody solution")
+    ax[2].plot(y, z, "grey", alpha=0.5, lw=5)
+    ax[3].plot(t, r, "grey", alpha=0.5, lw=5)
+
+    ax[1].legend()
+
+    plt.savefig("orbits" + save_fig_name_suffix + ".png", bbox_inches="tight")
+    plt.close()
+
+
+#%%Plots the orbits, the deviation from the circular orbit and the deviation from the original precomputed data
+# Notice that in this examples, the ouputs are set in suitable units in the parameters files.
+
+# General parameters
+N_snapshots = 1001
+N_part = 1
+boxsize = 1000.0  # kpc
+col = ["b", "r", "c", "y", "k"]
+
+# First type of units (kpc)
+output_dir = "output_1"
+save_fig_name_suffix = "_simulation_kpc"
+x_1, y_1, z_1, time_1, pot_1 = get_positions_and_time(
+    N_snapshots, N_part, output_dir, boxsize
+)
+plot_orbits(x_1, y_1, z_1, time_1, col, save_fig_name_suffix)
+
+
+# Second type of units (Mpc) (no need for units conversion to kpc, this is already done by swift in the snapshots)
+output_dir = "output_2"
+save_fig_name_suffix = "_simulation_Mpc"
+x_2, y_2, z_2, time_2, pot_2 = get_positions_and_time(
+    N_snapshots, N_part, output_dir, boxsize
+)
+plot_orbits(x_2, y_2, z_2, time_2, col, save_fig_name_suffix)

examples/GravityTests/MWPotential2014_df/params_unit_1.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.988409870698051e+43 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.0856775814913673e+21 # kpc
+  UnitVelocity_in_cgs: 1e5       # 1 km / s in cm/s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            10.0     # 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:              1e0    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output           # Common part of the name of output files
+  subdir:              output_1         # (Optional) Sub-directory in which to write the snapshots. Defaults to "" (i.e. the directory where SWIFT is run).
+  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          1e-2    # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.98848e+43
+  UnitLength_in_cgs:   3.086e+21
+  UnitVelocity_in_cgs: 1e5
+  UnitCurrent_in_cgs:  1
+  UnitTemp_in_cgs:     1
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.0    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          IC.hdf5 # The file to read
+  shift:              [500,500,500]
+  periodic:           0
+
+# NFW_MN_PSC potential parameters
+MWPotential2014Potential:
+  useabspos:       0                   # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]          #Centre of the potential with respect to centre of the box
+  timestep_mult:   0.0005              # Dimensionless pre-factor for the time-step condition
+  epsilon:         0.001e-3            # Softening size (internal units)
+  with_dynamical_friction:1            # Are we running with dynamical friction ? 0 -> no, 1 -> yes
+  df_lnLambda:        5.0              # Coulomb logarithm
+  df_satellite_mass_in_Msun: 1e10      # Satellite mass in solar mass
+  df_core_radius_in_kpc: 10            # Radius below which the dynamical friction vanishes.  
+  df_timestep_mult:0.1                 # Dimensionless pre-factor for the time-step condition for the dynamical friction force
+  df_polyfit_coeffs00: -2.96536595e-31 # Polynomial fit coefficient for the velocity dispersion model (order 16)
+  df_polyfit_coeffs01:  8.88944631e-28 # Polynomial fit coefficient for the velocity dispersion model (order 15)
+  df_polyfit_coeffs02: -1.18280578e-24 # Polynomial fit coefficient for the velocity dispersion model (order 14)
+  df_polyfit_coeffs03:  9.29479457e-22 # Polynomial fit coefficient for the velocity dispersion model (order 13)
+  df_polyfit_coeffs04: -4.82805265e-19 # Polynomial fit coefficient for the velocity dispersion model (order 12)
+  df_polyfit_coeffs05:  1.75460211e-16 # Polynomial fit coefficient for the velocity dispersion model (order 11)
+  df_polyfit_coeffs06: -4.59976540e-14 # Polynomial fit coefficient for the velocity dispersion model (order 10)
+  df_polyfit_coeffs07:  8.83166045e-12 # Polynomial fit coefficient for the velocity dispersion model (order 9)
+  df_polyfit_coeffs08: -1.24747700e-09 # Polynomial fit coefficient for the velocity dispersion model (order 8)
+  df_polyfit_coeffs09:  1.29060404e-07 # Polynomial fit coefficient for the velocity dispersion model (order 7)
+  df_polyfit_coeffs10: -9.65315026e-06 # Polynomial fit coefficient for the velocity dispersion model (order 6)
+  df_polyfit_coeffs11:  5.10187806e-04 # Polynomial fit coefficient for the velocity dispersion model (order 5)
+  df_polyfit_coeffs12: -1.83800281e-02 # Polynomial fit coefficient for the velocity dispersion model (order 4)
+  df_polyfit_coeffs13:  4.26501444e-01 # Polynomial fit coefficient for the velocity dispersion model (order 3)
+  df_polyfit_coeffs14: -5.78038064e+00 # Polynomial fit coefficient for the velocity dispersion model (order 2)
+  df_polyfit_coeffs15:  3.57956721e+01 # Polynomial fit coefficient for the velocity dispersion model (order 1)
+  df_polyfit_coeffs16:  1.85478908e+02 # Polynomial fit coefficient for the velocity dispersion model (order 0)
+  
+

examples/GravityTests/MWPotential2014_df/params_unit_2.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.988409870698051e+33 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.0856775814913673e+24 # Mpc
+  UnitVelocity_in_cgs: 1e5       # 1 km / s in cm/s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            10.0e-3     # The end time of the simulation (in internal units).
+  dt_min:              1e-13    # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e-3    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output           # Common part of the name of output files
+  subdir:              output_2         # (Optional) Sub-directory in which to write the snapshots. Defaults to "" (i.e. the directory where SWIFT is run).  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          1e-5    # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.98848e+43
+  UnitLength_in_cgs:   3.086e+21
+  UnitVelocity_in_cgs: 1e5
+  UnitCurrent_in_cgs:  1
+  UnitTemp_in_cgs:     1
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.0    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          IC.hdf5 # The file to read
+  shift:              [0.5,0.5,0.5]
+  periodic:           0
+  
+# NFW_MN_PSC potential parameters
+MWPotential2014Potential:
+  useabspos:       0                 # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]        #Centre of the potential with respect to centre of the box
+  timestep_mult:   0.0005            # Dimensionless pre-factor for the time-step condition
+  epsilon:         0.001e-3          # Softening size (internal units)
+  with_dynamical_friction:1          # Are we running with dynamical friction ? 0 -> no, 1 -> yes
+  df_lnLambda:        5.0              # Coulomb logarithm
+  df_satellite_mass_in_Msun: 1e10      # Satellite mass in solar mass
+  df_core_radius_in_kpc: 10            # Radius below which the dynamical friction vanishes.  
+  df_timestep_mult:0.1                 # Dimensionless pre-factor for the time-step condition for the dynamical friction force
+  df_polyfit_coeffs00: -2.96536595e-31 # Polynomial fit coefficient for the velocity dispersion model (order 16)
+  df_polyfit_coeffs01:  8.88944631e-28 # Polynomial fit coefficient for the velocity dispersion model (order 15)
+  df_polyfit_coeffs02: -1.18280578e-24 # Polynomial fit coefficient for the velocity dispersion model (order 14)
+  df_polyfit_coeffs03:  9.29479457e-22 # Polynomial fit coefficient for the velocity dispersion model (order 13)
+  df_polyfit_coeffs04: -4.82805265e-19 # Polynomial fit coefficient for the velocity dispersion model (order 12)
+  df_polyfit_coeffs05:  1.75460211e-16 # Polynomial fit coefficient for the velocity dispersion model (order 11)
+  df_polyfit_coeffs06: -4.59976540e-14 # Polynomial fit coefficient for the velocity dispersion model (order 10)
+  df_polyfit_coeffs07:  8.83166045e-12 # Polynomial fit coefficient for the velocity dispersion model (order 9)
+  df_polyfit_coeffs08: -1.24747700e-09 # Polynomial fit coefficient for the velocity dispersion model (order 8)
+  df_polyfit_coeffs09:  1.29060404e-07 # Polynomial fit coefficient for the velocity dispersion model (order 7)
+  df_polyfit_coeffs10: -9.65315026e-06 # Polynomial fit coefficient for the velocity dispersion model (order 6)
+  df_polyfit_coeffs11:  5.10187806e-04 # Polynomial fit coefficient for the velocity dispersion model (order 5)
+  df_polyfit_coeffs12: -1.83800281e-02 # Polynomial fit coefficient for the velocity dispersion model (order 4)
+  df_polyfit_coeffs13:  4.26501444e-01 # Polynomial fit coefficient for the velocity dispersion model (order 3)
+  df_polyfit_coeffs14: -5.78038064e+00 # Polynomial fit coefficient for the velocity dispersion model (order 2)
+  df_polyfit_coeffs15:  3.57956721e+01 # Polynomial fit coefficient for the velocity dispersion model (order 1)
+  df_polyfit_coeffs16:  1.85478908e+02 # Polynomial fit coefficient for the velocity dispersion model (order 0)

examples/GravityTests/MWPotential2014_df/run.sh

+#!/bin/bash
+
+#Clears the previous figures
+echo "Clearing existing figures."
+if [ -f "orbits_simulation_kpc.png" ]; then
+    rm orbits_simulation_kpc.png
+fi
+
+if [ -f "orbits_simulation_Mpc.png" ]; then
+    rm orbits_simulation_Mpc.png
+fi
+
+if [ ! -e orbit.csv ]; then
+    echo "Fetching solution..."
+    wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/MWPotential_2014/orbit.csv
+fi
+
+#Clears the IC file
+if [ -f "circular_orbits_MW.hdf5" ]; then
+    rm circular_orbits_MW.hdf5
+fi
+
+#Generates the initial conditions
+echo "Generate initial conditions for circular orbits"
+if command -v python3 &>/dev/null; then
+    python3 makeIC.py
+else
+    python3 makeIC.py
+fi
+
+#Runs the simulation
+# self gravity G, external potential g, hydro s, threads t and high verbosity v
+echo "Starting the simulation with the first type of units (kpc)... Look at output_1.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_1.yml 2>&1 >output_1.log
+echo "Simulation ended."
+
+echo "Starting the simulation with the second type of units (Mpc)... Look at output_2.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_2.yml 2>&1 >output_2.log
+echo "Simulation ended."
+
+#Saves the plots
+echo "Save plots of orbits."
+if command -v python3 &>/dev/null; then
+    python3 makePlots.py
+else
+    python3 makePlots.py
+fi

examples/GravityTests/NFW_Halo/run.sh

@@ -6,15 +6,15 @@ then
     if command -v python3 &>/dev/null; then
         python3 makeIC.py
     else 
-        python makeIC.py
+        python3 makeIC.py
     fi
 fi
 
 # self gravity G, external potential g, hydro s, threads t and high verbosity v
-../../swift --external-gravity --threads=6 test.yml 2>&1 | tee output.log
+../../../swift --external-gravity --threads=6 test.yml 2>&1 | tee output.log
 
 if command -v python3 &>/dev/null; then
     python3 makePlots.py
 else 
-    python makePlots.py
+    python3 makePlots.py
 fi

examples/GravityTests/NFW_Halo/test.yml

@@ -37,5 +37,6 @@ NFWPotential:
   position:           [0.0,0.0,0.0]      # Location of centre of potential with respect to centre of the box (internal units)
   concentration:      8.
   M_200:              2.0e+12  # Virial mass (internal units)
-  critical_density:   140      # Critical density (internal units)
+  epsilon:            0.8
+  h:                  0.7
   timestep_mult:      0.01     # Dimensionless pre-factor for the time-step condition

examples/GravityTests/Plummer_Selfgravitating/README

+This example generates a self consistent realization of an anisotropic Plummer model, by sampling
+the distribution function derived by (Dejonghe 1987).
+
+The parameters that control the ICs are essentially:
+*  a: 			The Plummer softening length
+*  M:			The total mass of the system
+*  q:			A parameter that controls the radial anisotropy of the system,
+			ranging from -inf to 2, where -inf corresponds to the "Einstein Sphere",
+			a model with only circular orbits and 2 to a model with only radial orbits.
+			q is related to the traditional beta parameter (which measures anisotropy) via
+			beta(r) = q/2 * r^2/(1+r^2)
+*  N:			The Number of particles in the system
+
+Some additional parameters are:
+*  bound:		A given cut radius above which particles are not included in the IC file
+			(this may reduce the effective number of particles to N_eff < N). N_eff is
+			printed out when the script is executed.
+*  use_parallel:	Use python's multiprocessing library for sampling
+*  nb_threads:		Number of threads to use for the above
+*  pickle_ics:		Also output .pkl files of the pos, vel and mass arrays
+
+All other variables should be self-explanatory.
+
+The radii of the particles are renerated through inverse transform sampling, exploiting
+the analiticity of the inverse of M(<r) for the Plummer model. (Masses are drawn uniformly
+between 0 and M, and their corresponding radii are computed via r(M))
+
+The radial and tangential components of the velocities are generated trough rejection
+sampling on the DF of the model. For each radius sampled above, a max on the DF is computed
+numerically to achieve this.
+
+The masses of all the particles are simply given by M/N.
+
+The script will also print out a guiding estimate of a good softening length to use, given by the radius
+of a sphere which contains on average one particle at the number density at the softening length a.
+
+NOTE: the positions of the particles produced by this script are centered at the origin, and must
+be shifted by some distance in x,y,z to be accepted by SWIFT. This will depend on the box size and
+the bound parameter described above.
+
+## USAGE ##
+Set the parameters described above as desired in plummerIC.py
+Run plummerIC.py to produce .hdf5 initial conditions
+Set suitable parameters in params.yml
+Run SWIFT with --self-gravity on
+
+To display results, another script (plotdensity.py) is provided. The basic usage is
+python3 plotdensity.py [path to snapshot files to plot (max 20)] -a [softening value] -M [mass value] -shift [shift applied in params.yml]
+
+A run with default values, 10k particles and subsequent density plotting can be automated with the provided run.sh script.

examples/GravityTests/Plummer_Selfgravitating/params.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.988e+43 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.086e+21 # kpc
+  UnitVelocity_in_cgs: 1e5       # km / s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            2.0     # The end time of the simulation (in internal units).
+  dt_min:              1e-11   # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e-3    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                           0.002      # Constant dimensionless multiplier for time integration.
+  MAC:                           geometric  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+  epsilon_fmm:                   0.001      # Tolerance parameter for the adaptive multipole acceptance criterion.
+  theta_cr:                      0.7        # Opening angle for the purely gemoetric criterion.
+  max_physical_DM_softening:     0.005      # Physical softening length (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snapshot  # Common part of the name of output files
+  time_first:          0.        # Time of the first output (in internal units)
+  delta_time:          2e-1      # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          plummer.hdf5 # The file to read
+  shift:              [2.,2.,2.]
+  periodic:           0

examples/GravityTests/Plummer_Selfgravitating/plotdensity.py

+import numpy as np
+from matplotlib import pyplot as plt
+import argparse
+from scipy.optimize import curve_fit
+import h5py
+
+# Parse user input
+parser = argparse.ArgumentParser(
+    description="Plot multiple density profiles against theoretical prediction"
+)
+parser.add_argument("files", nargs="+", help="snapshot files to be imaged")
+parser.add_argument("--notex", action="store_true", help="Flag to not use LaTeX markup")
+parser.add_argument(
+    "-a", type=float, default=0.1, help="Softening length of theoretical model"
+)
+parser.add_argument(
+    "-M", type=float, default=1.0e-5, help="Total Mass of theoretical model"
+)
+parser.add_argument("-Rmin", type=float, default=0.01, help="Min Radius")
+parser.add_argument("-Rmax", type=float, default=2.1, help="Max Radius")
+parser.add_argument(
+    "-nbsamples", type=int, default=64, help="Number of radii to sample (bins)"
+)
+parser.add_argument(
+    "-shift", type=float, default=2.0, help="Shift applied to particles in params.yml"
+)
+
+args = parser.parse_args()
+fnames = args.files
+
+# Limit number of snapshots to plot
+if len(fnames) > 20:
+    raise ValueError(
+        "Too many ({:d}) files provided (cannot plot more than 20).".format(len(fnames))
+    )
+
+# Set parameters
+tex = not args.notex
+if tex:
+    plt.rcParams.update({"text.usetex": tex, "font.size": 14, "font.family": "serif"})
+else:
+    plt.rcParams.update({"font.size": 12})
+figsize = 7
+
+# Model Parameters (Mestel surface density)
+G = 4.299581e04
+rsp = np.logspace(np.log10(args.Rmin), np.log10(args.Rmax), args.nbsamples)
+
+
+def plummer_analytical(r):
+    return (
+        3.0
+        * args.M
+        / (4.0 * np.pi * args.a ** 3)
+        * (1.0 + r ** 2 / args.a ** 2) ** (-2.5)
+    )
+
+
+# Plot densities
+fig, ax = plt.subplots(figsize=(1.2 * figsize, figsize))
+for fname in fnames:
+    print(fname)
+    f = h5py.File(fname, "r")
+    pos = np.array(f["DMParticles"]["Coordinates"]) - args.shift
+    time = (
+        f["Header"].attrs["Time"][0]
+        * f["Units"].attrs["Unit time in cgs (U_t)"][0]
+        / 31557600.0e9
+    )
+    mass = np.array(f["DMParticles"]["Masses"])
+    x = pos[:, 0]
+    y = pos[:, 1]
+    r = np.sqrt(np.sum(pos ** 2, 1))
+
+    # Methods to compute density profile
+    def mass_ins(R):
+        return ((r < R) * mass).sum()
+
+    mass_ins_vect = np.vectorize(mass_ins)
+
+    def density(R):
+        return np.diff(mass_ins_vect(R)) / np.diff(R) / (4.0 * np.pi * R[1:] ** 2)
+
+    dens = density(rsp)
+    rs = rsp[1:]
+
+    # remove empty bins
+    c = dens > 0
+    dens = np.compress(c, dens)
+    rs = np.compress(c, rs)
+
+    # Plot
+    # ax.loglog(rsp[1:], density(rsp), "o", ms=1.7, label=r"$t=$ {:.3f} Gyr".format(time))
+    ax.plot(rs, dens, label=r"$t=$ {:.3f} Gyr".format(time))
+
+ax.plot(rsp, plummer_analytical(rsp), c="black", label="Analytical Prediction")
+ax.set_xlabel("r [kpc]")
+ax.legend()
+ax.loglog()
+ax.set_title(
+    r"Plummer Density Profile: $a = {:.1e}$ kpc, $M = {:.1e}$ M$_{{\odot}}$".format(
+        args.a, args.M * 1e10
+    )
+)
+plt.tight_layout()
+ax.set_ylabel(r"$\rho(r)$ [$M_{\odot}$ kpc$^{-3}$]")
+plt.savefig("density_profile.png")

examples/GravityTests/Plummer_Selfgravitating/plummerIC.py

+################################################################################
+# Copyright (c) 2022 Patrick Hirling (patrick.hirling@epfl.ch)
+#
+# 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 numpy as np
+import scipy.special as sci
+from scipy.optimize import minimize
+import time
+import argparse
+from swiftsimio import Writer
+import unyt
+
+parser = argparse.ArgumentParser(
+    description="Generate discrete realization of Anisotropic Plummer Model"
+)
+
+# Parse Main Arguments
+parser.add_argument("-a", type=float, default=0.05, help="Softening Length (in kpc)")
+parser.add_argument(
+    "-M", type=float, default=1.0e-5, help="Total Mass of Model (in 10^10 solar masses)"
+)
+parser.add_argument(
+    "-q", type=float, default=1.0, help="Anisotropy parameter (between -inf and +2)"
+)
+parser.add_argument(
+    "-N", type=int, default=1e4, help="Number of particles in the simulation"
+)
+parser.add_argument(
+    "-bound", type=float, default=2.0, help="Upper bound of radius sampling"
+)
+parser.add_argument("-boxsize", type=float, default=4.0, help="Box Size of simulation")
+parser.add_argument(
+    "--noparallel",
+    action="store_false",
+    help="Do not use multiprocessing library (default: false)",
+)
+parser.add_argument(
+    "-nbthreads",
+    type=int,
+    default=6,
+    help="Number of threads to use for multiprocessing (default: 6)",
+)
+
+args = parser.parse_args()
+if args.q > 2.0:
+    raise ValueError(
+        "Anisotropy parameter must be smaller than 2.0 (q = {:.2f} given)".format(
+            args.q
+        )
+    )
+
+#### Parameters
+# Plummer Model
+G = 4.299581e04  # Gravitational constant [kpc / 10^10 M_s * (kms)^2]
+a = args.a  # Plummer softening length [kpc]
+M = args.M  # Total Mass [10^10 M_s]
+q = args.q  # Anisotropy Parameter (-inf,2]
+
+# IC File
+N = int(args.N)  # Number of Particles
+fname = "plummer.hdf5"  # Name of the ic file (dont forget .hdf5)
+pickle_ics = 0  # Optional: Pickle ics (pos,vel,mass) for future use
+periodic_bdry = 0
+
+# Parallelism for velocity sampling
+use_parallel = args.noparallel
+nb_threads = int(args.nbthreads)  # Set to None to use os.cpucount()
+
+### Print parameters
+print(
+    "\n############################################################################################# \n"
+)
+print("Generating Plummer ICs with the following parameters::\n")
+print(
+    "Softening length:                                  a = " + "{:.3e} kpc".format(a)
+)
+print(
+    "Total Mass:                                        M = "
+    + "{:.3e} Solar Masses".format(M * 1e10)
+)
+print("Anisotropy parameter:                              q = " + str(q))
+print("Output file:                                       " + fname)
+print("Number of particles:                               N = " + "{:.1e}".format(N))
+print(
+    "\n############################################################################################# \n"
+)
+
+# Estimate good softening length (density at origin, sphere s.t. contains on avg 1 particle)
+epsilon0 = a * N ** (-1.0 / 3.0)
+print(
+    "Recommended Softening length (times conventional factor): "
+    + "{:.4e}".format(epsilon0)
+    + " kpc"
+)
+
+### Generate Positions
+# Inverse transform sampling from analytical inverse of M(<r)
+def r_of_m(m, a, M):
+    return a * ((M / m) ** (2.0 / 3.0) - 1) ** (-1.0 / 2.0)
+
+
+m_rand = M * np.random.uniform(0.0, 1.0, N)
+r_rand = r_of_m(m_rand, a, M)
+phi_rand = np.random.uniform(0.0, 2 * np.pi, N)
+theta_rand = np.arccos(np.random.uniform(-1.0, 1.0, N))
+
+x = r_rand * np.sin(theta_rand) * np.cos(phi_rand)
+y = r_rand * np.sin(theta_rand) * np.sin(phi_rand)
+z = r_rand * np.cos(theta_rand)
+
+X = np.array([x, y, z]).transpose()
+
+### Generate Velocities
+# Sample the DF of the Plummer model at the radii generated above (give vr,vt)
+# Dejonghe (1987) Anisotropic Distribution Functionl
+normalization = 3.0 * sci.gamma(6.0 - q) / (2.0 * (2.0 * np.pi) ** (2.5))
+units = G ** (q - 5.0) * M ** (q - 4.0) * a ** (2.0 - q)
+
+
+def H(E, L):
+    x = L ** 2 / (2.0 * E * a ** 2)
+    if x <= 1.0:
+        return 1.0 / (sci.gamma(4.5 - q)) * sci.hyp2f1(0.5 * q, q - 3.5, 1.0, x)
+    else:
+        return (
+            1.0
+            / (sci.gamma(1.0 - 0.5 * q) * sci.gamma(4.5 - 0.5 * q) * (x ** (0.5 * q)))
+            * sci.hyp2f1(0.5 * q, 0.5 * q, 4.5 - 0.5 * q, 1.0 / x)
+        )
+
+
+def Fq(E, L):
+    if E < 0.0:
+        return 0.0
+    else:
+        return normalization * units * E ** (3.5 - q) * H(E, L)
+
+
+# Helper Functions
+# Total energy: E = phi - 1/2v^2. relative potential: psi = -phi. Relative energy: -E = psi - 1/2v^2
+def relative_potential(r):  # = - phi
+    return G * M / np.sqrt(r ** 2 + a ** 2)
+
+
+def relative_energy(r, vr, vt):
+    return relative_potential(r) - 0.5 * (vr ** 2 + vt ** 2)
+
+
+# N.B: Angular momentum: L = r * vt
+
+# Convenience Function for scipy.minimize negative of Fq*vt, indep. vars passed as array
+def Fq_tomin(v, r):
+    return -Fq(relative_energy(r, v[0], v[1]), r * v[1]) * v[1]
+
+
+# Find max of DF at given radius
+def fmax(r, vmax):
+    # Constraint function (defined locally, dep on r)
+    def vel_constr2(v):
+        return vmax ** 2 - v[0] ** 2 - v[1] ** 2
+
+    # Initial Guess
+    v0 = [0.1 * vmax, 0.2 * vmax]
+
+    # Constraint Object
+    cons = {"type": "ineq", "fun": vel_constr2}
+
+    # Args
+    args = (r,)
+
+    # Minimize through scipy.optimize.minimize
+    # fm = minimize(Fq_tomin,v0,constraints=cons,method = 'COBYLA',args=args)
+    fm = minimize(
+        Fq_tomin,
+        v0,
+        constraints=cons,
+        method="SLSQP",
+        args=args,
+        bounds=[(0, vmax), (0, vmax)],
+    )
+
+    # Min of negative df == max of df
+    return -fm.fun
+
+
+# Sample vr,vt from DF at given Radius through rejection sampling
+def sample_vel(r):
+    # Compute max velocity (equiv. condition for E>=0)
+    vmax = np.sqrt(2.0 * relative_potential(r))
+    # Compute max of DF at this radius
+    fm = 1.1 * fmax(r, vmax)  # 1.1 factor to be sure to include max
+    while True:
+        # Select candidates for vr,vt based on max full velocity
+        while True:
+            vr = np.random.uniform(0.0, vmax)
+            vt = np.random.uniform(0.0, vmax)
+            if vr ** 2 + vt ** 2 <= vmax ** 2:
+                break
+        # Rejection Sampling on Fq
+        f = np.random.uniform(0.0, fm)
+        if Fq(relative_energy(r, vr, vt), r * vt) * vt >= f:
+            return vr, vt
+
+
+print("Sampling velocities...")
+ti = time.time()
+if use_parallel:
+    from multiprocessing import Pool
+
+    with Pool(nb_threads) as p:
+        vels = np.array(p.map(sample_vel, r_rand)).transpose()
+else:
+    from tqdm import tqdm
+
+    vels = np.empty((2, N))
+    for j, r in enumerate(tqdm(r_rand)):
+        vels[:, j] = sample_vel(r)
+tf = time.time()
+print("Sampling took " + "{:.2f}".format(tf - ti) + " seconds.")
+
+# Convert to Cartesian
+# First: project vt on e_theta, e_phi with random orientation
+alph = np.random.uniform(0, 2 * np.pi, N)
+sgn = np.random.choice([-1, 1], size=N)
+vphi = vels[1] * np.cos(alph)
+vtheta = vels[1] * np.sin(alph)
+# project vr on e_r (random sign)
+vr = sgn * vels[0]
+
+# Convert Spherical to cartesian coordinates
+v_x = (
+    np.sin(theta_rand) * np.cos(phi_rand) * vr
+    + np.cos(theta_rand) * np.cos(phi_rand) * vtheta
+    - np.sin(phi_rand) * vphi
+)
+v_y = (
+    np.sin(theta_rand) * np.sin(phi_rand) * vr
+    + np.cos(theta_rand) * np.sin(phi_rand) * vtheta
+    + np.cos(phi_rand) * vphi
+)
+v_z = np.cos(theta_rand) * vr - np.sin(theta_rand) * vtheta
+
+# Create velocity array
+V = np.array([v_x, v_y, v_z]).transpose()
+
+### Generate masses
+
+m = M / N * np.ones(N)
+
+### Exclude extreme outliers
+idx = r_rand < args.bound
+X = X[idx]
+V = V[idx]
+m = m[idx]
+new_N = len(m)
+
+### Write to hdf5
+galactic_units = unyt.UnitSystem(
+    "galactic",
+    unyt.kpc,
+    unyt.unyt_quantity(1e10, units=unyt.Solar_Mass),
+    unyt.unyt_quantity(1.0, units=unyt.s * unyt.Mpc / unyt.km).to(unyt.Gyr),
+)
+wrt = Writer(galactic_units, args.boxsize * unyt.kpc)
+wrt.dark_matter.coordinates = X * unyt.kpc
+wrt.dark_matter.velocities = V * (unyt.km / unyt.s)
+wrt.dark_matter.masses = m * 1e10 * unyt.msun
+wrt.dark_matter.particle_ids = np.arange(new_N)
+wrt.write(fname)
+print("Writing IC file...")
+
+### Optional: Pickle ic arrays for future use
+if pickle_ics:
+    import pickle as pkl
+
+    with open("X.pkl", "wb") as f:
+        pkl.dump(X[:], f)
+
+    with open("V.pkl", "wb") as f:
+        pkl.dump(V[:], f)
+
+    with open("M.pkl", "wb") as f:
+        pkl.dump(m[:], f)

examples/GravityTests/Plummer_Selfgravitating/run.sh

+#!/bin/bash
+
+if [ ! -e plummer.hdf5 ]
+then
+    echo "Generating initial conditions for Plummer example..."
+    python3 plummerIC.py -a 0.1 -N 65536
+fi
+
+# create output directory
+if [ ! -e snap ]
+then
+  mkdir snap
+fi
+
+../../../swift --self-gravity --threads=8 params.yml 2>&1 | tee output.log
+
+echo "Plotting results..."
+# If params.yml is left at default values, should produce 10 snapshots
+python3 plotdensity.py snapshot_*.hdf5

examples/HydroTests/BlobTest_3D/makeIC.py

@@ -33,9 +33,11 @@ def generate_cube(num_on_side, side_length=1.0):
 
 
 def generate_bcc_lattice(num_on_side, side_length=1.0):
-    cube = generate_cube(num_on_side // 2, side_length)
+    num_per_cube = num_on_side // 2
 
-    mips = side_length / num_on_side
+    cube = generate_cube(num_per_cube, side_length)
+
+    mips = side_length / num_per_cube
 
     positions = np.concatenate([cube, cube + mips * 0.5])
 
@@ -166,4 +168,3 @@ def write_out_ics(
 
 if __name__ == "__main__":
     write_out_ics(filename="blob.hdf5", num_on_side=64)
-

examples/HydroTests/BlobTest_3D/run.sh

@@ -8,6 +8,6 @@ then
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=2 blob.yml
+../../../swift --hydro --threads=2 blob.yml
 
 python3 makeMovie.py

examples/HydroTests/BlobTest_Braspenning2022/README.md

+Blob Test
+=========
+
+The test in this folder is identical to the one presented in Braspenning+ 2022
+(https://ui.adsabs.harvard.edu/abs/2022arXiv220313915B/abstract)
+and can be used to reproduce those results
+
+You can use the provided ICs, where the first value indicates the particle number
+and the second the density contrast. Switching ICs requires indicating the new ICs
+in the parameter file 'blob.yml'
+You can run this test with any of the hydro schemes supported by SWIFT.
+
+Figures similar to those in Braspenning+ 2022 can be made using 'plot_cloud_evolution.py',
+this produces:
++ Evolution of the mass of dense gas
++ Evolution of the mass of intermediate-temperature gas
+The second and third row show the same evolutions by downsampled to a lower resolution grid.
+This is identical to Fig. 7 in Braspenning+ 2022
+
+The simulation can be run using this command:
+./swift --pin --hydro --limiter --threads=28 blob.yml

examples/HydroTests/BlobTest_Braspenning2022/blob.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e24 # 10^-9 Msun in grams
+  UnitLength_in_cgs:   3.08567758e18 # pc in centimeters
+  UnitVelocity_in_cgs: 1e5   # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   .1   # The end time of the simulation (in internal units).
+  dt_min:     1e-9  # 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:            blob  # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          0.001 # Time difference between consecutive outputs (in internal units)
+  compression:         1
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+
+Scheduler:
+  tasks_per_cell:      100
+
+# 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
+  H_mass_fraction:       0.76
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./blob_16_100.hdf5       # The file to read
+  periodic:   1
+  cleanup_smoothing_lengths: 1

examples/HydroTests/BlobTest_Braspenning2022/getICs.sh

+#!/bin/bash
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_16_100.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_16_10.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_32_100.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_32_10.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_64_100.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_64_10.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_128_100.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/BlobTest_Braspenning22/blob_128_10.hdf5