diff --git a/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml b/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml new file mode 100644 index 0000000000000000000000000000000000000000..5e31dd03d6bc21c8d8429ab375c83e0a2df38546 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml @@ -0,0 +1,58 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams + UnitLength_in_cgs: 3.08567758e21 # kpc 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.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: 0.1 # The maximal time-step size of the simulation (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-2 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-3 # Time between statistics output + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.05 # Constant dimensionless multiplier for time integration. + MAC: geometric + theta_cr: 0.7 + comoving_DM_softening: 0.01 # Comoving softening length (in internal units). + max_physical_DM_softening: 0.01 # Physical softening length (in internal units). + comoving_baryon_softening: 0.01 # Comoving softening length (in internal units). + max_physical_baryon_softening: 0.01 # Physical softening length (in internal units). + +# 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. # Kelvin + +# Parameters related to the initial conditions +InitialConditions: + file_name: ./nfw.hdf5 # The file to read + periodic: 0 # Non-periodic BCs + shift: [0,0,0] # Centre the box + +NFWPotential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) + M_200: 9.5 # M200 of the galaxy disk + h: 0.72 # reduced Hubble constant (value does not specify the used units!) + concentration: 17.0 # concentration of the Halo + diskfraction: 0.15 # Disk mass fraction (here this is the gas) + bulgefraction: 0.0 # Bulge mass fraction + timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.001 # Softening size (internal units) diff --git a/examples/HydroTests/NFW_Hydrostatic/README b/examples/HydroTests/NFW_Hydrostatic/README new file mode 100644 index 0000000000000000000000000000000000000000..17655f5c5c83a3bd0894caa79c6d11e36451d67e --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/README @@ -0,0 +1,13 @@ + +Evolve gas at the static equilibrium in an NFW potential for several dynamical times. +This test allows to compare the impact of different hydro solvers on the initial +hydrostatic equilibrium. + +To run SWIFT configured with the sphenix hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=sphenix + +To run SWIFT configured with the gadget2 hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=gadget2 --disable-hand-vec + +To run SWIFT configured with the gizmo (MFV) hydro solver: + ./configure --with-ext-potential=nfw --with-hydro=gizmo-mfv --with-riemann-solver=exact diff --git a/examples/HydroTests/NFW_Hydrostatic/getICs.sh b/examples/HydroTests/NFW_Hydrostatic/getICs.sh new file mode 100755 index 0000000000000000000000000000000000000000..4f132dc2c82b01e80a0e051137d4bf52c3b91fe5 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/getICs.sh @@ -0,0 +1,3 @@ +#!/bin/bash +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/NFW_Hydrostatic/nfw.hdf5 + diff --git a/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py b/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py new file mode 100755 index 0000000000000000000000000000000000000000..beb677f5127f98626b3dbd2e8f544bbb78c8de56 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py @@ -0,0 +1,126 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2023 Yves Revaz (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 matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from glob import glob +import h5py + + +plt.style.use("../../../tools/stylesheets/mnras.mplstyle") + +MyrInSec = 31557600000000.0 +gcm3InAcc = 1 / 5.978637406556783e23 + + +def ComputeDensity(snap): + + # Read the initial state of the gas + f = h5py.File(snap, "r") + + # Read the units parameters from the snapshot + units = f["InternalCodeUnits"] + unit_mass = units.attrs["Unit mass in cgs (U_M)"] + unit_length = units.attrs["Unit length in cgs (U_L)"] + unit_time = units.attrs["Unit time in cgs (U_t)"] + unit_density = unit_mass / unit_length ** 3 + + # Header + header = f["Header"] + BoxSize = header.attrs["BoxSize"] + Time = header.attrs["Time"] + + # Read data + pos = f["/PartType0/Coordinates"][:] + rho = f["/PartType0/Densities"][:] + mass = f["/PartType0/Masses"][:] + ids = f["/PartType0/ParticleIDs"][:] + + # Center the model and compute particle radius + pos = pos - BoxSize / 2 + x = pos[:, 0] + y = pos[:, 1] + z = pos[:, 2] + r = np.sqrt(x * x + y * y + z * z) + n = len(r) + + # sort particles according to their distance + idx = np.argsort(r) + r = r[idx] + mass = mass[idx] + ids = ids[idx] + + nparts_per_bin = 100 + + # loop over bins containing nparts_per_bin particles + idx = np.arange(n) + + # bins radius and mass + rs_beg = np.array([]) + rs_end = np.array([]) + ms = np.array([]) + + i = 0 + + while i + nparts_per_bin < n: + + rs_beg = np.concatenate((rs_beg, [r[i]])) + rs_end = np.concatenate((rs_end, [r[i + nparts_per_bin]])) + m_this_bin = np.sum(mass[i : i + nparts_per_bin]) + ms = np.concatenate((ms, [np.sum(m_this_bin)])) + + # shift + i = i + nparts_per_bin + + # compute density + vol = 4 / 3 * np.pi * (rs_end ** 3 - rs_beg ** 3) + rho = ms / vol + + # compute radius, we use the mean + rs = 0.5 * (rs_beg + rs_end) + + # convert rho to acc + rho = rho * unit_density / gcm3InAcc + + # convert time to Myr + Time = Time * unit_time / MyrInSec + + return rs, rho, Time + + +# Do the plot + +plt.figure() + +rs, rho, Time = ComputeDensity("snapshot_0000.hdf5") +plt.plot(rs, rho, c="b", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1) + +rs, rho, Time = ComputeDensity("snapshot_0050.hdf5") +plt.plot(rs, rho, c="r", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1) + +plt.loglog() + +plt.legend() +plt.xlabel("${\\rm{Radius~[kpc]}}$", labelpad=0) +plt.ylabel("${\\rm{Density~[atom/cm^3]}}$", labelpad=0) + + +plt.savefig("GasDensityProfile.png", dpi=200) diff --git a/examples/HydroTests/NFW_Hydrostatic/run.sh b/examples/HydroTests/NFW_Hydrostatic/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..713957f30caaa90782902fcccab4d44b2ba366e4 --- /dev/null +++ b/examples/HydroTests/NFW_Hydrostatic/run.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +# Generate the initial conditions if they are not present. +if [ ! -e nfw.hdf5 ] +then + echo "Fetching initial conditions file for the example..." + ./getICs.sh +fi + + +# Run SWIFT +../../../swift --hydro --external-gravity --self-gravity --threads=14 NFW_Hydrostatic.yml + +# Compute gas density profile +python3 plotGasDensityProfile.py