working on newton's method integration scheme for cooling

examples/CoolingBox/analytical_test.py

@@ -86,7 +86,7 @@ temp_0 = 1.0e7
 rho = rho*unit_mass/(unit_length**3)
 
 # Read snapshots
-nsnap = 150
+nsnap = 110
 npart = 4096
 u_snapshots_cgs = zeros(nsnap)
 u_part_snapshots_cgs = zeros((nsnap,npart))
@@ -137,6 +137,7 @@ p = figure()
 p1, = plot(t_sol,temp_sol,linewidth = 0.5,color = 'k',label = 'analytical')
 p2, = plot(t_snapshots_cgs,mean_temp,linewidth = 0.5,color = 'r',label = 'swift mean')
 l = legend(handles = [p1,p2])
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
 xlabel("${\\rm{Time~[s]}}$", labelpad=0)
 ylabel("${\\rm{Temperature~[K]}}$")
 

examples/CoolingBox/makeIC.py

@@ -27,7 +27,7 @@ from numpy import *
 # Parameters
 periodic= 1           # 1 For periodic box
 boxSize = 1           # 1 kiloparsec    
-rho = 3.2e3           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+rho = 3.2e6           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
 P = 4.5e7             # Pressure in code units (at 10^6K)
 gamma = 5./3.         # Gas adiabatic index
 eta = 1.2349          # 48 ngbs with cubic spline kernel

examples/CoolingBox/run.sh

@@ -21,7 +21,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -C -t 8 coolingBox.yml
+../swift -s -C -t 16 coolingBox.yml
 
 # Check energy conservation and cooling rate
 python energy_plot.py

examples/EAGLE_12/eagle_12.yml

@@ -54,5 +54,22 @@ SPH:
 InitialConditions:
   file_name:  ./EAGLE_ICs_12.hdf5    # The file to read
   cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
-  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+
+EAGLEChemistry:
+  InitMetallicity:         0.014
+  InitAbundance_Hydrogen:  0.70649785
+  InitAbundance_Helium:    0.28055534
+  InitAbundance_Carbon:    2.0665436e-3
+  InitAbundance_Nitrogen:  8.3562563e-4
+  InitAbundance_Oxygen:    5.4926244e-3
+  InitAbundance_Neon:      1.4144605e-3
+  InitAbundance_Magnesium: 5.907064e-4
+  InitAbundance_Silicon:   6.825874e-4
+  InitAbundance_Iron:      1.1032152e-3
+  CalciumOverSilicon:      0.0941736
+  SulphurOverSilicon:      0.6054160
+
+EagleCooling:
+  filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/
+  reionisation_redshift:   8.898
 

examples/EAGLE_12/run.sh

@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -G -S -t 16 eagle_12.yml 2>&1 | tee output.log
+../swift -s -c -C -t 16 eagle_12.yml 2>&1 | tee output.log
 

examples/EAGLE_6/eagle_6.yml

@@ -70,3 +70,40 @@ InitialConditions:
   cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
 
 
+EAGLEChemistry:
+  InitMetallicity:         0.014
+  InitAbundance_Hydrogen:  0.70649785
+  InitAbundance_Helium:    0.28055534
+  InitAbundance_Carbon:    2.0665436e-3
+  InitAbundance_Nitrogen:  8.3562563e-4
+  InitAbundance_Oxygen:    5.4926244e-3
+  InitAbundance_Neon:      1.4144605e-3
+  InitAbundance_Magnesium: 5.907064e-4
+  InitAbundance_Silicon:   6.825874e-4
+  InitAbundance_Iron:      1.1032152e-3
+  CalciumOverSilicon:      0.0941736
+  SulphurOverSilicon:      0.6054160
+
+EagleCooling:
+  filename:                /cosma5/data/Eagle/BG_Tables/CoolingTables/
+  reionisation_redshift:   8.898
+
+#Cosmology:
+#  Omega_m:        0.307
+#  Omega_lambda:   0.693
+#  Omega_b:        0.0455
+#  h:              0.6777
+#  #a_begin:        0.0078125
+#  a_begin:        1.0
+#  #a_begin:        0.5
+#  a_end:          1.0
+
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+

