Added radiative cooling to energy conversion test

examples/CoolingBox/coolingBox.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   4.    # 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_min:     1e-4  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots
@@ -36,7 +36,7 @@ InitialConditions:
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
-  lambda_cgs:                  0.0    # Cooling rate (in cgs units)
+  lambda_cgs:                 1.0e-22    # Cooling rate (in cgs units)
   minimum_temperature:         1.0e4  # Minimal temperature (Kelvin)
   mean_molecular_weight:       0.59   # Mean molecular weight
   hydrogen_mass_abundance:     0.75   # Hydrogen mass abundance (dimensionless)

examples/CoolingBox/energy_plot.py

@@ -11,8 +11,8 @@ snap_filename = "coolingBox_000.hdf5"
 k_b = 1.38E-16 #boltzmann
 m_p = 1.67e-24 #proton mass
 #initial conditions set in makeIC.py
-rho = 3.2e2
-P = 4.5e5
+rho = 3.2e3
+P = 4.5e6
 #n_H_cgs = 0.0001
 gamma = 5./3.
 T_init = 1.0e5
@@ -32,51 +32,65 @@ X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
 #get number of particles
 header = f["Header"]
 n_particles = header.attrs["NumPart_ThisFile"][0]
+
 #read energy and time arrays
 array = np.genfromtxt(stats_filename,skip_header = 1)
 time = array[:,0]
-total_energy = array[:,2]
+kin_plus_therm = array[:,2]
+radiated = array[:,6]
 total_mass = array[:,1]
 
+#ignore first row where there are just zeros
 time = time[1:]
-total_energy = total_energy[1:]
+kin_plus_therm = kin_plus_therm[1:]
+radiated = radiated[1:]
 total_mass = total_mass[1:]
 
+total_energy = kin_plus_therm + radiated
+initial_energy = total_energy[0]
 #conversions to cgs
 rho_cgs = rho * unit_mass / (unit_length)**3
 time_cgs = time * unit_time
-u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2 
+initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2 
 n_H_cgs = X_H * rho_cgs / m_p
+
 #find the energy floor
 u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
+
 #find analytic solution
-analytic_time_cgs = np.linspace(time_cgs[0],time_cgs[-1],1000)
+analytic_time_cgs = np.linspace(0,time_cgs[-1],1000)
 du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
-u_analytic = du_dt_cgs*(analytic_time_cgs - analytic_time_cgs[0]) + u_init_cgs
-cooling_time = u_init_cgs/(-du_dt_cgs)
+u_analytic_cgs = du_dt_cgs*analytic_time_cgs + initial_energy_cgs
+cooling_time_cgs = initial_energy_cgs/(-du_dt_cgs)
+
+for i in range(u_analytic_cgs.size):
+    if u_analytic_cgs[i]<u_floor_cgs:
+        u_analytic_cgs[i] = u_floor_cgs
+
+#rescale analytic solution
+u_analytic = u_analytic_cgs/initial_energy_cgs
 
 #put time in units of cooling_time
-time=time_cgs/cooling_time
-analytic_time = analytic_time_cgs/cooling_time
-#rescale energy to initial energy
-total_energy /= total_energy[0]
-u_analytic /= u_init_cgs
-u_floor_cgs /= u_init_cgs
-# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
-# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
-#analytic_solution = np.zeros(n_snaps-1)
-for i in range(u_analytic.size):
-    if u_analytic[i]<u_floor_cgs:
-        u_analytic[i] = u_floor_cgs
-plt.plot(time-time[0],total_energy,'kd',label = "Numerical solution")
-plt.plot(analytic_time-analytic_time[0],u_analytic,'r',lw = 2.0,label = "Analytic Solution")
+time=time_cgs/cooling_time_cgs
+analytic_time = analytic_time_cgs/cooling_time_cgs
+
+#rescale (numerical) energy by initial energy
+radiated /= initial_energy
+kin_plus_therm /= initial_energy
+total_energy = kin_plus_therm + radiated
+plt.plot(time,kin_plus_therm,'kd',label = "Kinetic + thermal energy")
+plt.plot(time,radiated,'bo',label = "Radiated energy")
+plt.plot(time,total_energy,'g',label = "Total energy")
+plt.plot(analytic_time,u_analytic,'r',lw = 2.0,label = "Analytic Solution")
+#plt.plot(analytic_time,1-u_analytic,'k',lw = 2.0)
 #plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
 #plt.plot((time[1]-time_cgs[0],time_cgs[1]-time_cgs[0]),(0,1),'m',label = "First output")
-plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
-plt.xlabel("Time (cooling time)")
-plt.ylabel("Energy/Initial energy")
-plt.ylim(0,1)
-plt.xlim(0,min(10,time[-1]))
+#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
+plt.xlabel("Time / cooling time")
+plt.ylabel("Energy / Initial energy")
+#plt.ylim(0,1.1)
+plt.ylim(0.999,1.001)
+#plt.xlim(0,min(10,time[-1]))
 plt.legend(loc = "upper right")    
 if (int(sys.argv[1])==0):
     plt.show()

examples/CoolingBox/makeIC.py

@@ -29,8 +29,8 @@ from numpy import *
 periodic= 1           # 1 For periodic box
 boxSize = 1           # 1 kiloparsec    
 L = int(sys.argv[1])  # Number of particles along one axis
-rho = 3.2e6           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 4.5e9             # Pressure in code units (at 10^5K)
+rho = 3.2e3           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+P = 4.5e6             # Pressure in code units (at 10^5K)
 gamma = 5./3.         # Gas adiabatic index
 eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "coolingBox.hdf5" 

examples/CoolingBox/run.sh

@@ -9,6 +9,6 @@ python makeIC.py 10
 
 #-C 2>&1 | tee output.log
 
-#python energy_plot.py 0
+python energy_plot.py 0
 
-python test_energy_conservation.py 0
+#python test_energy_conservation.py 0

src/const.h

@@ -90,14 +90,14 @@
 #define const_gravity_eta 0.025f
 
 /* External gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
-//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+//#define EXTERNAL_POTENTIAL_POINTMASS
+#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 //#define EXTERNAL_POTENTIAL_DISK_PATCH
 
 /* Cooling properties */
-#define COOLING_NONE
+//#define COOLING_NONE
 //#define COOLING_CONST_DU
-//#define COOLING_CONST_LAMBDA
+#define COOLING_CONST_LAMBDA
 //#define COOLING_GRACKLE
 
 /* Are we debugging ? */

src/space.c

@@ -1619,8 +1619,10 @@ void space_init(struct space *s, const struct swift_params *params,
     } else {
       for (size_t k = 0; k < Npart; k++)
         for (int j = 0; j < 3; j++)
-          if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
+          if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j]){
+	    printf("parts[%lld].x[%d] = %f , dim[%d] = %f\n" , k , j , parts[k].x[j] , j , dim[j]);
             error("Not all particles are within the specified domain.");
+	  }
     }
 
     /* Same for the gparts */