Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
f83f665c
Commit
f83f665c
authored
Oct 25, 2016
by
Stefan Arridge
Browse files
Made changes to Hydrostatic halo scripts and included a README file
parent
6578b2fa
Changes
13
Hide whitespace changes
Inline
Side-by-side
examples/CoolingHalo/run.sh
View file @
f83f665c
...
...
@@ -6,8 +6,8 @@ python makeIC.py 10000
../swift
-g
-s
-C
-t
16 cooling_halo.yml 2>&1 |
tee
output.log
#
python radial_profile.py 10
python radial_profile.py
2. 200
10
0
#
python internal_energy_profile.py 10
python internal_energy_profile.py
2. 200
10
0
#
python test_energy_conservation.py
python test_energy_conservation.py
2. 200 100
examples/HydrostaticHalo/README
View file @
f83f665c
;
;In this example we distribute a set of gas particles in sphere with an r^-2
;density profile, with an external isothermal potential. The particles are
;stationary with a pressure gradient such that we have hydrostatic equilibrium.
;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 mass and radius of the halo is set by 'v_rot'
\ No newline at end of file
;The parameter v_rot (in makeIC.py and cooling.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.
;
;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.
;
;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 with time. These can be used to check if the example
;has run properly.
;
\ No newline at end of file
examples/HydrostaticHalo/
radial
_profile.py
→
examples/HydrostaticHalo/
density
_profile.py
View file @
f83f665c
...
...
@@ -28,7 +28,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
[
"
Softened
IsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...
...
@@ -114,7 +114,7 @@ for i in range(n_snaps):
#plt.xscale('log')
#plt.yscale('log')
plt
.
legend
(
loc
=
"upper right"
)
plot_filename
=
"density_profile_%03d.png"
%
i
plot_filename
=
"
./plots/density_profile/
density_profile_%03d.png"
%
i
plt
.
savefig
(
plot_filename
,
format
=
"png"
)
plt
.
close
()
examples/HydrostaticHalo/hydrostatic.yml
View file @
f83f665c
...
...
@@ -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
.0
# The end time of the simulation (in internal units).
time_end
:
1.0
# 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).
...
...
@@ -38,10 +38,11 @@ InitialConditions:
shift_z
:
0.
# External potential parameters
IsothermalPotential
:
Softened
IsothermalPotential
:
position_x
:
0.
# location of centre of isothermal potential in internal units
position_y
:
0.
position_z
:
0.
vrot
:
200.
# rotation speed of isothermal potential in internal units
epsilon
:
0.1
timestep_mult
:
0.03
# controls time step
examples/HydrostaticHalo/internal_energy_profile.py
View file @
f83f665c
...
...
@@ -46,7 +46,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
[
"
Softened
IsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...
...
@@ -102,7 +102,7 @@ for i in range(n_snaps):
plt
.
ylabel
(
r
"$u / (v_c^2 / (2(\gamma - 1)) $"
)
plt
.
title
(
r
"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$"
%
(
snap_time_cgs
,
N
,
v_c
))
plt
.
ylim
((
0
,
2
))
plot_filename
=
"internal_energy_profile_%03d.png"
%
i
plot_filename
=
"
./plots/internal_energy/
internal_energy_profile_%03d.png"
%
i
plt
.
savefig
(
plot_filename
,
format
=
"png"
)
plt
.
close
()
...
...
examples/HydrostaticHalo/run.sh
View file @
f83f665c
...
...
@@ -6,8 +6,10 @@ python makeIC.py 10000
../swift
-g
-s
-t
16 hydrostatic.yml 2>&1 |
tee
output.log
python
radial
_profile.py 10
python
density
_profile.py
2. 200
10
0
python internal_energy_profile.py 10
python internal_energy_profile.py 2. 200 100
python velocity_profile.py 2. 200 100
python test_energy_conservation.py
examples/HydrostaticHalo/test_energy_conservation.py
View file @
f83f665c
...
...
@@ -91,6 +91,5 @@ plt.xlabel(r"$t / t_{dyn}$")
plt
.
ylabel
(
r
"$E / v_c^2$"
)
plt
.
title
(
r
"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$"
%
(
N
,
v_c
))
plt
.
ylim
((
-
2
,
2
))
#plot_filename = "density_profile_%03d.png" %i
plt
.
show
()
plt
.
savefig
(
"energy_conservation.png"
,
format
=
'png'
)
examples/HydrostaticHalo/velocity_profile.py
0 → 100644
View file @
f83f665c
import
numpy
as
np
import
h5py
as
h5
import
matplotlib.pyplot
as
plt
import
sys
def
do_binning
(
x
,
y
,
x_bin_edges
):
#x and y are arrays, where y = f(x)
#returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
n_bins
=
x_bin_edges
.
size
-
1
count
=
np
.
zeros
(
n_bins
)
y_totals
=
np
.
zeros
(
n_bins
)
for
i
in
range
(
n_bins
):
ind
=
np
.
intersect1d
(
np
.
where
(
x
>
bin_edges
[
i
])[
0
],
np
.
where
(
x
<=
bin_edges
[
i
+
1
])[
0
])
count
[
i
]
=
ind
.
size
binned_y
=
y
[
ind
]
y_totals
[
i
]
=
np
.
sum
(
binned_y
)
return
(
count
,
y_totals
)
#for the plotting
max_r
=
float
(
sys
.
argv
[
1
])
n_radial_bins
=
int
(
sys
.
argv
[
2
])
n_snaps
=
int
(
sys
.
argv
[
3
])
#some constants
OMEGA
=
0.3
# Cosmological matter fraction at z = 0
PARSEC_IN_CGS
=
3.0856776e18
KM_PER_SEC_IN_CGS
=
1.0e5
CONST_G_CGS
=
6.672e-8
CONST_m_H_CGS
=
1.67e-24
h
=
0.67777
# hubble parameter
gamma
=
5.
/
3.
eta
=
1.2349
H_0_cgs
=
100.
*
h
*
KM_PER_SEC_IN_CGS
/
(
1.0e6
*
PARSEC_IN_CGS
)
#read some header/parameter information from the first snapshot
filename
=
"Hydrostatic_000.hdf5"
f
=
h5
.
File
(
filename
,
'r'
)
params
=
f
[
"Parameters"
]
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
[
"SoftenedIsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
header
=
f
[
"Header"
]
N
=
header
.
attrs
[
"NumPart_Total"
][
0
]
box_centre
=
np
.
array
(
header
.
attrs
[
"BoxSize"
])
#calculate r_vir and M_vir from v_c
r_vir_cgs
=
v_c_cgs
/
(
10.
*
H_0_cgs
*
np
.
sqrt
(
OMEGA
))
M_vir_cgs
=
r_vir_cgs
*
v_c_cgs
**
2
/
CONST_G_CGS
for
i
in
range
(
n_snaps
):
filename
=
"Hydrostatic_%03d.hdf5"
%
i
f
=
h5
.
File
(
filename
,
'r'
)
coords_dset
=
f
[
"PartType0/Coordinates"
]
coords
=
np
.
array
(
coords_dset
)
#translate coords by centre of box
header
=
f
[
"Header"
]
snap_time
=
header
.
attrs
[
"Time"
]
snap_time_cgs
=
snap_time
*
unit_time_cgs
coords
[:,
0
]
-=
box_centre
[
0
]
/
2.
coords
[:,
1
]
-=
box_centre
[
1
]
/
2.
coords
[:,
2
]
-=
box_centre
[
2
]
/
2.
radius
=
np
.
sqrt
(
coords
[:,
0
]
**
2
+
coords
[:,
1
]
**
2
+
coords
[:,
2
]
**
2
)
radius_cgs
=
radius
*
unit_length_cgs
radius_over_virial_radius
=
radius_cgs
/
r_vir_cgs
#get the internal energies
vel_dset
=
f
[
"PartType0/Velocities"
]
vel
=
np
.
array
(
vel_dset
)
#make dimensionless
vel
/=
v_c
r
=
radius_over_virial_radius
#find radial component of velocity
v_r
=
np
.
zeros
(
r
.
size
)
for
j
in
range
(
r
.
size
):
v_r
[
j
]
=
-
np
.
dot
(
coords
[
j
,:],
vel
[
j
,:])
/
radius
[
j
]
bin_edges
=
np
.
linspace
(
0
,
max_r
,
n_radial_bins
+
1
)
(
hist
,
v_r_totals
)
=
do_binning
(
r
,
v_r
,
bin_edges
)
bin_widths
=
bin_edges
[
1
]
-
bin_edges
[
0
]
radial_bin_mids
=
np
.
linspace
(
bin_widths
/
2.
,
max_r
-
bin_widths
/
2.
,
n_radial_bins
)
binned_v_r
=
v_r_totals
/
hist
#calculate cooling radius
#r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
plt
.
plot
(
radial_bin_mids
,
binned_v_r
,
'ko'
,
label
=
"Average radial velocity in shell"
)
#plt.plot((0,1),(1,1),label = "Analytic Solution")
#plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
plt
.
legend
(
loc
=
"upper right"
)
plt
.
xlabel
(
r
"$r / r_{vir}$"
)
plt
.
ylabel
(
r
"$v_r / v_c$"
)
plt
.
title
(
r
"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$"
%
(
snap_time_cgs
,
N
,
v_c
))
plt
.
ylim
((
0
,
2
))
plot_filename
=
"./plots/radial_velocity_profile/velocity_profile_%03d.png"
%
i
plt
.
savefig
(
plot_filename
,
format
=
"png"
)
plt
.
close
()
src/const.h
View file @
f83f665c
...
...
@@ -92,10 +92,10 @@
/* External gravity properties */
#define EXTERNAL_POTENTIAL_NONE
//
#define EXTERNAL_POTENTIAL_NONE
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//
#define EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL
#define EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL
//#define EXTERNAL_POTENTIAL_DISC_PATCH
/* Source terms */
...
...
src/potential/disc_patch/potential.h
View file @
f83f665c
...
...
@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_get_potential_energy
(
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
p
)
{
const
struct
phys_const
*
const
phys_const
,
const
struct
g
part
*
p
)
{
return
0
.
f
;
}
...
...
src/potential/none/potential.h
View file @
f83f665c
...
...
@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_get_potential_energy
(
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
g
)
{
const
struct
phys_const
*
const
phys_const
,
const
struct
g
part
*
g
)
{
return
0
.
f
;
}
...
...
src/potential/point_mass/potential.h
View file @
f83f665c
...
...
@@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_get_potential_energy
(
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
g
)
{
const
struct
phys_const
*
const
phys_const
,
const
struct
g
part
*
g
)
{
const
float
dx
=
g
->
x
[
0
]
-
potential
->
x
;
const
float
dy
=
g
->
x
[
1
]
-
potential
->
y
;
...
...
src/potential/softened_isothermal/potential.h
View file @
f83f665c
...
...
@@ -134,7 +134,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_get_potential_energy
(
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
part
*
g
)
{
const
struct
phys_const
*
const
phys_const
,
const
struct
g
part
*
g
)
{
const
float
dx
=
g
->
x
[
0
]
-
potential
->
x
;
const
float
dy
=
g
->
x
[
1
]
-
potential
->
y
;
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment