(DRAFT) feat: introduce 3-velocity grid variable to support relativistic transport

Jenkinsfile

@@ -44,7 +44,7 @@ pipeline {
 
 		stage('Fast Flavor'){ steps{
 				dir('Exec'){
-				sh 'python ../Scripts/initial_conditions/st2_2beam_fast_flavor.py'
+					sh 'python ../Scripts/initial_conditions/st2_2beam_fast_flavor.py'
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor'
 					sh 'python ../Scripts/tests/fast_flavor_test.py'
 					sh 'rm -rf plt*'
@@ -54,7 +54,7 @@ pipeline {
 
 		stage('Fast Flavor k'){ steps{
 				dir('Exec'){
-				sh 'python ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py'
+					sh 'python ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py'
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor_nonzerok'
 					sh 'python ../Scripts/tests/fast_flavor_k_test.py'
 					sh 'rm -rf plt*'
@@ -80,26 +80,13 @@ pipeline {
 			}
 		}
 
-		stage('Fiducial 3F CPU HDF5'){ steps{
-				dir('Exec'){
-					sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5 GNUmakefile'
-					sh 'make realclean; make generate; make -j'
-					sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py'
-					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial'
-					/*sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'*/
-					sh 'rm -rf plt*'
-				}
-			}
-		}
-
 		stage('Collisions flavor instability'){ steps{
 				dir('Exec'){
-					sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile'
-					sh 'make realclean; make generate NUM_FLAVORS=2; make -j NUM_FLAVORS=2'
 					sh 'python ../Scripts/initial_conditions/st8_coll_inst_test.py'
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_collisional_instability_test'
 					sh 'python ../Scripts/data_reduction/reduce_data.py'
 					sh 'python ../Scripts/tests/coll_inst_test.py'
+					archiveArtifacts artifacts: '*.pdf'
 					sh 'rm -rf plt* *pdf'
 				}
 			}
@@ -109,26 +96,40 @@ pipeline {
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_bc_periodic_init'
 					sh 'python ../Scripts/collisions/writeparticleinfohdf5.py'
 					sh 'python ../Scripts/tests/bc_empty_init_test.py'
+					archiveArtifacts artifacts: '*.pdf'
 					sh 'rm -rf plt* *pdf'
 				}
 			}
 		}
+
+		stage('Fiducial 3F CPU HDF5'){ steps{
+				dir('Exec'){
+					sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5 GNUmakefile'
+					sh 'make realclean; make generate; make -j'
+					sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py'
+					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial'
+					/*sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'*/
+					sh 'rm -rf plt*'
+				}
+			}
+		}
+
 		stage('Collisions to equilibrium'){ steps{
 				dir('Exec'){
 					sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile'
 					sh 'make realclean; make generate NUM_FLAVORS=3; make -j NUM_FLAVORS=3'
 					sh 'python ../Scripts/initial_conditions/st7_empty_particles.py'
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_coll_equi_test'
 					sh 'python ../Scripts/tests/coll_equi_test.py'
-					sh 'rm -rf plt* *pdf'
+					sh 'rm -rf plt*'
 				}
 			}
 		}
 
 		stage('Fermi-Dirac test'){ steps{
 				dir('Exec'){
 					sh 'python ../Scripts/initial_conditions/st9_empty_particles_multi_energy.py'
-					sh 'python ../Scripts/collisions/nsm_constant_background_rho_Ye_T__writer.py'
+					sh 'python ../Scripts/collisions/nsm_constant_background_rho_Ye_T_writer.py'
 					sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fermi_dirac_test'
 					sh 'python ../Scripts/collisions/writeparticleinfohdf5.py'
 					sh 'python ../Scripts/tests/fermi_dirac_test.py'

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/convert_cellh5_inidat.py

@@ -0,0 +1,79 @@
+'''
+Created by Erick Urquilla, Department of Physics and Astronomy, University of Tennessee, Knoxville.
+This script convert the cell hdf5 files into a particle_input.dat file that can be used as an initial condition
+for the Emu code.
+The cell hdf5 files are generated by the script write_particles_per_cell.py
+The particle data is stored in hdf5 files with the following labels: cell_i_j_k.h5
+where i, j, k are the cell indexes in the x, y, z directions respectively.
+The ddf5 files has to be in the directory this script is run.
+'''
+
+import glob
+import h5py
+import numpy as np
+import sys
+import os
+importpath = os.path.dirname(os.path.realpath(__file__))
+sys.path.append(importpath)
+sys.path.append(importpath+"/../initial_conditions")
+from initial_condition_tools import uniform_sphere, write_particles
+import amrex_plot_tools as amrex
+
+def read_hdf5_file(file_name):
+    """
+    Reads an HDF5 file and returns its datasets as a dictionary.
+
+    Parameters:
+        file_name (str): Path to the HDF5 file.
+
+    Returns:
+        dict: A dictionary where keys are dataset names and values are NumPy arrays.
+    """
+    with h5py.File(file_name, 'r') as file:
+        return {dataset: np.array(file[dataset]) for dataset in file.keys()}
+
+# Find all files matching the pattern
+file_pattern = 'cell_*_*_*.h5'
+files = glob.glob(file_pattern)
+
+data = read_hdf5_file(files[0])
+
+labels_nf2=['pos_x','pos_y','pos_z', 'time', 'x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N11_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N11_Rebar', 'TrHN', 'Vphase']
+labels_nf3=['pos_x','pos_y','pos_z','time','x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N02_Re', 'N02_Im', 'N11_Re', 'N12_Re', 'N12_Im' ,'N22_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N02_Rebar', 'N02_Imbar', 'N11_Rebar', 'N12_Rebar' ,'N12_Imbar', 'N22_Rebar', 'TrHN', 'Vphase']
+
+NF = 0
+if len(data.keys()) == len(labels_nf2):
+    labels = labels_nf2
+    NF = 2
+elif len(data.keys()) == len(labels_nf3):
+    labels = labels_nf3
+    NF = 3
+else:
+    print(f'len(data.keys()) = {len(data.keys())} does not match any of the known labels')
+    exit(1)
+
+# Get variable keys
+rkey, ikey = amrex.get_particle_keys(NF, ignore_pos=True)
+
+# Generate the number of variables that describe each particle
+n_variables = len(rkey)
+
+# Generate the number of particles
+n_particles = len(data[labels[0]])
+
+# Initialize a NumPy array to store all particles
+particles = np.zeros((n_particles, n_variables))
+
+# Volume of the simulation box
+volume = 1e5**3 # ccm
+
+for label in rkey:
+
+    print(f'Writing {label} to particles')
+    particles[:, rkey[label]] = data[label]
+
+    if (label.startswith('N') | label.startswith('Vphase')):
+        particles[:, rkey[label]] /= volume
+
+# Write particles initial condition file
+write_particles(np.array(particles), NF, "particle_input.dat")

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()