Erick/production 1.0

Scripts/collisions/compute_dE_coll_inst_test.py

@@ -18,16 +18,19 @@
 hc = h * c # erg cm
 
 # Simulation parameters 
-V = 1 # Volume of a cell ( ccm ) 
+V = (1.0e3)**3 # Volume of a cell ( ccm ) 
 Ndir = 92 # Number of momentum beams isotropically distributed per cell 
 E = 20.0 # Neutrinos and antineutrinos energy bin center ( Mev )
 T = 7.0 # Background matter temperature ( Mev )
 
-N_eq_electron_neutrino = 3.260869565e+31 # Number of electron neutrinos at equilibrium
+N_eq_electron_neutrino_density = 3.0e33 # Density of electron neutrinos at equilibrium  (ccm)
+N_eq_electron_neutrino = N_eq_electron_neutrino_density * V / Ndir # Number of electron neutrinos at equilibrium per beam
+print(f'Number of electron neutrinos at equilibrium per beam = {N_eq_electron_neutrino}')
 u_electron_neutrino = 20.0 # Electron neutrino chemical potential ( Mev )
 
 # Fermi-dirac distribution factor for electron neutrinos
 f_eq_electron_neutrinos = 1 / ( 1 + np.exp( ( E - u_electron_neutrino ) / T ) ) # adimentional
+print(f'Fermi-Dirac distribution factor for electron neutrinos = {f_eq_electron_neutrinos}')
 
 # We know : 
 #   dE^3      = 3 *          Neq           * ( hc )^ 3 / ( dV *        dOmega       *        feq              )
@@ -47,9 +50,13 @@
         dE=np.real(complex_deltaE)
 
 # Electron neutrino flavor
-N_eq_electron_antineutrino = 2.717391304e+31 # Number of electron antineutrinos at equilibrium
+N_eq_electron_antineutrino_density = 2.5e33 # Density of electron antineutrinos at equilibrium (ccm)
+N_eq_electron_antineutrino = N_eq_electron_antineutrino_density * V / Ndir # Number of electron antineutrinos at equilibrium per beam
+print(f'Number of electron antineutrinos at equilibrium per beam = {N_eq_electron_antineutrino}')
+Vphase = V * ( 4 * np.pi / Ndir ) * delta_E_cubic  / 3
+print(f'Phase space volume in ccm MeV^3 = {Vphase}')
 
 # Computing electron antineutrino chemical potential
 f_eq_electron_antineutrino = 3 * N_eq_electron_antineutrino * ( hc )**3 / ( V * ( 4 * np.pi / Ndir ) * delta_E_cubic )
 u_electron_antineutrino = E - T * np.log( 1 / f_eq_electron_antineutrino - 1 )
-print(f'Electron neutrino chemical potential in MeV = {u_electron_antineutrino}')
+print(f'Electron antineutrino chemical potential in MeV = {u_electron_antineutrino}')

Scripts/collisions/single_particle_plot.py

