Commit 2adc1473 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Expanded the HydrostaticHalo a little bit and added it to the 'make dist' target.

parent 893bbb76
Hydrostatic halo in equilibrium in an isothermal potential. Running
for 10 dynamical times.
To make the initial conditions we distribute gas particles randomly in
a cube with a side length twice that of the virial radius. The density
profile of the gas is proportional to r^(-2) where r is the distance
from the centre of the cube.
The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
The parameter v_rot (in makeIC.py and hydrostatic.yml) sets the circular
velocity of the halo, and by extension, the viral radius, viral mass,
and the internal energy of the gas such that hydrostatic equilibrium
is achieved.
......@@ -12,10 +14,12 @@ is achieved.
To run this example, make such that the code is compiled with either
the isothermal potential or softened isothermal potential set in
src/const.h. In the latter case, a (small) value of epsilon needs to
be set in cooling.yml. 0.1 kpc should work well.
be set in hydrostatic.yml. ~1 kpc should work well.
The plotting scripts produce a plot of the density, internal energy
and radial velocity profile for each
snapshot. test_energy_conservation.py shows the evolution of energy
and radial velocity profile for each snapshot and divides the profile
by the expected profile.
The script test_energy_conservation.py shows the evolution of energy
with time. These can be used to check if the example has run properly.
......@@ -9,7 +9,7 @@ InternalUnitSystem:
# 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).
time_end: 30. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
......@@ -20,21 +20,18 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
basename: Hydrostatic # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 1. # The tolerance for the targetted number of neighbours.
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: Hydrostatic.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
SoftenedIsothermalPotential:
......@@ -42,6 +39,6 @@ SoftenedIsothermalPotential:
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
epsilon: 0.1
epsilon: 1.0
timestep_mult: 0.03 # controls time step
......@@ -116,9 +116,9 @@ print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
print np.mean(coords[:,0])
print np.mean(coords[:,1])
print np.mean(coords[:,2])
#print np.mean(coords[:,0])
#print np.mean(coords[:,1])
#print np.mean(coords[:,2])
#now find the particles which are within the box
......
......@@ -2,14 +2,23 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 10000
python makeIC.py 100000
../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
# Run for 10 dynamical times
../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
python density_profile.py 2. 200 100
echo "Plotting density profiles"
mkdir plots
mkdir plots/density_profile
python density_profile.py 2. 200 300
python internal_energy_profile.py 2. 200 100
echo "Plotting internal energy profiles"
mkdir plots/internal_energy
python internal_energy_profile.py 2. 200 300
python velocity_profile.py 2. 200 100
echo "Plotting radial velocity profiles"
mkdir plots/radial_velocity_profile
python velocity_profile.py 2. 200 300
python test_energy_conservation.py
echo "Plotting energy as a function of time"
python test_energy_conservation.py 300
......@@ -3,7 +3,7 @@ import h5py as h5
import matplotlib.pyplot as plt
import sys
n_snaps = 5
n_snaps = int(sys.argv[1])
#some constants
OMEGA = 0.3 # Cosmological matter fraction at z = 0
......@@ -24,7 +24,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
......
......@@ -65,6 +65,8 @@ EXTRA_DIST = BigCosmoVolume/makeIC.py \
EAGLE_50/eagle_50.yml EAGLE_50/getIC.sh EAGLE_50/README EAGLE_50/run.sh \
ExternalPointMass/externalPointMass.yml ExternalPointMass/makeIC.py ExternalPointMass/run.sh ExternalPointMass/test.pro \
GreshoVortex_2D/getGlass.sh GreshoVortex_2D/gresho.yml GreshoVortex_2D/makeIC.py GreshoVortex_2D/plotSolution.py GreshoVortex_2D/run.sh \
HydrostaticHalo/README HydrostaticHalo/hydrostatic.yml HydrostaticHalo/makeIC.py HydrostaticHalo/run.sh \
HydrostaticHalo/density_profile.py HydrostaticHalo/velocity_profile.py HydrostaticHalo/internal_energy_profile.py HydrostaticHalo/test_energy_conservation.py
KelvinHelmholtz_2D/kelvinHelmholtz.yml KelvinHelmholtz_2D/makeIC.py KelvinHelmholtz_2D/plotSolution.py KelvinHelmholtz_2D/run.sh \
MultiTypes/makeIC.py MultiTypes/multiTypes.yml MultiTypes/run.sh \
PerturbedBox_2D/makeIC.py PerturbedBox_2D/perturbedPlane.yml \
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment