diff --git a/examples/CoolingHalo/run.sh b/examples/CoolingHalo/run.sh
index 3a0d9c02000e760b030a96107038d3c6163f3227..60ceae649d183dce3a7e5019a1ff94ce7bc4f08d 100755
--- a/examples/CoolingHalo/run.sh
+++ b/examples/CoolingHalo/run.sh
@@ -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 100
-# python internal_energy_profile.py 10
+python internal_energy_profile.py 2. 200 100
-# python test_energy_conservation.py
+python test_energy_conservation.py 2. 200 100
diff --git a/examples/HydrostaticHalo/README b/examples/HydrostaticHalo/README
index 8644cf0c7bb3a796246edfa70936a99b44770061..fefc6a910a9afb8e57d39dc2309c5d0e6d268b50 100644
--- a/examples/HydrostaticHalo/README
+++ b/examples/HydrostaticHalo/README
@@ -1,6 +1,21 @@
;
-;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
diff --git a/examples/HydrostaticHalo/radial_profile.py b/examples/HydrostaticHalo/density_profile.py
similarity index 96%
rename from examples/HydrostaticHalo/radial_profile.py
rename to examples/HydrostaticHalo/density_profile.py
index fd37df88048fbab997d6f3893254371c1c9ae022..52bebb9ffefa77dae66af155fb31fed539dcde13 100644
--- a/examples/HydrostaticHalo/radial_profile.py
+++ b/examples/HydrostaticHalo/density_profile.py
@@ -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["SoftenedIsothermalPotential: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()
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index 601210deac10350bee475b6935745cef7b7e4f23..21d1a79dd1a63a3a2ce3ff300632d6c79cb9a973 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -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: 10.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:
+SoftenedIsothermalPotential:
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
diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py
index 65ba231ed7732bd386f663a364c053b5ec9a124e..a1b2bda314a66eb965974d34519f66c544ee8aed 100644
--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
@@ -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["SoftenedIsothermalPotential: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()
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index aa3f2ae3074c88c0b462b1b0185e355aae319cc5..f3d3694760ced9e41719de64229f0b5565adf5d4 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -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 100
-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
diff --git a/examples/HydrostaticHalo/test_energy_conservation.py b/examples/HydrostaticHalo/test_energy_conservation.py
index 6fd0a9b411528105723a9223d73ffed09dcbf5c0..26041890c0e5c8b80abda4448f0ed55d8e3b98d4 100644
--- a/examples/HydrostaticHalo/test_energy_conservation.py
+++ b/examples/HydrostaticHalo/test_energy_conservation.py
@@ -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')
diff --git a/examples/HydrostaticHalo/velocity_profile.py b/examples/HydrostaticHalo/velocity_profile.py
new file mode 100644
index 0000000000000000000000000000000000000000..f6a7350b9731d660b2092266d4d6ad3730bab48c
--- /dev/null
+++ b/examples/HydrostaticHalo/velocity_profile.py
@@ -0,0 +1,111 @@
+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()
diff --git a/src/const.h b/src/const.h
index 3fba98c02c9d68b32cb99f89655798551ebe63bb..333274ee1502a4380473dd0810a9b15940f833c4 100644
--- a/src/const.h
+++ b/src/const.h
@@ -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 */
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
index 60c5b230bc2229eff38384fc3d56f3161cada2c4..fe1df8796f046edded0c5b1779859a1c6fffffc0 100644
--- a/src/potential/disc_patch/potential.h
+++ b/src/potential/disc_patch/potential.h
@@ -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 gpart* p) {
return 0.f;
}
diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h
index a17932c06d7b5784c94da4f06aa2394307a5740b..8248b64678e28e06b9df4aab375cde0b5ed5281b 100644
--- a/src/potential/none/potential.h
+++ b/src/potential/none/potential.h
@@ -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 gpart* g) {
return 0.f;
}
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index 16aaf1330daee68b8746feb188ba3a35ef543e1d..7b8092d666168d3e207fe5cf94dc08b90b56e5d3 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -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 gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
diff --git a/src/potential/softened_isothermal/potential.h b/src/potential/softened_isothermal/potential.h
index b28e5f7a17e0c7fd1956e099a8d4e21a5ec7ff20..24e59b12a5745728fb1189fbbfbc7cc3c06fbfa6 100644
--- a/src/potential/softened_isothermal/potential.h
+++ b/src/potential/softened_isothermal/potential.h
@@ -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 gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;