diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index 5a907333043951fea0e20748ef6ce4ffc99ec1d7..1484f60596e68734f0f98685ab2ab845f2e0b407 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -27,10 +27,10 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
-boxSize = 1.0
+boxSize = 1.
L = int(sys.argv[1]) # Number of particles along one axis
-rho = 3.3e6 # Density in solar masses per cubic kiloparsec (0.1 hydrogen atoms per cm^3)
-P = 1.4e8 # Pressure in code units (at 10^5K)
+rho = 2. # Density
+P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "uniformBox.hdf5"
@@ -38,7 +38,6 @@ fileName = "uniformBox.hdf5"
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
-print mass
internalEnergy = P / ((gamma - 1.)*rho)
#--------------------------------------------------
diff --git a/examples/UniformBox/plots/energy_plot.py b/examples/UniformBox/plots/energy_plot.py
deleted file mode 100644
index ccfecfc92c40c5c7a064764892f5c181aa63bfe5..0000000000000000000000000000000000000000
--- a/examples/UniformBox/plots/energy_plot.py
+++ /dev/null
@@ -1,49 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-import h5py as h5
-
-n_snaps = 101
-directory = '../'
-snap_name_base = "uniformBox_"
-plot_dir = "./"
-cooling_lambda = 1.0
-u_floor = 0.0
-delta_time = 0.1
-plot_title = r"$\Lambda \, = \, %1.3f \, \, \, u_{floor} = %1.3f$" %(cooling_lambda,u_floor)
-plot_filename = "energy_plot_creasey_lambda_1p0eminus50.png"
-e_kin_array=[]
-e_therm_array=[]
-e_total_array = []
-time = np.zeros(n_snaps-1)
-analytic_solution = np.zeros(n_snaps-1)
-for i in range(1,n_snaps):
- snap_number = "%03d" %i
- filename = directory + snap_name_base + snap_number + ".hdf5"
- f = h5.File(filename)
- print "Reading snap number %d" %i
- header = f["Header"]
- n_parts = header.attrs["NumPart_ThisFile"][0]
- time[i-1] = i*delta_time
- u = np.array(f["PartType0/InternalEnergy"])
- u_total = np.sum(u)
- e_therm_array = np.append(e_therm_array,u_total)
- v = np.array(f["PartType0/Velocities"])
- e_kin_particles = 0.5*v**2
- e_kin_total = np.sum(e_kin_particles)
- e_kin_array = np.append(e_kin_array,e_kin_total)
- e_total_array = np.append(e_total_array,e_kin_total+u_total)
- analytic_solution[i-1] = e_total_array[0] - n_parts*(time[i-1] - time[0])*cooling_lambda
-
-
-# plt.plot(e_kin_array,label = "Kinetic Energy")
-# plt.plot(e_therm_array,label = "Internal Energy")
-plt.plot(time,e_total_array,'r',label = "Total Energy")
-plt.plot(time,analytic_solution,'--k',label = "Analytic Solution")
-plt.xlabel("Time (seconds)")
-plt.ylabel("Energy (erg/g)")
-plt.ylim((0,1000))
-plt.legend(loc = "upper right")
-plt.title(plot_title)
-full_plot_filename = "%s/%s" %(plot_dir,plot_filename)
-plt.savefig(full_plot_filename,format = "png")
-plt.close()
diff --git a/examples/UniformBox/plots/plot_particle_positions.py b/examples/UniformBox/plots/plot_particle_positions.py
deleted file mode 100644
index 97c8cca439c973b219aef5cbff4673cae85d6b53..0000000000000000000000000000000000000000
--- a/examples/UniformBox/plots/plot_particle_positions.py
+++ /dev/null
@@ -1,56 +0,0 @@
-import h5py as h5
-import numpy as np
-import matplotlib.pyplot as plt
-
-n_snaps = 100
-directory = '../'
-snap_name_base = "uniformBox_"
-plot_directory = "./"
-plot_filename_base = "snap"
-for i in range(n_snaps):
- snap_number = "%03d" %i
- filename = directory + snap_name_base + snap_number + ".hdf5"
- f = h5.File(filename)
- #find the box size
-
- header = f["Header"]
- box_size = header.attrs["BoxSize"][0]
- pos = np.array(f["PartType0/Coordinates"])
- x = pos[:,0]
- y = pos[:,1]
- z = pos[:,2]
-
- #x_projection
- plt.plot(y,z,'bo')
- plt.xlabel("y")
- plt.ylabel("z")
- plt.xlim((0.0,box_size))
- plt.ylim((0.0,box_size))
- plot_dir = "%s/projection_x" %plot_directory
- plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
- plt.savefig(plot_filename,format = "png")
- plt.close()
-
- #y_projection
- plt.plot(x,z,'bo')
- plt.xlabel("x")
- plt.ylabel("z")
- plt.xlim((0.0,box_size))
- plt.ylim((0.0,box_size))
- plot_dir = "%s/projection_y" %plot_directory
- plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
- plt.savefig(plot_filename,format = "png")
- plt.close()
-
- #x_projection
- plt.plot(x,y,'bo')
- plt.xlabel("x")
- plt.ylabel("y")
- plt.xlim((0.0,box_size))
- plt.ylim((0.0,box_size))
- plot_dir = "%s/projection_z" %plot_directory
- plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
- plt.savefig(plot_filename,format = "png")
- plt.close()
-
-
diff --git a/examples/UniformBox/run.sh b/examples/UniformBox/run.sh
index 36ef1d6af1cb97de659d04a8d3cd9090b027ac15..0cb0a505915be47bafb99ed7531685bfeb3dc829 100755
--- a/examples/UniformBox/run.sh
+++ b/examples/UniformBox/run.sh
@@ -4,7 +4,7 @@
if [ ! -e uniformBox.hdf5 ]
then
echo "Generating initial conditions for the uniform box example..."
- python makeIC.py 10
+ python makeIC.py 100
fi
-../swift -s -C -t 16 uniformBox.yml
+../swift -s -t 16 uniformBox.yml
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
index 670bbd483b9fa68f56e26a675df7134114f807b6..7c9c74e1342bffb939131a265188cae269cc773f 100644
--- a/examples/UniformBox/uniformBox.yml
+++ b/examples/UniformBox/uniformBox.yml
@@ -1,15 +1,15 @@
# Define the system of units to use internally.
InternalUnitSystem:
- UnitMass_in_cgs: 2.0e33 # Solar masses
- UnitLength_in_cgs: 3.01e21 # Kilparsecs
- UnitVelocity_in_cgs: 1.0e5 # Kilometres per second
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # 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: 10. # The end time of the simulation (in internal units).
+ time_end: 1. # 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).
@@ -17,7 +17,7 @@ TimeIntegration:
Snapshots:
basename: uniformBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
+ delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
@@ -41,12 +41,3 @@ PointMass:
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass in internal units
-
- # Cooling parameters (Creasey cooling) (always in cgs)
-
-CreaseyCooling:
- Lambda: 1.0e-22
- minimum_temperature: 1.0e4
- mean_molecular_weight: 0.59
- hydrogen_mass_abundance: 0.75
- cooling_tstep_mult: 1.0
\ No newline at end of file
diff --git a/src/cooling.c b/src/cooling.c
index 3edb7fcf7af8f4ec83ad16594dc887f7dd176e61..6e4aa6314b7b51030faf16ac78d9fb6d9b39fce1 100644
--- a/src/cooling.c
+++ b/src/cooling.c
@@ -133,7 +133,7 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
float u_new;
float m_p = phys_const->const_proton_mass;
float X_H = cooling->creasey_cooling.hydrogen_mass_abundance;
- float lambda = cooling->creasey_cooling.lambda; //this is always in cgs
+ float lambda_cgs = cooling->creasey_cooling.lambda; //this is always in cgs
float u_floor = cooling->creasey_cooling.min_internal_energy;
/*convert from internal code units to cgs*/
@@ -143,19 +143,19 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
float n_H_cgs = rho_cgs / m_p_cgs;
float u_old_cgs = u_old * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
float u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
- float du_dt = -lambda * n_H_cgs * n_H_cgs / rho;
+ float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
float u_new_cgs;
- if (u_old_cgs - du_dt*dt_cgs > u_floor_cgs){
- u_new_cgs = u_old_cgs + du_dt*dt_cgs;
+ if (u_old_cgs + du_dt_cgs*dt_cgs > u_floor_cgs){
+ u_new_cgs = u_old_cgs + du_dt_cgs*dt_cgs;
}
else{
u_new_cgs = u_floor_cgs;
}
-
/*convert back to internal code units when returning new internal energy*/
u_new = u_new_cgs / units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
+
#endif /*CREASEY_COOLING*/
return u_new;
}