Erick/production 1.0

Scripts/initial_conditions/st9_empty_particles_multi_energy.py

@@ -16,12 +16,22 @@
 nphi_equator = 16 # number of direction in equator
 NF = 3 # number of flavors
 
+'''
+--------- 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]
+'''
+
+# 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

Scripts/symbolic_hermitians/generate_code.py

@@ -526,4 +526,23 @@ def sgn(t,var):
         code.append(dNdt.code())
 
     code = [line for sublist in code for line in sublist]
-    write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_dfdt_fill_zeros"))
+    write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_dfdt_fill_zeros"))
+
+    #========================#
+    # Evolve.cpp_compute_trace #
+    #========================#
+
+    # List that will store the code for setting the derivative of the matrices N and Nbar to zero.
+    code = []
+
+    # Looping over neutrinos(tail: no tail) and antineutrinos(tail: bar)
+    for t in tails:
+
+        # Define n and nbar matrix
+        N = HermitianMatrix(args.N, "fab(i\,j\,k\,GIdx::N{}{}_{}"+t+")")
+        # Assign the trace of n and nbar to a variable
+        trace_N = N.trace()
+        code.append(["trace_n"+t+" = " + sympy.cxxcode(sympy.simplify(trace_N)) + ";"])
+
+    code = [line for sublist in code for line in sublist]
+    write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_compute_trace"))

Scripts/tests/coll_inst_test.py

@@ -43,8 +43,8 @@
             t[i]            = np.array(hf['t(s)'][:][0])
 
     # Fit the exponential function ( y = a e ^ ( b x ) ) to the data
-    l1 = 10 # initial item for fit
-    l2 = 40 # last item for fit
+    l1 = 30 # initial item for fit
+    l2 = 70 # last item for fit
     coefficients = np.polyfit(t[l1:l2], np.log(N_avg_mag[:,0,1][l1:l2]), 1)
     coefficients_bar = np.polyfit(t[l1:l2], np.log(Nbar_avg_mag[:,0,1][l1:l2]), 1)
     a = np.exp(coefficients[1])

Scripts/tests/fermi_dirac_test.py

@@ -19,11 +19,11 @@
 directories = sorted(directories, key=lambda x: int(x.split("plt")[1].split(".")[0]))
 
 # Energy bin centers extracted from NuLib table
-energies_center_Mev = np.array([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
+energies_center_Mev = np.array([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 = np.array([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])
+energies_bottom_Mev = np.array([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 = np.array([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])
+energies_top_Mev = np.array([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 * eV # Energy in ergs

Source/Evolve.H

@@ -26,7 +26,40 @@ namespace GIdx
     void Initialize();
 }
 
-amrex::Real compute_dt(const amrex::Geometry& geom, const MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, const TestParams* parms);
+/**
+ * @brief Computes the maximum Inverse Mean Free Path (IMFP) for absorption and emission processes.
+ *
+ * This function calculates the maximum IMFP for absorption processes based on the specified method in the parameters.
+ * It supports two methods:
+ * 1. Using IMFP values from the input file.
+ * 2. Using a NuLib table to interpolate IMFP values.
+ *
+ * @param geom The geometry of the computational domain.
+ * @param state The MultiFab containing the state variables.
+ * @param parms Pointer to the TestParams structure containing the parameters for the computation.
+ * @return A MultiFab containing the maximum IMFP values.
+ */
+MultiFab compute_max_IMFP(const Geometry& geom, const MultiFab& state, const TestParams* parms);
+
+/**
+ * @brief Computes the time step size for the simulation based on various CFL factors and conditions.
+ *
+ * This function calculates the time step size considering translation, flavor, and collision CFL factors.
+ * It ensures that the time step is limited by the smallest of these factors to maintain stability.
+ *
+ * @param geom The geometry of the simulation domain.
+ * @param state The state MultiFab containing the simulation data.
+ * @param neutrinos The container for flavored neutrinos.
+ * @param parms Pointer to the structure containing simulation parameters.
+ * @param maximum_IMFP_abs The MultiFab containing the maximum inverse mean free path for absorption.
+ *
+ * @return The computed time step size.
+ *
+ * @note At least one of cfl_factor, flavor_cfl_factor, or collision_cfl_factor must be greater than 0.0.
+ * @note The function handles periodic boundary conditions and black hole regions if specified in the parameters.
+ * @note The function performs a reduction operation to find the minimum time step across all cells and MPI ranks.
+ */
+amrex::Real compute_dt(const amrex::Geometry& geom, const MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, const TestParams* parms, const MultiFab& maximum_IMFP_abs);
 
 void deposit_to_mesh(const FlavoredNeutrinoContainer& neutrinos, amrex::MultiFab& state, const amrex::Geometry& geom);