diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index 1501c1e1dc00987dee7bdf25b008f393affe8cff..01b958860257cb714f551608b6f91bd2dbf46e95 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/examples/CoolingHaloWithSpin/cooling_halo.yml
@@ -40,7 +40,7 @@ IsothermalPotential:
position_z: 0.
vrot: 200. # Rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # Controls time step
- epsilon: 1.0 # Softening for the isothermal potential
+ epsilon: 0.1 # Softening for the isothermal potential
# Cooling parameters
LambdaCooling:
@@ -48,4 +48,4 @@ LambdaCooling:
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
- cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
+ cooling_tstep_mult: 0.1 # Dimensionless pre-factor for the time-step condition
diff --git a/examples/CoolingHaloWithSpin/density_profile.py b/examples/CoolingHaloWithSpin/density_profile.py
index ea282328e5b75530a128eab2dec5f065e46cf819..fb88ddd6aea71603a6f6fcb36b13771106737e6a 100644
--- a/examples/CoolingHaloWithSpin/density_profile.py
+++ b/examples/CoolingHaloWithSpin/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["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["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"])
@@ -101,18 +101,18 @@ for i in range(n_snaps):
rho_0 = density[0]
rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
- plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
- #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
+ plt.plot(radial_bin_mids,density,'ko',label = "Average density of shell")
+ plt.plot(radial_bin_mids,rho_analytic_init,label = "Initial analytic density profile")
plt.xlabel(r"$r / r_{vir}$")
- plt.ylabel(r"$\rho / \rho_{init})$")
+ plt.ylabel(r"$(\rho / \rho_{init})$")
plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
#plt.ylim((1.e-2,1.e1))
- plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
+ plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(1.0e-4,1.0e4),'r',label = "Cooling radius")
plt.xlim((radial_bin_mids[0],max_r))
- plt.ylim((0,20))
- plt.plot((0,max_r),(1,1))
+ plt.ylim((1.0e-4,1.0e4))
+ #plt.plot((0,max_r),(1,1))
#plt.xscale('log')
- #plt.yscale('log')
+ plt.yscale('log')
plt.legend(loc = "upper right")
plot_filename = "density_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
diff --git a/examples/CoolingHaloWithSpin/internal_energy_profile.py b/examples/CoolingHaloWithSpin/internal_energy_profile.py
index a3e470cc24a939c9bc915371e927d9bd39196bff..5f71d69ca7a978de242559f84ec390faa86a27f0 100644
--- a/examples/CoolingHaloWithSpin/internal_energy_profile.py
+++ b/examples/CoolingHaloWithSpin/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["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["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"])
diff --git a/examples/CoolingHaloWithSpin/makeIC.py b/examples/CoolingHaloWithSpin/makeIC.py
index 8970fbaa70578532a4f41bab7a096d8fa3565d26..a6d57868ad7542498b27007a5c3ef9234b9feb84 100644
--- a/examples/CoolingHaloWithSpin/makeIC.py
+++ b/examples/CoolingHaloWithSpin/makeIC.py
@@ -36,6 +36,7 @@ h = 0.67777 # hubble parameter
gamma = 5./3.
eta = 1.2349
spin_lambda = 0.05 #spin parameter
+f_b = 0.2 #baryon fraction
# First set unit velocity and then the circular velocity parameter for the isothermal potential
const_unit_velocity_in_cgs = 1.e5 #kms^-1
@@ -99,6 +100,8 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
# set seed for random number
np.random.seed(1234)
+gas_mass = f_b * np.sqrt(3.) / 2. #virial mass of halo is 1, virial radius is 1, enclosed mass scales with r
+gas_particle_mass = gas_mass / float(N)
# Positions
# r^(-2) distribution corresponds to uniform distribution in radius
@@ -164,12 +167,12 @@ N = x_coords.size
print "Number of particles in the box = " , N
#make the coords and radius arrays again
-coords_2 = np.zeros((N,3))
-coords_2[:,0] = x_coords
-coords_2[:,1] = y_coords
-coords_2[:,2] = z_coords
+coords= np.zeros((N,3))
+coords[:,0] = x_coords
+coords[:,1] = y_coords
+coords[:,2] = z_coords
-radius = np.sqrt((coords_2[:,0]-boxSize/2.)**2 + (coords_2[:,1]-boxSize/2.)**2 + (coords_2[:,2]-boxSize/2.)**2)
+radius = np.sqrt((coords[:,0]-boxSize/2.)**2 + (coords[:,1]-boxSize/2.)**2 + (coords[:,2]-boxSize/2.)**2)
#now give particle's velocities
v = np.zeros((N,3))
@@ -184,7 +187,7 @@ print "J =", J
omega = np.zeros((N,3))
for i in range(N):
omega[i,2] = 3.*J / radius[i]
- v[i,:] = np.cross(omega[i,:],(coords_2[i,:]-boxSize/2.))
+ v[i,:] = np.cross(omega[i,:],(coords[i,:]-boxSize/2.))
# Header
grp = file.create_group("/Header")
@@ -202,16 +205,15 @@ grp.attrs["Dimension"] = 3
grp = file.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (N, 3), 'd')
-ds[()] = coords_2
-coords_2 = np.zeros(1)
+ds[()] = coords
+coords = np.zeros(1)
ds = grp.create_dataset('Velocities', (N, 3), 'f')
ds[()] = v
v = np.zeros(1)
# All particles of equal mass
-mass = 1. / N
-m = np.full((N,),mass)
+m = np.full((N,),gas_particle_mass)
ds = grp.create_dataset('Masses', (N, ), 'f')
ds[()] = m
m = np.zeros(1)
diff --git a/examples/CoolingHaloWithSpin/velocity_profile.py b/examples/CoolingHaloWithSpin/velocity_profile.py
index d64d255b18482bc26578f21f46199aa3540ae7b5..07df8e1b0751307513c30a5b128773b193c3a9cd 100644
--- a/examples/CoolingHaloWithSpin/velocity_profile.py
+++ b/examples/CoolingHaloWithSpin/velocity_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["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index ae13b49c9561096b11d021f710c687d153f2d269..fb51ea0af1b04125079c0788104c82e2e59577e9 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -89,8 +89,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
float cooling_du_dt = cooling_rate(phys_const, us, cooling, p);
/* Integrate cooling equation to enforce energy floor */
- if (u_old + (hydro_du_dt + cooling_du_dt) * dt < u_floor) {
- cooling_du_dt = -(u_old + dt * hydro_du_dt - u_floor) / dt;
+ /* Factor of 1.5 included since timestep could potentially double */
+ if (u_old + (hydro_du_dt + cooling_du_dt) * 1.5f * dt < u_floor) {
+ cooling_du_dt = -(u_old + 1.5f *dt * hydro_du_dt - u_floor) /(1.5f * dt);
}
/* Update the internal energy time derivative */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 202099839227dffb0244dee09af0c87be14d6ab2..3c49f8a288e37859a088b0852e983a72c5f80017 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -332,11 +332,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->rho *= expf(w2);
- if(p->entropy + p->entropy_dt*dt < 0.0)
- printf("Negative entropy: Old entropy = %g, New entropy = %g ", p->entropy,p->entropy + p->entropy_dt*dt);
/* Predict the entropy */
p->entropy += p->entropy_dt * dt;
-
/* Re-compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
@@ -379,7 +376,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt) {
/* Do not decrease the entropy by more than a factor of 2 */
- p->entropy_dt = max(-0.5f * xp->entropy_full / dt , p->entropy_dt);
+
+ if (p->entropy_dt < -0.5f * xp->entropy_full / dt){
+ message("Warning! Limiting entropy_dt. Possible cooling error.\n entropy_full = %g \n entropy_dt * dt =%g \n",
+ xp->entropy_full,p->entropy_dt * dt);
+ p->entropy_dt = -0.5f * xp->entropy_full / dt;
+ }
xp->entropy_full += p->entropy_dt * dt;
/* Compute the pressure */
@@ -409,7 +411,6 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* We read u in the entropy field. We now get S from u */
xp->entropy_full = gas_entropy_from_internal_energy(p->rho, p->entropy);
p->entropy = xp->entropy_full;
-
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);