Skip to content
Snippets Groups Projects
Commit 7e0656ff authored by Yves Revaz's avatar Yves Revaz Committed by Matthieu Schaller
Browse files

NFW hydrostatic

parent 017c0f72
No related branches found
No related tags found
1 merge request!1807NFW hydrostatic
# 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)
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
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/NFW_Hydrostatic/nfw.hdf5
###############################################################################
# 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)
#!/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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment