Commit 22970bd2 authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Cooling actually works now

parent 6e2c2eac
...@@ -2,46 +2,40 @@ import numpy as np ...@@ -2,46 +2,40 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import h5py as h5 import h5py as h5
n_snaps = 102 n_snaps = 101
directory = '../' directory = '../'
snap_name_base = "Isothermal_" snap_name_base = "uniformBox_"
plot_dir = "random_orbits_1d" plot_dir = "./"
plot_filename_base = "snap" plot_title = r"$\Lambda \, = \, 1.0$"
plot_filename = "energy_plot_lambda_1.png"
v_rot = 200.
r_0 = 100.
e_kin_array = [] e_kin_array = []
e_pot_array = [] e_therm_array = []
e_total_array = [] e_total_array = []
snap_numbers = np.linspace(0,n_snaps) snap_numbers = np.linspace(0,n_snaps)
for i in range(n_snaps): for i in range(n_snaps):
snap_number = "%03d" %i snap_number = "%03d" %i
filename = directory + snap_name_base + snap_number + ".hdf5" filename = directory + snap_name_base + snap_number + ".hdf5"
f = h5.File(filename) f = h5.File(filename)
print "Reading snap number %d" %i
header = f["Header"] header = f["Header"]
n_parts = header.attrs["NumPart_ThisFile"][1] n_parts = header.attrs["NumPart_ThisFile"][1]
pos = np.array(f["PartType1/Coordinates"]) u = np.array(f["PartType0/InternalEnergy"])
x = pos[:,0] u_total = np.sum(u)
y = pos[:,1] e_therm_array = np.append(e_therm_array,u_total)
z = pos[:,2] v = np.array(f["PartType0/Velocities"])
r = np.sqrt(x**2 + y**2 + z**2)
e_pot_particles = v_rot**2 * np.log(r/r_0)
e_pot_total = np.sum(e_pot_particles)
e_pot_array = np.append(e_pot_array,e_pot_total)
v = np.array(f["PartType1/Velocities"])
e_kin_particles = 0.5*v**2 e_kin_particles = 0.5*v**2
e_kin_total = np.sum(e_kin_particles) e_kin_total = np.sum(e_kin_particles)
e_kin_array = np.append(e_kin_array,e_kin_total) e_kin_array = np.append(e_kin_array,e_kin_total)
e_total_array = np.append(e_total_array,e_kin_total+e_pot_total) e_total_array = np.append(e_total_array,e_kin_total+u_total)
print snap_numbers.shape
print e_kin_array.shape
plt.plot(e_kin_array,label = "Kinetic Energy") plt.plot(e_kin_array,label = "Kinetic Energy")
plt.plot(e_pot_array,label = "Potential Energy") plt.plot(e_therm_array,label = "Internal Energy")
plt.plot(e_total_array,label = "Total Energy") plt.plot(e_total_array,label = "Total Energy")
plt.legend(loc = "upper right") plt.xlabel("Snapshot number")
plot_filename = "%s/energy_plot.png" %plot_dir plt.ylabel("Energy (code units)")
plt.savefig(plot_filename,format = "png") plt.legend(loc = "lower right")
plt.title(plot_title)
full_plot_filename = "%s/%s" %(plot_dir,plot_filename)
plt.savefig(full_plot_filename,format = "png")
plt.close() plt.close()
...@@ -41,3 +41,10 @@ PointMass: ...@@ -41,3 +41,10 @@ PointMass:
position_y: 50. position_y: 50.
position_z: 50. position_z: 50.
mass: 1e10 # mass of external point mass in internal units mass: 1e10 # mass of external point mass in internal units
# Cooling parameters (Constant cooling)
Cooling:
lambda: 0.1
min_energy: 0.0
cooling_tstep_mult: 1.0
\ No newline at end of file
...@@ -88,11 +88,11 @@ float calculate_new_thermal_energy(float u_old, float dt, const struct cooling_d ...@@ -88,11 +88,11 @@ float calculate_new_thermal_energy(float u_old, float dt, const struct cooling_d
//This function integrates the cooling equation, given the initial thermal energy and the timestep dt. //This function integrates the cooling equation, given the initial thermal energy and the timestep dt.
//Returns 0 if successful and 1 if not //Returns 0 if successful and 1 if not
int status = 0; int status = 0;
float du_dt = cooling->const_cooling.lambda; float du_dt = -cooling->const_cooling.lambda;
float u_floor = cooling->const_cooling.min_energy; float u_floor = cooling->const_cooling.min_energy;
float u_new; float u_new;
if (u_old - du_dt*dt > u_floor){ if (u_old - du_dt*dt > u_floor){
u_new = u_old - du_dt*dt; u_new = u_old + du_dt*dt;
} }
else{ else{
u_new = u_floor; u_new = u_floor;
......
...@@ -167,14 +167,14 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) { ...@@ -167,14 +167,14 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
OUT; OUT;
#endif #endif
/* Loop over the gparts in this cell. */ /* Loop over the parts in this cell. */
for (int i = 0; i < count; i++) { for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */ /* Get a direct pointer on the part. */
struct part *restrict p = &parts[i]; struct part *restrict p = &parts[i];
/* Is this part within the time step? */ /* Kick has already updated ti_end, so need to check ti_begin */
if (p->ti_end <= ti_current) { if (p->ti_begin == ti_current) {
dt = (p->ti_end - p->ti_begin)*timeBase; dt = (p->ti_end - p->ti_begin)*timeBase;
update_entropy(cooling, constants, p, dt); update_entropy(cooling, constants, p, dt);
} }
......
...@@ -101,6 +101,7 @@ enum task_actions task_acts_on(const struct task *t) { ...@@ -101,6 +101,7 @@ enum task_actions task_acts_on(const struct task *t) {
case task_type_sort: case task_type_sort:
case task_type_ghost: case task_type_ghost:
case task_type_cooling:
return task_action_part; return task_action_part;
break; break;
......
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