examples/GEAR/AgoraDisk/agora_disk.yml

@@ -83,14 +83,15 @@ GrackleCooling:
   max_steps: 1000
   convergence_limit: 1e-2
   thermal_time_myr: 0
+  maximal_density_Hpcm3: -1 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
 
 
 GEARStarFormation:
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
-  maximal_temperature:  1e10         # Upper limit to the temperature of a star forming particle
+  maximal_temperature_K:     1e10   # Upper limit to the temperature of a star forming particle
+  density_threshold_Hpcm3:   10     # Density threshold in Hydrogen atoms/cm3
   n_stars_per_particle: 1
   min_mass_frac: 0.5
-  density_threshold:   1.67e-23   # Density threashold in g/cm3
 
 
 GEARPressureFloor:

examples/GEAR/AgoraDisk/getChemistryTable.sh

 #!/bin/bash
 
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5
+wget https://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
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/GEAR/AgoraDisk/getIC.sh

@@ -6,6 +6,5 @@ if [ "$#" -ne 1 ]; then
     exit
 fi
 
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
 mv snap.$1.hdf5 agora_disk.hdf5
-

examples/GEAR/AgoraDisk/getSolution.sh

 #!/bin/bash
 
-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
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0000.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/Gear/snapshot_0050.hdf5

examples/GEAR/ZoomIn/getChemistryTable.sh

 #!/bin/bash
 
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/FeedbackTables/POPIIsw.h5
+wget https://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
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/GEAR/ZoomIn/getIC.sh

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

examples/GEAR/ZoomIn/zoom_in.yml

@@ -91,12 +91,13 @@ GrackleCooling:
   max_steps: 1000
   convergence_limit: 1e-2
   thermal_time_myr: 5
+  maximal_density_Hpcm3: -1 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
 
 
 GEARStarFormation:
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
-  maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
-  density_threshold:   1.67e-25    # Density threashold in g/cm3
+  maximal_temperature_K:     3e4    # Upper limit to the temperature of a star forming particle
+  density_threshold_Hpcm3:   0.1    # Density threshold in Hydrogen atoms/cm3
   n_stars_per_particle: 4
   min_mass_frac: 0.5
 

examples/GravityTests/DiscPatch/HydroStatic/getGlass.sh

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

examples/GravityTests/MWPotential2014_circularorbit/makeIC.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch)
+# 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

examples/GravityTests/MWPotential2014_circularorbit/makePlots.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch)
+# 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

examples/GravityTests/MWPotential2014_circularorbit/params_unit_1.yml

 # Define the system of units to use internally.
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.98848e+43 # 10^10 Solar masses
-  UnitLength_in_cgs:   3.086e+21 # kpc
+  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
@@ -46,7 +46,7 @@ MWPotential2014Potential:
   # Mdisk_Msun:       6.8e10
   # Rdisk_kpc:        3.0
   # Zdisk_kpc:        0.280
-  # amplitude:        1.0
+  # amplitude_Msun_per_kpc3: 1.0e10
   # r_1_kpc:          1.0
   # alpha:            1.8
   # r_c_kpc:          1.9

examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml

 # Define the system of units to use internally.
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.98848e+43 # 10^10 Solar masses
-  UnitLength_in_cgs:   3.086e+24 # Mpc
+  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
@@ -45,7 +45,7 @@ MWPotential2014Potential:
   # Mdisk_Msun:       6.8e10
   # Rdisk_kpc:        3.0
   # Zdisk_kpc:        0.280
-  # amplitude:        1.0
+  # amplitude_Msun_per_kpc3: 1e10
   # r_1_kpc:          1.0
   # alpha:            1.8
   # r_c_kpc:          1.9

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/Planetary/SquareTest_2D/makeICDifferentMasses.py → examples/GravityTests/MWPotential2014_df/makeIC.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
-#               2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+# 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
@@ -16,96 +15,65 @@
 # 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
 
-import h5py
-from numpy import *
+box_size = 1000.0
+N_PARTICLES = 1
+print("Initial conditions written to 'IC.hdf5'")
 
-# Generates a swift IC file for the Square test in a periodic box
+pos = np.zeros((1, 3))
+pos[0, 0] = 0.1
+pos[0, 1] = 0.1
+pos[0, 2] = 100
 
-# Parameters
-L = 2 * 64  # Number of particles on the side
-gamma = 5.0 / 3.0  # Gas adiabatic index
-rho0 = 4.0  # Gas central density
-rho1 = 1.0  # Gas outskirt density
-P0 = 2.5  # Gas central pressure
-P1 = 2.5  # Gas central pressure
-vx = 0.0  # Random velocity for all particles
-vy = 0.0
-fileOutputName = "square.hdf5"
-# ---------------------------------------------------
 
-vol = 1.0
+# pos += np.array(
+#    [box_size / 2, box_size / 2, box_size / 2]
+# )  # shifts the particles to the center of the box
 
-numPart = L * L
+vel = np.zeros((1, 3))
+vel[0, 0] = 80
+vel[0, 1] = 0
+vel[0, 2] = 0
 
-pos_x = arange(0, 1, 1.0 / L)
-xv, yv = meshgrid(pos_x, pos_x)
-pos = zeros((numPart, 3), dtype=float)
-pos[:, 0] = xv.flatten()
-pos[:, 1] = yv.flatten()
-
-# Now we can get 2d masks!
-inside = logical_and.reduce([xv < 0.75, xv > 0.25, yv < 0.75, yv > 0.25])
-
-mass_in = rho0 / numPart
-mass_out = rho1 / numPart
-
-m = ones_like(xv) * mass_out
-m[inside] = mass_in
-m = m.flatten()
-
-h = ones_like(m) / L
-
-u_in = P0 / ((gamma - 1) * rho0)
-u_out = P1 / ((gamma - 1) * rho1)
-u = ones_like(xv) * u_out
-u[inside] = u_in
-u = u.flatten()
-vel = zeros_like(pos)
-ids = arange(numPart)
-mat = zeros(numPart)
 
+ids = np.array([1.0])
+mass = np.array([1.0]) * 1e-10
 
 # File
-fileOutput = h5py.File(fileOutputName, "w")
+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 = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [vol, vol, 0.2]
-grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+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"] = [numPart, 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"] = 2
+grp.attrs["Dimension"] = 3
 
-# Units
-grp = fileOutput.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.0
-grp.attrs["Unit mass in cgs (U_M)"] = 1.0
-grp.attrs["Unit time in cgs (U_t)"] = 1.0
-grp.attrs["Unit current in cgs (U_I)"] = 1.0
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
 
 # Particle group
-grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
-ds[()] = pos
-ds = grp.create_dataset("Velocities", (numPart, 3), "f")
-ds[()] = vel
-ds = grp.create_dataset("Masses", (numPart, 1), "f")
-ds[()] = m.reshape((numPart, 1))
-ds = grp.create_dataset("SmoothingLengths", (numPart, 1), "f")
-ds[()] = h.reshape((numPart, 1))
-ds = grp.create_dataset("InternalEnergies", (numPart, 1), "f")
-ds[()] = u.reshape((numPart, 1))
-ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
-ds[()] = ids.reshape((numPart, 1))
-ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
-ds[()] = mat.reshape((numPart, 1))
-
+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)
 
-fileOutput.close()
+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