examples/EAGLE_6/run.sh

@@ -7,5 +7,5 @@ then
     ./getIC.sh
 fi
 
-../swift -c -s -G -S -t 16 eagle_6.yml 2>&1 | tee output.log
+../swift -c -C -s -t 16 eagle_6.yml 2>&1 | tee output.log
 

examples/main.c

@@ -51,6 +51,12 @@
 /* Global profiler. */
 struct profiler prof;
 
+/* number of calls to eagle cooling rate */
+int n_eagle_cooling_rate_calls_1 = 0;
+int n_eagle_cooling_rate_calls_2 = 0;
+int n_eagle_cooling_rate_calls_3 = 0;
+int n_eagle_cooling_rate_calls_4 = 0;
+
 /**
  * @brief Help messages for the command line parameters.
  */
@@ -125,7 +131,7 @@ int main(int argc, char *argv[]) {
 
   struct clocks_time tic, toc;
   struct engine e;
-
+  
   /* Structs used by the engine. Declare now to make sure these are always in
    * scope.  */
   struct chemistry_global_data chemistry;

src/cooling/EAGLE/cooling_struct.h

@@ -120,6 +120,8 @@ struct cooling_function_data {
   float calcium_over_silicon_ratio;
   float sulphur_over_silicon_ratio;
 
+  double delta_u;
+
 };
 
 /**

src/cooling/EAGLE/eagle_cool_tables.h

@@ -610,7 +610,7 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
       for (i = 0; i < cooling->N_nH; i++){
         for (j = 0; j < cooling->N_Temp; j++){
           table_index = row_major_index_2d(j,i,cooling->N_Temp,cooling->N_nH);
-          cooling_index = row_major_index_4d(z_index,specs,j,i,cooling->N_Redshifts,cooling->N_Elements,cooling->N_Temp,cooling->N_nH); 
+          cooling_index = row_major_index_4d(z_index,specs,i,j,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); 
           //cooling_index = row_major_index_3d(specs,j,i,cooling->N_Elements,cooling->N_Temp,cooling->N_nH); 
           cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
         }
@@ -640,7 +640,7 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
       for (j = 0; j < cooling->N_Temp; j++){
         for (k = 0; k < cooling->N_nH; k++) {
           table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
-          cooling_index = row_major_index_4d(z_index,i,j,k,cooling->N_Redshifts,cooling->N_He,cooling->N_Temp,cooling->N_nH); 
+          cooling_index = row_major_index_4d(z_index,k,i,j,cooling->N_Redshifts,cooling->N_nH,cooling->N_He,cooling->N_Temp); 
           //cooling_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH); 
           cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
           cooling_table.H_plus_He_electron_abundance[cooling_index] =
@@ -659,7 +659,7 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
     for (i = 0; i < cooling->N_Temp; i++){
       for (j = 0; j < cooling->N_nH; j++){
         table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
-        cooling_index = row_major_index_3d(z_index,i,j,cooling->N_Redshifts,cooling->N_Temp,cooling->N_nH); 
+        cooling_index = row_major_index_3d(z_index,j,i,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); 
         //cooling_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH); 
         cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
       }

tests/plots.py

@@ -3,6 +3,8 @@ import matplotlib.pyplot as plt
 elements = 11
 temperature = []
 cooling_rate = [[] for i in range(elements+1)]
+u = []
+fu = []
 length = 0
 file_in = open('cooling_output.dat', 'r')
 for line in file_in:
@@ -19,22 +21,38 @@ for elem in range(elements):
         	cooling_rate[elem+1].append(-float(data[0]))
 	file_in.close()
 
+file_in = open('newton_output.dat', 'r')
+for line in file_in:
+	data = line.split()
+	u.append(float(data[0]))
+	fu.append(float(data[1]))
 
-p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
-p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
-p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
-p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
-p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
-p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
-p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
-p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
-p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
-p10, = plt.loglog(temperature, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
-plt.ylim([1e-24,1e-21])
-plt.xlim([1e4,1e8])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
-#plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
+file_in.close()
+
+p0, = plt.plot(u,fu, color = 'k')
+#p1, = plt.plot(u,[0 for x in u], color = 'r')
+#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
+plt.xlim([1e13,2.0e14])
+plt.ylim([0e13,2e14])
+plt.xlabel('u')
+plt.ylabel('f(u)')
 plt.show()
+
+#p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+#p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
+#p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
+#p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+#p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
+#p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
+#p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
+#p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
+#p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+#p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
+#p10, = plt.loglog(temperature, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
+#plt.ylim([1e-24,1e-21])
+#plt.xlim([1e4,1e8])
+#plt.xlabel('Temperature (K)')
+#plt.ylabel('Cooling rate (eagle units...)')
+#plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
+##plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
+#plt.show()

tests/testCooling.c

@@ -49,9 +49,10 @@ int main() {
   units_init(&us, params, "InternalUnitSystem");
   phys_const_init(&us, params, &internal_const);
 
-  double hydrogen_number_density_cgs = 1e-4;
-  double u, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs;
+  double hydrogen_number_density_cgs = 1e-1;
+  double u,u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
   int n_t_i = 2000;
+  dt_cgs = 4.0e-8*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
 
   char *output_filename = "cooling_output.dat";
   FILE *output_file = fopen(output_filename, "w");
@@ -67,6 +68,13 @@ int main() {
       printf("Error opening file!\n");
       exit(1);
   }
+  char *output_filename3 = "newton_output.dat";
+  FILE *output_file3 = fopen(output_filename3, "w");
+  if (output_file3 == NULL)
+  {
+      printf("Error opening file!\n");
+      exit(1);
+  }
 
   gamma = 5.0/3.0;
 
@@ -77,12 +85,9 @@ int main() {
   cosmology_init(params, &us, &internal_const, &cosmo);
   cosmology_print(&cosmo);
 
-  //float scale_factor = 1.0/(1.0 + cosmo.z);
-  //float scaled_hydrogen_number_density_cgs = hydrogen_number_density_cgs/pow(scale_factor,3);
     
   u = 1.0*pow(10.0,11)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
   pressure = u*p.rho*(gamma -1.0);
-  //hydrogen_number_density = scaled_hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
   hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
   p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
   p.entropy = pressure/(pow(p.rho,gamma));
@@ -91,22 +96,43 @@ int main() {
   cooling_init(params, &us, &internal_const, &cooling);
   cooling_print(&cooling);
 
+  // construct 1d table of cooling rates wrt temperature
+  float H_plus_He_heat_table[176];                      // WARNING sort out how it is declared/allocated
+  float H_plus_He_electron_abundance_table[176];        // WARNING sort out how it is declared/allocated
+  float element_cooling_table[9*176];                   // WARNING sort out how it is declared/allocated
+  float element_electron_abundance_table[176];          // WARNING sort out how it is declared/allocated
+  construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
+  construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_electron_abundance,H_plus_He_electron_abundance_table);
+  construct_1d_table_from_4d_elements(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.metal_heating,element_cooling_table);
+  construct_1d_table_from_3d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.electron_abundance,element_electron_abundance_table);
+
+  float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
+  float inn_h = chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
+  ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+  u_ini_cgs = pow(10.0,14);
 
   for(int t_i = 0; t_i < n_t_i; t_i++){
     
-    u = 1.0*pow(10.0,11 + t_i*6.0/n_t_i)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
+    u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
+    u = u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
     pressure = u*p.rho*(gamma -1.0);
     hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
     p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
     p.entropy = pressure/(pow(p.rho,gamma));
 
-    cooling_du_dt = eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
-    temperature_cgs = eagle_convert_u_to_temp(&p,&cooling,&cosmo,&internal_const);
+    //cooling_du_dt = eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
+    cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,&p,&cooling,&cosmo,&internal_const);
+    temperature_cgs = eagle_convert_u_to_temp(u_cgs,&p,&cooling,&cosmo,&internal_const);
     fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
     fprintf(output_file2,"%.5e %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)), temperature_cgs);
+
+    newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
+    fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
+    //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
   }
   fclose(output_file);
   fclose(output_file2);
+  fclose(output_file3);
 
   free(params);