@@ -0,0 +1,100 @@
+'''
+Created by Erick Urquilla. University of Tennessee Knoxville, USA.
+This script is used to plot the single particle polarization vector in the x, y, and z directions.
+This script works exclusively for the two flavor case.
+The script reads the plt files generated with the writeparticleinfohdf5.py script.
+'''
+
+import numpy as np
+import glob
+import h5py
+import matplotlib.pyplot as plt
+
+file_names = np.array(glob.glob('plt*.h5'))
+file_names = [file_name for file_name in file_names if '_reduced_data' not in file_name]
+file_names.sort(key=lambda x: int(x.split('plt')[1].split('.h5')[0]))
+
+particle_momentum = [0,0,0,0]
+with h5py.File(file_names[0], 'r') as hdf:
+    data = {}
+    for key in hdf.keys():
+        data[key] = np.array(hdf.get(key))
+    particle_momentum[0] = data['pupx'][0]
+    particle_momentum[1] = data['pupy'][0]
+    particle_momentum[2] = data['pupz'][0]
+    particle_momentum[3] = data['pupt'][0]
+
+time = np.zeros(len(file_names))
+N00_Re = np.zeros(len(file_names))
+N00_Rebar = np.zeros(len(file_names))
+N01_Im = np.zeros(len(file_names))
+N01_Imbar = np.zeros(len(file_names))
+N01_Re = np.zeros(len(file_names))
+N01_Rebar = np.zeros(len(file_names))
+N11_Re = np.zeros(len(file_names))
+N11_Rebar = np.zeros(len(file_names))
+
+for i, file in enumerate(file_names):
+    with h5py.File(file, 'r') as hdf:
+
+        # ['N00_Re', 'N00_Rebar', 'N01_Im', 'N01_Imbar', 'N01_Re', 'N01_Rebar', 'N11_Re', 'N11_Rebar', 'TrHN', 'Vphase', 'pos_x', 'pos_y', 'pos_z', 'pupt', 'pupx', 'pupy', 'pupz', 'time', 'x', 'y', 'z']
+
+        data = {}
+        for key in hdf.keys():
+            data[key] = np.array(hdf.get(key))
+
+        for j in range(len(data['pupx'])):
+            if (data['pupx'][j] == particle_momentum[0] and
+            data['pupy'][j] == particle_momentum[1] and
+            data['pupz'][j] == particle_momentum[2] and
+            data['pupt'][j] == particle_momentum[3]):
+                time[i] = data['time'][j]
+                N00_Re[i] = data['N00_Re'][j]
+                N00_Rebar[i] = data['N00_Rebar'][j]
+                N01_Im[i] = data['N01_Im'][j]
+                N01_Imbar[i] = data['N01_Imbar'][j]
+                N01_Re[i] = data['N01_Re'][j]
+                N01_Rebar[i] = data['N01_Rebar'][j]
+                N11_Re[i] = data['N11_Re'][j]
+                N11_Rebar[i] = data['N11_Rebar'][j]
+                break
+
+P_x = 0.5 * N01_Re
+P_y = 0.5 * N01_Im
+P_z = 0.5 * (N00_Re - N11_Re)
+
+# Plot P_x vs P_y
+plt.plot(P_x, P_y, color='gray', alpha=0.5)
+sc = plt.scatter(P_x, P_y, c=time, cmap='viridis')
+plt.colorbar(sc, label=r'$t \, (s)$')
+plt.xlabel(r'$P_x$')
+plt.ylabel(r'$P_y$')
+plt.savefig('particle_momentum_plot_Px_Py.pdf')
+plt.close()
+
+# Plot log(P_x) vs log(P_y)
+plt.plot(np.log10(np.abs(P_x)), np.log10(np.abs(P_y)), color='gray', alpha=0.5)
+sc = plt.scatter(np.log10(np.abs(P_x)), np.log10(np.abs(P_y)), c=time, cmap='viridis')
+plt.colorbar(sc, label=r'$t \, (s)$')
+plt.xlabel(r'$\log_{10}(|P_x|)$')
+plt.ylabel(r'$\log_{10}(|P_y|)$')
+plt.savefig('particle_momentum_plot_log_Px_Py.pdf')
+plt.close()
+
+# Plot P_x vs P_z
+plt.plot(P_x, P_z, color='gray', alpha=0.5)
+sc = plt.scatter(P_x, P_z, c=time, cmap='viridis')
+plt.colorbar(sc, label=r'$t \, (s)$')
+plt.xlabel(r'$P_x$')
+plt.ylabel(r'$P_z$')
+plt.savefig('particle_momentum_plot_Px_Pz.pdf')
+plt.close()
+
+# Plot P_y vs P_z
+plt.plot(P_y, P_z, color='gray', alpha=0.5)
+sc = plt.scatter(P_y, P_z, c=time, cmap='viridis')
+plt.colorbar(sc, label=r'$t \, (s)$')
+plt.xlabel(r'$P_y$')
+plt.ylabel(r'$P_z$')
+plt.savefig('particle_momentum_plot_Py_Pz.pdf')
+plt.close()

Scripts/data_reduction/convertToHDF5.py

@@ -49,14 +49,19 @@ def convert_to_HDF5(sim_directory, DELETE_ALL_BUT_LAST_RESTART=False):
         if d==fluid_directories[0]:
             NF = eds.get_num_flavors()
             allData = h5py.File(sim_directory+"/allData.h5","w")
+            allData["dx(cm)"] = eds.dx
+            allData["dy(cm)"] = eds.dy
             allData["dz(cm)"] = eds.dz
+            allData["Nx"] = eds.Nx
+            allData["Ny"] = eds.Ny
+            allData["Nz"] = eds.Nz
             allData.create_dataset("t(s)", data=np.zeros(0), maxshape=(None,), dtype=datatype)
             allData.create_dataset("it", data=np.zeros(0), maxshape=(None,), dtype=int)
 
             # create fluid data sets
-            maxshape = (None, eds.Nz)
-            chunkshape = (1, eds.Nz)
-            zeros = np.zeros((0,eds.Nz))
+            maxshape = (None,eds.Nx,eds.Ny,eds.Nz)
+            chunkshape = (1,eds.Nx,eds.Ny,eds.Nz)
+            zeros = np.zeros((0,eds.Nx,eds.Ny,eds.Nz))
             varlist = []
             for v in fluid_vars:
                 for f1 in range(NF):

Scripts/initial_conditions/st9_empty_particles_multi_energy.py

@@ -16,12 +16,28 @@
 nphi_equator = 16 # number of direction in equator
 NF = 3 # number of flavors
 
+'''
+# Energy bins from NuLib table (No cutted)
+--------- NuLib table energy bins ---------
 # Energy bin centers extracted from NuLib table
 energies_center_Mev = [1, 3, 5.23824, 8.00974, 11.4415, 15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542] # Energy in Mev
 # Energy bin bottom extracted from NuLib table
 energies_bottom_Mev = [0, 2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250]
 # Energy bin top extracted from NuLib table
 energies_top_Mev = [2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085]
+'''
+
+# The following energy bins are an extraction of the NuLib table
+# I have delete the first five energy bin to make the fermi-dirac test facter
+# If I simulation want to run with all energy bins the user must change the energy bins manually
+
+# Energy bins from NuLib table (Cutted)
+# Energy bin centers extracted from NuLib table
+energies_center_Mev = [15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542] # Energy in Mev
+# Energy bin bottom extracted from NuLib table
+energies_bottom_Mev = [13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250]
+# Energy bin top extracted from NuLib table
+energies_top_Mev = [18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085]
 
 # Energies in ergs
 energies_center_erg = np.array(energies_center_Mev) * 1e6*amrex.eV # Energy in ergs