(DRAFT) Relativistic transport dev

Scripts/symbolic_hermitians/generate_code.py

@@ -408,8 +408,27 @@ def sgn(t,var):
 
         # matter potential
         if("V00" in Vlist[icomp]):
-            line = line + " + " + rhoye
+
+            # We will introduce the relativistic correction in a separate variable, the frame-invariant correction [ -p^a u_a / (p^t) ]
+            lorentz_line = "double lorentz_factor = 1.0 / sqrt(1.0 - "
+            for i in range(len(direction)):
+                lorentz_line = lorentz_line + "std::pow(" + string_interp + "vup"+direction[i]+"), 2)"
+                if direction[i] != direction[-1]:
+                    lorentz_line = lorentz_line + " + "
+            lorentz_line = lorentz_line + ");"
+            code.append(lorentz_line)
+            velocity_line = "double relativistic_correction = (-1/p.rdata(PIdx::pupt)) * lorentz_factor * ("
 
+            for i in range(len(direction)):
+                # using (-1)vupt = vdownt for now; will need to implement contraction with the metric once the metric is stored 
+                # use sympy to create a matrix for the metric tensor and use symbolic matrix operations built in
+                velocity_line = velocity_line + "p.rdata(PIdx::pup"+direction[i]+") * " + string_interp + "vup"+direction[i]+") + "
+            # Again, using (-1)*vup for now; will need to implement contraction with the metric once metric is stored
+            velocity_line = velocity_line + " p.rdata(PIdx::pupt) * (-1)"
+            velocity_line = velocity_line + ");"
+            code.append(velocity_line)
+            # We can easily comment out the relativistic correction below, for testing purposes 
+            line = line + " + " +  rhoye  + " * relativistic_correction "
         line = line + ";"
         code.append(line)
         code.append("")

Scripts/tests/msw_test.py

@@ -7,6 +7,7 @@
 importpath = os.path.dirname(os.path.realpath(__file__))+"/../data_reduction/"
 sys.path.append(importpath)
 import amrex_plot_tools as amrex
+import matplotlib.pyplot as plt 
 
 parser = argparse.ArgumentParser()
 parser.add_argument("-na", "--no_assert", action="store_true", help="If --no_assert is supplied, do not raise assertion errors if the test error > tolerance.")
@@ -111,3 +112,10 @@ def myassert(condition):
     print("conservation_errorbar:", conservation_errorbar)
     myassert(conservation_errorbar < tolerance)
 
+plt.plot(t, Nee_analytic, linestyle='-', label='Analytic')
+plt.plot(t, Nee, linestyle='--', label="EMU")
+plt.title("Difference in MSW Analytic Solution and Relativistic EMU Solution")
+plt.xlabel("Time (s)")
+plt.ylabel("N_ee (cm^-3)")
+plt.legend()
+plt.savefig("msw.pdf")

Source/Evolve.H

@@ -18,7 +18,10 @@ namespace GIdx
     enum {
         rho, T, Ye, // g/ccm, MeV, unitless
         #include "generated_files/Evolve.H_fill"
-        ncomp
+        vupx,
+	vupy,
+	vupz,
+	ncomp
     };
 
     extern amrex::Vector<std::string> names;

Source/Parameters.H

@@ -25,6 +25,9 @@ struct TestParams : public amrex::Gpu::Managed
     Real rho_in; // Density in g/ccm
     Real Ye_in; // Electron fraction (dimensionless quantity)
     Real kT_in; // Temperature in erg
+    Real vupx_in; // X-component of velocity, cm/s
+    Real vupy_in; // Y-component, cm/s
+    Real vupz_in ; // Z-component, cm/s
     Real cfl_factor, flavor_cfl_factor, collision_cfl_factor;
     Real max_adaptive_speedup; 
     bool do_restart; // Flag to restart from a previous simulation
@@ -143,6 +146,10 @@ struct TestParams : public amrex::Gpu::Managed
             pp.get("Ye", Ye_in); // Electron fraction (dimensionless quantity)
             pp.get("T_MeV", kT_in); // Temperature in erg
             kT_in *= 1e6 * CGSUnitsConst::eV; // convert to cgs : erg
+            pp.get("vupx", vupx_in); // Velocity cm s^-1
+            pp.get("vupy", vupy_in); // Velocity cm s^-1
+            pp.get("vupz", vupz_in); // Velocity cm s^-1
+
         } else if(read_rho_T_Ye_from_table==1){ 
             // Read background matter from a table
             pp.get("background_rho_Ye_T_table_name", background_rho_Ye_T_table_name);

Source/main.cpp

@@ -117,10 +117,16 @@ void evolve_flavor(const TestParams* parms)
     //Else set rho, T and Ye to constant value throughout the grid using values from parameter file.
     if (parms->read_rho_T_Ye_from_table){
         set_rho_T_Ye(state, geom, parms);
+	state.setVal(parms->vupx_in, GIdx::vupx,1); // cm/s 
+	state.setVal(parms->vupy_in, GIdx::vupy,1); // cm/s
+	state.setVal(parms->vupz_in, GIdx::vupz,1); // cm/s
     } else {      
         state.setVal(parms->rho_in,GIdx::rho,1); // g/ccm
-        state.setVal(parms->Ye_in,GIdx::Ye,1);
-        state.setVal(parms->kT_in,GIdx::T,1); // erg
+        state.setVal(parms->Ye_in,GIdx::Ye,1); 
+        state.setVal(parms->kT_in,GIdx::T,1); // erg/s
+	state.setVal(parms->vupx_in, GIdx::vupx,1); // cm/s
+	state.setVal(parms->vupy_in, GIdx::vupy,1); // cm/s
+	state.setVal(parms->vupz_in, GIdx::vupz,1); // cm/s
     }
 
     state.FillBoundary(geom.periodicity());

sample_inputs/inputs_msw_test

@@ -39,6 +39,10 @@ amrex.fpe_trap_invalid=1
 rho_g_ccm = 1346160134.709288
 T_MeV = 10
 Ye = 1
+vupx = 0
+vupy = 0
+vupz = 0
+
 
 # Write plotfiles
 write_plot_every = 1
@@ -102,4 +106,4 @@ do_blackhole = 0
 # Boundary conditions
 #################
 boundary_condition_type = 0
-do_periodic_empty_bc = 0
+do_periodic_empty_bc = 0