Development (#66)

* Added descriptive amrex::Error messages in generate_code.py (#53) * Added descriptive amrex::Error messages in generate_code.py Swtiched amrex::Abort to amrex::Error with descriptive error message * added more to amrex::Error messages * added grid variable to use in timestep calculation * add ability to do scalar math with HermitianMatrix-es * added 3 flavor functionality to sample_inputs test simulations. Note: mass3_eV = 0 for these to run. (#56) * update readme for code generation and fix commutator language (#55) * rewrite how the magnitude of a complex number is written so sympy simplifies the result without extraneous Is in real numbers * use amrex reduction functions to evaluate and reduce Hamiltonian flavor vector length in place * reduce the flavor cfl factor in the input files to account for the fact that the timestep is more liberal now * separate the SU vector magnitude function into one that returns the magnitude and one that returns the magnitude squared. This will allow us to add the flux components without repeatedly square rooting and squaring. * added flux term to max hamiltonian for timestep calculation * make the timestep account for the vacuum potential as well. Compute the minimum particle energy in initializing step to avoid doing a reduction over particles at every timestep under the assumption that particle energies do not change * fixed unit issue with max vacuum potential * override copyParticles function to also copy the new Vvac_max member variable * fixed broken implementation of conjugate * fixed error in generate code - the length of the Hamiltonian vector needs antineutrinos to be conjugated * pushed flavor cfl factors back up as high as they can go without runing into the maxError limit * set the reduction operator to return a second quantity. Will be used to limit timestep based on total lepton density * account for vacuum potential. Also cover case where total potential is 0 because of equal numbers of neutrinos and antineutrinos by introducing a parameter (max_adaptive_speedup). When set to 0 it behaves exactly as before, except for the inclusion of the vacuum potential in the timestep calculation. When set towards infinity, it always allows the "optimized" timestep, which may itself become infinite if terms happen to cancel by chance, but otherwise is a much better estimate of the required timestep * forgot to update the 1d_fiducial inputs file * addressing PR comments. Use static inline variable rather than overriding the ParticleContainer copy operator, fix mistake in conjugation of Hermitian matrices. * Minerbo initial conditions (#59) * added vscode and generated_files directories to the gitignore file * added Minerbo closure initial conditions * added references for the minerbo equations * added a tad more explanatory text * improved st5 code comments * fixed coefficient in front of minerbo closure to make it have an expectation value of 1. Was getting results too small by a factor of 4pi before. Co-authored-by: Sherwood Richers <[email protected]> * Output build info (#63) * Copied Castro function for printing build information to screen * Copy castro function for outputting information into plotfiles, simply deleting lines that don't work in emu Co-authored-by: Sherwood Richers <[email protected]> * im stoopid - fixing timestep logic Co-authored-by: Nicole F <[email protected]> Co-authored-by: Sherwood Richers <[email protected]> Co-authored-by: Don E. Willcox <[email protected]>
AMReX-Astro · Sep 17, 2021 · d25f5bb · d25f5bb
1 parent 79a3fb5
commit d25f5bb
Show file tree

Hide file tree

Showing 20 changed files with 739 additions and 86 deletions.
diff --git a/.gitignore b/.gitignore
@@ -37,4 +37,6 @@ Scripts/symbolic_hermitians/__pycache__
 .cproject
 .project
 .pydevproject
-.settings
+.settings
+.vscode
+Source/generated_files
diff --git a/Make.Emu b/Make.Emu
@@ -24,10 +24,48 @@ include $(Ppack)
 
 DEFINES += -DNUM_FLAVORS=$(NUM_FLAVORS) -DSHAPE_FACTOR_ORDER=$(SHAPE_FACTOR_ORDER)
 
-all: generate $(executable)
+all: generate $(objEXETempDir)/AMReX_buildInfo.o $(executable)
 	@echo SUCCESS
 
 generate:
 	python3 $(EMU_HOME)/Scripts/symbolic_hermitians/generate_code.py $(NUM_FLAVORS) --emu_home $(EMU_HOME)
 
+#------------------------------------------------------------------------------
+# build info (from Castro/Exec/Make.auto_source)
+#------------------------------------------------------------------------------
+CEXE_headers += $(AMREX_HOME)/Tools/C_scripts/AMReX_buildInfo.H
+INCLUDE_LOCATIONS +=  $(AMREX_HOME)/Tools/C_scripts
+
+# we make AMReX_buildInfo.cpp as we make the .o file, so we can delete
+# it immediately.  this way if the build is interrupted, we are
+# guaranteed to remake it
+
+objForExecs += $(objEXETempDir)/AMReX_buildInfo.o
+
+.FORCE:
+.PHONE: .FORCE
+
+# set BUILD_GIT_NAME and BUILD_GIT_DIR if you are building in a
+# git-controlled dir not under Castro/
+EXTRA_BUILD_INFO :=
+ifdef BUILD_GIT_NAME
+   EXTRA_BUILD_INFO := --build_git_name "$(BUILD_GIT_NAME)" \
+                       --build_git_dir "$(BUILD_GIT_DIR)"
+endif
+
+$(objEXETempDir)/AMReX_buildInfo.o: .FORCE
+	echo $(objEXETempDir)
+	$(AMREX_HOME)/Tools/C_scripts/makebuildinfo_C.py \
+          --amrex_home "$(AMREX_HOME)" \
+          --COMP "$(COMP)" --COMP_VERSION "$(COMP_VERSION)" \
+          --CXX_comp_name "$(CXX)" --CXX_flags "$(CXXFLAGS) $(CPPFLAGS) $(includes)" \
+          --F_comp_name "$(F90)" --F_flags "$(F90FLAGS)" \
+          --link_flags "$(LDFLAGS)" --libraries "$(libraries)" \
+          --MODULES "$(MNAMES)" $(EXTRA_BUILD_INFO) \
+          --GIT "$(TOP) $(AMREX_HOME) $(MICROPHYSICS_HOME)"
+	$(SILENT) $(CCACHE) $(CXX) $(CXXFLAGS) $(CPPFLAGS) -c $(CXXEXEFLAGS) AMReX_buildInfo.cpp -o $(objEXETempDir)/AMReX_buildInfo.o
+	$(SILENT) $(RM) AMReX_buildInfo.cpp
+
+
+
 include $(AMREX_HOME)/Tools/GNUMake/Make.rules
diff --git a/Scripts/symbolic_hermitians/HermitianUtils.py b/Scripts/symbolic_hermitians/HermitianUtils.py
@@ -37,18 +37,41 @@ def __init__(self, size, entry_template = "H{}{}_{}"):
         self.construct()
 
     def __mul__(self, other):
-        result = copy.deepcopy(other)
-        result.H = self.H * other.H
+        result = copy.deepcopy(self)
+        if isinstance(other, self.__class__):
+            result.H = self.H * other.H
+        else:
+            for i in range(self.size):
+                for j in range(self.size):
+                    result.H[i,j] *= other
+        return result
+
+    def __truediv__(self, other):
+        result = copy.deepcopy(self)
+        if isinstance(other, self.__class__):
+            raise ArithmeticError
+        else:
+            for i in range(self.size):
+                for j in range(self.size):
+                    result.H[i,j] /= other
         return result
 
     def __add__(self, other):
-        result = copy.deepcopy(other)
-        result.H = self.H + other.H
+        result = copy.deepcopy(self)
+        if isinstance(other, self.__class__):
+            result.H = self.H + other.H
+        else:
+            for i in range(self.size):
+                result.H[i,i] += other
         return result
 
     def __sub__(self, other):
-        result = copy.deepcopy(other)
-        result.H = self.H - other.H
+        result = copy.deepcopy(self)
+        if isinstance(other, self.__class__):
+            result.H = self.H - other.H
+        else:
+            for i in range(self.size):
+                result.H[i,i] -= other
         return result
 
     def construct(self):
@@ -70,16 +93,19 @@ def anticommutator(self, HermitianA, HermitianB):
         return self
 
     def conjugate(self):
-        self.H = Conjugate(self.H)
+        for i in range(self.size):
+            for j in range(self.size):
+                self.H[i,j] = conjugate(self.H[i,j])
         return self
 
     # return the length of the SU(n) vector
-    def SU_vector_magnitude(self):
+    def SU_vector_magnitude2(self):
         # first get the sum of the square of the off-diagonal elements
         mag2 = 0
         for i in range(self.size):
             for j in range(i+1,self.size):
-                mag2 += self.H[i,j]*self.H[j,i]
+                re,im = self.H[i,j].as_real_imag()
+                mag2 += re**2 + im**2
 
         # Now get the contribution from the diagonals
         # See wolfram page for generalization of Gell-Mann matrices
@@ -91,7 +117,10 @@ def SU_vector_magnitude(self):
             basis_coefficient *= sympy.sqrt(2./(l*(l+1.)))
             mag2 += (basis_coefficient/2.)**2
 
-        return sympy.sqrt(mag2)
+        return mag2
+
+    def SU_vector_magnitude(self):
+        return sympy.sqrt(self.SU_vector_magnitude2())
 
     def trace(self):
         result = 0

diff --git a/Scripts/symbolic_hermitians/README.md b/Scripts/symbolic_hermitians/README.md
@@ -5,42 +5,19 @@ expresses a symbolic Hermitian matrix in terms of its
 independent real-valued Real and Imaginary components.
 
 This class also provides functions for calculating an
-anticommutator scaled by I and generating code to implement this.
+commutator scaled by I and generating code to implement this.
 
-## Using HermitianUtils
+For a self-contained demonstration, see the Jupyter notebook
+`Symbolic-Hermitian-Commutator.ipynb`.
 
-The python script `IAntiCommutator.py` gives an example
-of how to use the `HermitianMatrix` class to calculate `C=i*[A,B]`.
+# Using HermitianUtils to generate Emu code
 
-The script generalizes to any desired matrix sizes and
-determines symbol names with runtime arguments.
+We generate all the computational kernels for Emu using the HermitianUtils
+python module in the `generate_code.py` script.
 
-The only required runtime argument is an integer matrix size.
+The generated source files are placed in `Emu/Source/generated_files` and are
+inserted into the Emu source code using compile-time `#include` statements.
 
-To see all optional runtime arguments:
-
-```
-$ python3 IAntiCommutator.py -h
-```
-
-## Generating code with Particle data indexing
-
-To generate sample code with indexing into real Particle data:
-
-```
-$ python3 IAntiCommutator.py 2 -o code.cpp -a "p.rdata(PIdx::H{}{}_{})" -b "p.rdata(PIdx::f{}{}_{})" -c "p.rdata(PIdx::dfdt{}{}_{})"
-```
-
-And the file `code.cpp` will contain:
-
-```
-{
-p.rdata(PIdx::dfdt00_Re) = -2*p.rdata(PIdx::H01_Im)*p.rdata(PIdx::f01_Re) + 2*p.rdata(PIdx::H01_Re)*p.rdata(PIdx::f01_Im);
-
-p.rdata(PIdx::dfdt01_Re) = -p.rdata(PIdx::H00_Re)*p.rdata(PIdx::f01_Im) + p.rdata(PIdx::H01_Im)*p.rdata(PIdx::f00_Re) - p.rdata(PIdx::H01_Im)*p.rdata(PIdx::f11_Re) + p.rdata(PIdx::H11_Re)*p.rdata(PIdx::f01_Im);
-
-p.rdata(PIdx::dfdt01_Im) = p.rdata(PIdx::H00_Re)*p.rdata(PIdx::f01_Re) - p.rdata(PIdx::H01_Re)*p.rdata(PIdx::f00_Re) + p.rdata(PIdx::H01_Re)*p.rdata(PIdx::f11_Re) - p.rdata(PIdx::H11_Re)*p.rdata(PIdx::f01_Re);
-
-p.rdata(PIdx::dfdt11_Re) = 2*p.rdata(PIdx::H01_Im)*p.rdata(PIdx::f01_Re) - 2*p.rdata(PIdx::H01_Re)*p.rdata(PIdx::f01_Im);
-}
-```
+To see options for `generate_code.py`, pass the `-h` option. You do not need to
+run this script manually, when building Emu, use the `make generate
+NUM_FLAVORS=N` command to generate these source files.
diff --git a/...s/Symbolic-Hermitian-Anticommutator.ipynb → ...tians/Symbolic-Hermitian-Commutator.ipynb b/...s/Symbolic-Hermitian-Anticommutator.ipynb → ...tians/Symbolic-Hermitian-Commutator.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Symbolic Hermitian Anticommutator"
+    "# Symbolic Hermitian Commutator"
    ]
   },
   {
@@ -113,7 +113,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Calculate anticommutator $[H,f] = H \\cdot f - f \\cdot H$"
+    "Calculate commutator $[H,f] = H \\cdot f - f \\cdot H$"
    ]
   },
   {
@@ -122,7 +122,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "AC = H*F - F*H"
+    "Commutator = H*F - F*H"
    ]
   },
   {
@@ -138,7 +138,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "dFdt = sympy.I * AC"
+    "dFdt = sympy.I * Commutator"
    ]
   },
   {
@@ -263,7 +263,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,

diff --git a/Scripts/symbolic_hermitians/generate_code.py b/Scripts/symbolic_hermitians/generate_code.py
@@ -7,6 +7,7 @@
 from sympy.codegen.ast import Assignment
 from HermitianUtils import HermitianMatrix,SU_vector_ideal_magnitude
 import shutil
+import math
 
 parser = argparse.ArgumentParser(description="Generates code for calculating C = i * [A,B] for symbolic NxN Hermitian matrices A, B, C, using real-valued Real and Imaginary components.")
 parser.add_argument("N", type=int, help="Size of NxN Hermitian matrices.")
@@ -215,16 +216,26 @@ def delete_generated_files():
         U = U23*U13*U12*P
 
     # create M2 matrix in Evolve.H
-    M2 = sympy.zeros(args.N,args.N)
+    M2_massbasis = sympy.zeros(args.N,args.N)
     for i in range(args.N):
-        M2[i,i] = sympy.symbols('parms->mass'+str(i+1),real=True)**2
-    M2 = U*M2*Dagger(U)
+        M2_massbasis[i,i] = sympy.symbols('parms->mass'+str(i+1),real=True)**2
+    M2 = U*M2_massbasis*Dagger(U)
     massmatrix = HermitianMatrix(args.N, "M2matrix{}{}_{}")
     massmatrix.H = M2
     code = massmatrix.code()
     code = ["double "+code[i] for i in range(len(code))]
     write_code(code, os.path.join(args.emu_home, "Source/generated_files","Evolve.H_M2_fill"))
 
+    #=============================================#
+    # FlavoredNeutrinoContainerInit.cpp_Vvac_fill #
+    #=============================================#
+    code = []
+    massmatrix_massbasis = HermitianMatrix(args.N, "M2massbasis{}{}_{}")
+    massmatrix_massbasis.H = M2_massbasis
+    M2length = massmatrix_massbasis.SU_vector_magnitude()
+    code.append("Vvac_max = "+sympy.cxxcode(sympy.simplify(M2length))+"*PhysConst::c4/pupt_min;")
+    write_code(code, os.path.join(args.emu_home,"Source/generated_files","FlavoredNeutrinoContainerInit.cpp_Vvac_fill"))
+
     #======================#
     # Evolve.cpp_Vvac_fill #
     #======================#
@@ -247,11 +258,39 @@ def delete_generated_files():
     # Evolve.cpp_compute_dt_fill #
     #============================#
     code = []
-    for t in tails:
-        for i in range(args.N):
-            line = "N_diag_max = max(N_diag_max, state.max(GIdx::N"+str(i)+str(i)+"_Re"+t+"));"
-            code.append(line)
-    code.append("N_diag_max *= 2*"+str(args.N)+";") # overestimate of net neutrino+antineutrino number density
+    length = sympy.symbols("length",real=True)
+    cell_volume = sympy.symbols("cell_volume",real=True)
+    rho = sympy.symbols("fab(i\,j\,k\,GIdx\:\:rho)",real=True)
+    Ye = sympy.symbols("fab(i\,j\,k\,GIdx\:\:Ye)",real=True)
+    mp = sympy.symbols("PhysConst\:\:Mp",real=True)
+    sqrt2GF = sympy.symbols("M_SQRT2*PhysConst\:\:GF",real=True)
+
+    # spherically symmetric part
+    N    = HermitianMatrix(args.N, "fab(i\,j\,k\,GIdx::N{}{}_{})")
+    Nbar = HermitianMatrix(args.N, "fab(i\,j\,k\,GIdx::N{}{}_{}bar)")
+    HSI  = (N-Nbar.conjugate())
+    HSI.H[0,0] += rho*Ye/mp * cell_volume
+    V_adaptive2 = HSI.SU_vector_magnitude2()
+    code.append("V_adaptive2 += "+sympy.cxxcode(sympy.simplify(V_adaptive2))+";")
+
+    # flux part
+    for component in ["x","y","z"]:
+        F    = HermitianMatrix(args.N, "fab(i\,j\,k\,GIdx::F"+component+"{}{}_{})")
+        Fbar = HermitianMatrix(args.N, "fab(i\,j\,k\,GIdx::F"+component+"{}{}_{}bar)")
+        HSI  = (F-Fbar)
+        V_adaptive2 = HSI.SU_vector_magnitude2()
+        code.append("V_adaptive2 += "+sympy.cxxcode(sympy.simplify(V_adaptive2))+";")
+
+    # put in the units
+    code.append("V_adaptive = sqrt(V_adaptive2)*"+sympy.cxxcode(sqrt2GF/cell_volume)+";")
+
+    # old "stupid" way of computing the timestep.
+    # the factor of 2 accounts for potential worst-case effects of neutrinos and antineutrinos
+    for i in range(args.N):
+        code.append("V_stupid = max(V_stupid,"+sympy.cxxcode(N.H[i,i])+");")
+        code.append("V_stupid = max(V_stupid,"+sympy.cxxcode(Nbar.H[i,i])+");")
+    code.append("V_stupid = max(V_stupid,"+sympy.cxxcode(rho*Ye/mp*cell_volume)+");")
+    code.append("V_stupid *= "+sympy.cxxcode(2.0*args.N*sqrt2GF/cell_volume)+";")
     write_code(code, os.path.join(args.emu_home,"Source/generated_files","Evolve.cpp_compute_dt_fill"))
 
     #=======================================#
@@ -362,7 +401,12 @@ def sgn(t,var):
         for fii in fdlist:
             code.append("sumP += " + fii + ";")
         code.append("error = sumP-1.0;")
-        code.append("if( std::abs(error) > 100.*parms->maxError) amrex::Abort();")
+        code.append('if( std::abs(error) > 100.*parms->maxError) {')
+        code.append("std::ostringstream Convert;")
+        code.append('Convert << "Matrix trace (SumP) is not equal to 1, trace error exceeds 100*maxError: " << std::abs(error) << " > " << 100.*parms->maxError;')
+        code.append("std::string Trace_Error = Convert.str();")
+        code.append('amrex::Error(Trace_Error);')
+        code.append("}")
         code.append("if( std::abs(error) > parms->maxError ) {")
         for fii in fdlist:
             code.append(fii + " -= error/"+str(args.N)+";")
@@ -371,7 +415,12 @@ def sgn(t,var):
 
         # make sure diagonals are positive
         for fii in fdlist:
-            code.append("if("+fii+"<-100.*parms->maxError) amrex::Abort();")
+            code.append('if('+fii+'<-100.*parms->maxError) {')
+            code.append("std::ostringstream Convert;")
+            code.append('Convert << "Diagonal element '+fii[14:20]+' is negative, less than -100*maxError: " << '+fii+' << " < " << -100.*parms->maxError;')
+            code.append("std::string Sign_Error = Convert.str();")
+            code.append('amrex::Error(Sign_Error);')
+            code.append("}")
             code.append("if("+fii+"<-parms->maxError) "+fii+"=0;")
         code.append("")
 
@@ -381,7 +430,12 @@ def sgn(t,var):
         target_length = "p.rdata(PIdx::L"+t+")"
         code.append("length = "+sympy.cxxcode(sympy.simplify(length))+";")
         code.append("error = length-"+str(target_length)+";")
-        code.append("if( std::abs(error) > 100.*parms->maxError) amrex::Abort();")
+        code.append('if( std::abs(error) > 100.*parms->maxError) {')
+        code.append("std::ostringstream Convert;")
+        code.append('Convert << "flavor vector length differs from target length by more than 100*maxError: " << std::abs(error) << " > " << 100.*parms->maxError;')
+        code.append("std::string Length_Error = Convert.str();")
+        code.append('amrex::Error(Length_Error);')
+        code.append("}")
         code.append("if( std::abs(error) > parms->maxError) {")
         for fii in flist:
             code.append(fii+" /= length/"+str(target_length)+";")
@@ -400,4 +454,4 @@ def sgn(t,var):
     for t in tails:
         f = HermitianMatrix(args.N, "p.rdata(PIdx::f{}{}_{}"+t+")")
         code.append("p.rdata(PIdx::L"+t+") = "+sympy.cxxcode(sympy.simplify(f.SU_vector_magnitude()))+";" )
-    write_code(code, os.path.join(args.emu_home, "Source/generated_files/FlavoredNeutrinoContainerInit.cpp_set_trace_length"))
+    write_code(code, os.path.join(args.emu_home, "Source/generated_files/FlavoredNeutrinoContainerInit.cpp_set_trace_length"))
diff --git a/Source/Evolve.H b/Source/Evolve.H
@@ -26,7 +26,7 @@ namespace GIdx
     void Initialize();
 };
 
-amrex::Real compute_dt(const amrex::Geometry& geom, const amrex::Real cfl_factor, const MultiFab& state, const Real flavor_cfl_factor);
+amrex::Real compute_dt(const amrex::Geometry& geom, const amrex::Real cfl_factor, const MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, const Real flavor_cfl_factor, const Real max_adaptive_speedup);
 
 void deposit_to_mesh(const FlavoredNeutrinoContainer& neutrinos, amrex::MultiFab& state, const amrex::Geometry& geom);