From d25f5bb69eff0ddd247e3b43445a91519382eb5c Mon Sep 17 00:00:00 2001 From: Sherwood Richers Date: Fri, 17 Sep 2021 09:35:58 -0700 Subject: [PATCH] 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 * 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 * im stoopid - fixing timestep logic Co-authored-by: Nicole F <73133869+nmford20@users.noreply.github.com> Co-authored-by: Sherwood Richers Co-authored-by: Don E. Willcox --- .gitignore | 4 +- Make.Emu | 40 ++- Scripts/symbolic_hermitians/HermitianUtils.py | 49 +++- Scripts/symbolic_hermitians/README.md | 45 +-- ...nb => Symbolic-Hermitian-Commutator.ipynb} | 10 +- Scripts/symbolic_hermitians/generate_code.py | 78 ++++- Source/Evolve.H | 2 +- Source/Evolve.cpp | 56 +++- Source/FlavoredNeutrinoContainer.H | 2 + Source/FlavoredNeutrinoContainer.cpp | 2 + Source/FlavoredNeutrinoContainerInit.cpp | 216 +++++++++++++- Source/IO.H | 6 + Source/IO.cpp | 267 ++++++++++++++++++ Source/Parameters.H | 25 ++ Source/main.cpp | 9 +- sample_inputs/inputs_1d_fiducial | 1 + sample_inputs/inputs_bipolar_test | 4 +- sample_inputs/inputs_fast_flavor | 4 +- sample_inputs/inputs_fast_flavor_nonzerok | 1 + sample_inputs/inputs_msw_test | 4 +- 20 files changed, 739 insertions(+), 86 deletions(-) rename Scripts/symbolic_hermitians/{Symbolic-Hermitian-Anticommutator.ipynb => Symbolic-Hermitian-Commutator.ipynb} (97%) diff --git a/.gitignore b/.gitignore index 64aa844c..1c8c4b5f 100644 --- a/.gitignore +++ b/.gitignore @@ -37,4 +37,6 @@ Scripts/symbolic_hermitians/__pycache__ .cproject .project .pydevproject -.settings \ No newline at end of file +.settings +.vscode +Source/generated_files \ No newline at end of file diff --git a/Make.Emu b/Make.Emu index 77e39393..24d9d145 100644 --- 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 index f0c47a36..434bb797 100644 --- 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 index 84cc6a6c..4b99b1f0 100644 --- 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); -} -``` \ No newline at end of file +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/Scripts/symbolic_hermitians/Symbolic-Hermitian-Anticommutator.ipynb b/Scripts/symbolic_hermitians/Symbolic-Hermitian-Commutator.ipynb similarity index 97% rename from Scripts/symbolic_hermitians/Symbolic-Hermitian-Anticommutator.ipynb rename to Scripts/symbolic_hermitians/Symbolic-Hermitian-Commutator.ipynb index ae39193a..1cd8f0db 100644 --- a/Scripts/symbolic_hermitians/Symbolic-Hermitian-Anticommutator.ipynb +++ b/Scripts/symbolic_hermitians/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 index a00e4636..397f40f8 100755 --- 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")) \ No newline at end of file + 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 index 409a897e..ed387f1a 100644 --- 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); diff --git a/Source/Evolve.cpp b/Source/Evolve.cpp index ae2c93b0..a1fbfe98 100644 --- a/Source/Evolve.cpp +++ b/Source/Evolve.cpp @@ -19,7 +19,7 @@ namespace GIdx } } -Real compute_dt(const Geometry& geom, const Real cfl_factor, const MultiFab& state, const Real flavor_cfl_factor) +Real compute_dt(const Geometry& geom, const Real cfl_factor, const MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, const Real flavor_cfl_factor, const Real max_adaptive_speedup) { AMREX_ASSERT(cfl_factor > 0.0 || flavor_cfl_factor > 0.0); @@ -33,24 +33,54 @@ Real compute_dt(const Geometry& geom, const Real cfl_factor, const MultiFab& sta dt_translation = std::min(std::min(dxi[0],dxi[1]), dxi[2]) / PhysConst::c * cfl_factor; } - Real dt_si_matter = 0.0; + Real dt_flavor = 0.0; if (flavor_cfl_factor > 0.0) { - // self-interaction and matter part of timestep limit - // NOTE: these currently over-estimate both potentials, but avoid reduction over all particles - // NOTE: the vacuum potential is currently ignored. This requires a min reduction over particle energies - Real N_diag_max = 0; - #include "generated_files/Evolve.cpp_compute_dt_fill" - Real Vmax = std::sqrt(2.) * PhysConst::GF * std::max(N_diag_max/cell_volume, state.max(GIdx::rho)/PhysConst::Mp); - if(Vmax>0) dt_si_matter = PhysConst::hbar/Vmax*flavor_cfl_factor; + // define the reduction operator to get the max contribution to + // the potential from matter and neutrinos + // compute "effective" potential (ergs) that produces characteristic timescale + // when multiplied by hbar + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + for (MFIter mfi(state); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.fabbox(); + auto const& fab = state.array(mfi); + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + Real V_adaptive=0, V_adaptive2=0, V_stupid=0; + #include "generated_files/Evolve.cpp_compute_dt_fill" + return {V_adaptive, V_stupid}; + }); + } + + // extract the reduced values from the combined reduced data structure + auto rv = reduce_data.value(); + Real Vmax_adaptive = amrex::get<0>(rv) + FlavoredNeutrinoContainer::Vvac_max; + Real Vmax_stupid = amrex::get<1>(rv) + FlavoredNeutrinoContainer::Vvac_max; + + // reduce across MPI ranks + ParallelDescriptor::ReduceRealMax(Vmax_adaptive); + ParallelDescriptor::ReduceRealMax(Vmax_stupid ); + + // define the dt associated with each method + Real dt_flavor_adaptive = PhysConst::hbar/Vmax_adaptive*flavor_cfl_factor; + Real dt_flavor_stupid = PhysConst::hbar/Vmax_stupid *flavor_cfl_factor; + + // pick the appropriate timestep + dt_flavor = dt_flavor_stupid; + if(max_adaptive_speedup>1) + dt_flavor = min(dt_flavor_stupid*max_adaptive_speedup, dt_flavor_adaptive); } Real dt = 0.0; - if (dt_translation != 0.0 && dt_si_matter != 0.0) { - dt = std::min(dt_translation, dt_si_matter); + if (dt_translation != 0.0 && dt_flavor != 0.0) { + dt = std::min(dt_translation, dt_flavor); } else if (dt_translation != 0.0) { dt = dt_translation; - } else if (dt_si_matter != 0.0) { - dt = dt_si_matter; + } else if (dt_flavor != 0.0) { + dt = dt_flavor; } else { amrex::Error("Timestep selection failed, try using both cfl_factor and flavor_cfl_factor"); } diff --git a/Source/FlavoredNeutrinoContainer.H b/Source/FlavoredNeutrinoContainer.H index 43d4da6f..cb6bcd68 100644 --- a/Source/FlavoredNeutrinoContainer.H +++ b/Source/FlavoredNeutrinoContainer.H @@ -68,6 +68,8 @@ class FlavoredNeutrinoContainer public: + inline static Real Vvac_max; + FlavoredNeutrinoContainer(const amrex::Geometry & a_geom, const amrex::DistributionMapping & a_dmap, const amrex::BoxArray & a_ba); diff --git a/Source/FlavoredNeutrinoContainer.cpp b/Source/FlavoredNeutrinoContainer.cpp index e53c3386..ee050e2d 100644 --- a/Source/FlavoredNeutrinoContainer.cpp +++ b/Source/FlavoredNeutrinoContainer.cpp @@ -1,5 +1,7 @@ #include "FlavoredNeutrinoContainer.H" #include "Constants.H" +#include +#include using namespace amrex; diff --git a/Source/FlavoredNeutrinoContainerInit.cpp b/Source/FlavoredNeutrinoContainerInit.cpp index 8117d11b..3014b623 100644 --- a/Source/FlavoredNeutrinoContainerInit.cpp +++ b/Source/FlavoredNeutrinoContainerInit.cpp @@ -36,6 +36,49 @@ Gpu::ManagedVector > uniform_sphere_xyz(int nphi_at_equator){ return xyz; } +// residual for the root finder +// Z needs to be bigger if residual is positive +// Minerbo (1978) (unfortunately under Elsevier paywall) +// Can also see Richers (2020) https://ui.adsabs.harvard.edu/abs/2020PhRvD.102h3017R +// Eq.41 (where a is Z), but in the non-degenerate limit +// k->0, eta->0, N->Z/(4pi sinh(Z)) (just to make it integrate to 1) +// minerbo_residual is the "f" equation between eq.42 and 43 +Real minerbo_residual(const Real fluxfac, const Real Z){ + return fluxfac - 1.0/std::tanh(Z) + 1.0 / Z; +} +Real minerbo_residual_derivative(const Real fluxfac, const Real Z){ + return 1.0/(std::sinh(Z)*std::sinh(Z)) - 1.0/(Z*Z); +} +Real minerbo_Z(const Real fluxfac){ + // hard-code in these parameters because they are not + // really very important... + Real maxresidual = 1e-6; + Real maxcount = 20; + Real minfluxfac = 1e-3; + + // set the initial conditions + Real Z = 1.0; + + // catch the small flux factor case to prevent nans + if(fluxfac < minfluxfac) + Z = 3.*fluxfac; + else{ + Real residual = 1.0; + int count = 0; + while(std::abs(residual)>maxresidual and countmaxresidual) + amrex::Error("Failed to converge on a solution."); + } + + amrex::Print() << "fluxfac="< Real { return p.rdata(PIdx::pupt); }); + ParallelDescriptor::ReduceRealMin(pupt_min); + #include "generated_files/FlavoredNeutrinoContainerInit.cpp_Vvac_fill" } diff --git a/Source/IO.H b/Source/IO.H index fd57f6cd..e851f0a2 100644 --- a/Source/IO.H +++ b/Source/IO.H @@ -17,4 +17,10 @@ RecoverParticles (const std::string& dir, FlavoredNeutrinoContainer& neutrinos, amrex::Real& time, int& step); +void +writeBuildInfo (); + +void +writeJobInfo (const std::string& dir, const amrex::Geometry& geom); + #endif diff --git a/Source/IO.cpp b/Source/IO.cpp index 7f291dff..a1b09cea 100644 --- a/Source/IO.cpp +++ b/Source/IO.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "FlavoredNeutrinoContainer.H" #include "IO.H" @@ -31,6 +32,9 @@ WritePlotFile (const amrex::MultiFab& state, auto neutrino_varnames = neutrinos.get_attribute_names(); neutrinos.Checkpoint(plotfilename, "neutrinos", true, neutrino_varnames); } + + // write job information + writeJobInfo (plotfilename, geom); } void @@ -57,3 +61,266 @@ RecoverParticles (const std::string& dir, // print the step/time for the restart amrex::Print() << "Restarting after time step: " << step-1 << " t = " << time << " s. ct = " << PhysConst::c * time << " cm" << std::endl; } + + +// writeBuildInfo and writeJobInfo are copied from Castro/Source/driver/Castro_io.cpp +// and modified by Sherwood Richers +/* +SOURCE CODE LICENSE AGREEMENT +Castro, Copyright (c) 2015, +The Regents of the University of California, +through Lawrence Berkeley National Laboratory +(subject to receipt of any required approvals from the U.S. +Dept. of Energy). All rights reserved." + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +(1) Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +(2) Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +(3) Neither the name of the University of California, Lawrence +Berkeley National Laboratory, U.S. Dept. of Energy nor the names of +its contributors may be used to endorse or promote products derived +from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +You are under no obligation whatsoever to provide any bug fixes, +patches, or upgrades to the features, functionality or performance of +the source code ("Enhancements") to anyone; however, if you choose to +make your Enhancements available either publicly, or directly to +Lawrence Berkeley National Laboratory, without imposing a separate +written license agreement for such Enhancements, then you hereby grant +the following license: a non-exclusive, royalty-free perpetual +license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such +enhancements or derivative works thereof, in binary and source code +form. +*/ +void +writeBuildInfo (){ + std::string PrettyLine = std::string(78, '=') + "\n"; + std::string OtherLine = std::string(78, '-') + "\n"; + std::string SkipSpace = std::string(8, ' '); + + // build information + std::cout << PrettyLine; + std::cout << " Emu Build Information\n"; + std::cout << PrettyLine; + + std::cout << "build date: " << buildInfoGetBuildDate() << "\n"; + std::cout << "build machine: " << buildInfoGetBuildMachine() << "\n"; + std::cout << "build dir: " << buildInfoGetBuildDir() << "\n"; + std::cout << "AMReX dir: " << buildInfoGetAMReXDir() << "\n"; + + std::cout << "\n"; + + std::cout << "COMP: " << buildInfoGetComp() << "\n"; + std::cout << "COMP version: " << buildInfoGetCompVersion() << "\n"; + + std::cout << "\n"; + + std::cout << "C++ compiler: " << buildInfoGetCXXName() << "\n"; + std::cout << "C++ flags: " << buildInfoGetCXXFlags() << "\n"; + + std::cout << "\n"; + + std::cout << "Link flags: " << buildInfoGetLinkFlags() << "\n"; + std::cout << "Libraries: " << buildInfoGetLibraries() << "\n"; + + std::cout << "\n"; + + for (int n = 1; n <= buildInfoGetNumModules(); n++) { + std::cout << buildInfoGetModuleName(n) << ": " << buildInfoGetModuleVal(n) << "\n"; + } + + std::cout << "\n"; + + const char* githash1 = buildInfoGetGitHash(1); + const char* githash2 = buildInfoGetGitHash(2); + if (strlen(githash1) > 0) { + std::cout << "Emu git describe: " << githash1 << "\n"; + } + if (strlen(githash2) > 0) { + std::cout << "AMReX git describe: " << githash2 << "\n"; + } + + const char* buildgithash = buildInfoGetBuildGitHash(); + const char* buildgitname = buildInfoGetBuildGitName(); + if (strlen(buildgithash) > 0){ + std::cout << buildgitname << " git describe: " << buildgithash << "\n"; + } + + std::cout << "\n"< 0) { + jobInfoFile << "Emu git describe: " << githash1 << "\n"; + } + if (strlen(githash2) > 0) { + jobInfoFile << "AMReX git describe: " << githash2 << "\n"; + } + + const char* buildgithash = buildInfoGetBuildGitHash(); + const char* buildgitname = buildInfoGetBuildGitName(); + if (strlen(buildgithash) > 0){ + jobInfoFile << buildgitname << " git describe: " << buildgithash << "\n"; + } + + jobInfoFile << "\n\n"; + + + // grid information + jobInfoFile << PrettyLine; + jobInfoFile << " Grid Information\n"; + jobInfoFile << PrettyLine; + + jobInfoFile << "geometry.is_periodic: "; + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + jobInfoFile << geom.isPeriodic(idir) << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << "geometry.coord_sys: " << geom.Coord() << "\n"; + + jobInfoFile << "geometry.prob_lo: "; + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + jobInfoFile << geom.ProbLo(idir) << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << "geometry.prob_hi: "; + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + jobInfoFile << geom.ProbHi(idir) << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << "amr.n_cell: "; + const int* domain_lo = geom.Domain().loVect(); + const int* domain_hi = geom.Domain().hiVect(); + for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { + jobInfoFile << domain_hi[idir] - domain_lo[idir] + 1 << " "; + } + jobInfoFile << "\n\n"; + + jobInfoFile.close(); + + // now the external parameters + const int jobinfo_file_length = FullPathJobInfoFile.length(); + Vector jobinfo_file_name(jobinfo_file_length); + + for (int i = 0; i < jobinfo_file_length; i++) { + jobinfo_file_name[i] = FullPathJobInfoFile[i]; + } + +} diff --git a/Source/Parameters.H b/Source/Parameters.H index 26fbb946..498b6d91 100644 --- a/Source/Parameters.H +++ b/Source/Parameters.H @@ -22,6 +22,7 @@ struct TestParams : public amrex::Gpu::Managed Real rho_in, Ye_in, T_in; // g/ccm, 1, MeV int simulation_type; Real cfl_factor, flavor_cfl_factor; + Real max_adaptive_speedup; bool do_restart; std::string restart_dir; Real maxError; @@ -41,6 +42,13 @@ struct TestParams : public amrex::Gpu::Managed Real st4_ndensbar, st4_thetabar, st4_phibar, st4_fluxfacbar; Real st4_amplitude; + // simulation_type==5 + Real st5_nnue , st5_nnua , st5_nnux ; + Real st5_fxnue, st5_fxnua, st5_fxnux; + Real st5_fynue, st5_fynua, st5_fynux; + Real st5_fznue, st5_fznua, st5_fznux; + Real st5_amplitude; + void Initialize(){ ParmParse pp; pp.get("simulation_type", simulation_type); @@ -58,6 +66,7 @@ struct TestParams : public amrex::Gpu::Managed pp.get("T_MeV", T_in); pp.get("cfl_factor", cfl_factor); pp.get("flavor_cfl_factor", flavor_cfl_factor); + pp.get("max_adaptive_speedup", max_adaptive_speedup); pp.get("write_plot_every", write_plot_every); pp.get("write_plot_particles_every", write_plot_particles_every); pp.get("do_restart", do_restart); @@ -103,6 +112,22 @@ struct TestParams : public amrex::Gpu::Managed pp.get("st4_fluxfacbar", st4_fluxfacbar); pp.get("st4_amplitude", st4_amplitude); } + + if(simulation_type==5){ + pp.get("st5_nnue",st5_nnue); + pp.get("st5_nnua",st5_nnua); + pp.get("st5_nnux",st5_nnux); + pp.get("st5_fxnue",st5_fxnue); + pp.get("st5_fxnua",st5_fxnua); + pp.get("st5_fxnux",st5_fxnux); + pp.get("st5_fynua",st5_fynua); + pp.get("st5_fynux",st5_fynux); + pp.get("st5_fynue",st5_fynue); + pp.get("st5_fznue",st5_fznue); + pp.get("st5_fznua",st5_fznua); + pp.get("st5_fznux",st5_fznux); + pp.get("st5_amplitude",st5_amplitude); + } } }; diff --git a/Source/main.cpp b/Source/main.cpp index 309d5e63..edd2d6c3 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -194,7 +194,7 @@ void evolve_flavor(const TestParams* parms) // Note: this won't be the same as the new-time grid data // because the last deposit_to_mesh call was at either the old time (forward Euler) // or the final RK stage, if using Runge-Kutta. - const Real dt = compute_dt(geom,parms->cfl_factor,state,parms->flavor_cfl_factor); + const Real dt = compute_dt(geom,parms->cfl_factor,state,neutrinos,parms->flavor_cfl_factor,parms->max_adaptive_speedup); integrator.set_timestep(dt); }; @@ -203,7 +203,7 @@ void evolve_flavor(const TestParams* parms) integrator.set_post_timestep(post_timestep_fun); // Get a starting timestep - const Real starting_dt = compute_dt(geom,parms->cfl_factor,state,parms->flavor_cfl_factor); + const Real starting_dt = compute_dt(geom,parms->cfl_factor,state,neutrinos_old,parms->flavor_cfl_factor, parms->max_adaptive_speedup); // Do all the science! amrex::Print() << "Starting timestepping loop... " << std::endl; @@ -230,6 +230,11 @@ int main(int argc, char* argv[]) { amrex::Initialize(argc,argv); + // write build information to screen + if (ParallelDescriptor::IOProcessor()) { + writeBuildInfo(); + } + // by default amrex initializes rng deterministically // this uses the time for a different run each time amrex::InitRandom(ParallelDescriptor::MyProc()+time(NULL), ParallelDescriptor::NProcs()); diff --git a/sample_inputs/inputs_1d_fiducial b/sample_inputs/inputs_1d_fiducial index 1af1dc46..72ee3bfb 100644 --- a/sample_inputs/inputs_1d_fiducial +++ b/sample_inputs/inputs_1d_fiducial @@ -11,6 +11,7 @@ st4_amplitude = 1e-6 cfl_factor = 0.5 flavor_cfl_factor = 0.5 +max_adaptive_speedup = 0 maxError = 1e-6 integration.type = 1 diff --git a/sample_inputs/inputs_bipolar_test b/sample_inputs/inputs_bipolar_test index e8f043fb..305ddf47 100644 --- a/sample_inputs/inputs_bipolar_test +++ b/sample_inputs/inputs_bipolar_test @@ -1,5 +1,6 @@ simulation_type = 1 cfl_factor = 0.5 +max_adaptive_speedup = 0 flavor_cfl_factor = .5 maxError = 1e-6 @@ -55,7 +56,8 @@ mass1_eV = -0.008596511 mass2_eV = 0 # mass state 3 mass in eV [NO:sqrt(2.449e-3) IO:-sqrt(2.509e-3)] -mass3_eV = 0.049487372 +#mass3_eV = 0.049487372 +mass3_eV = 0 # 1-2 mixing angle in degrees [NO/IO:33.82] theta12_degrees = 33.82 diff --git a/sample_inputs/inputs_fast_flavor b/sample_inputs/inputs_fast_flavor index dd994cba..adc998b1 100644 --- a/sample_inputs/inputs_fast_flavor +++ b/sample_inputs/inputs_fast_flavor @@ -1,6 +1,7 @@ simulation_type = 2 cfl_factor = 0.5 flavor_cfl_factor = .5 +max_adaptive_speedup = 0 maxError = 1e-6 integration.type = 1 @@ -55,7 +56,8 @@ mass1_eV = -0.008596511 mass2_eV = 0 # mass state 3 mass in eV [NO:sqrt(2.449e-3) IO:-sqrt(2.509e-3)] -mass3_eV = 0.049487372 +#mass3_eV = 0.049487372 +mass3_eV = 0 # 1-2 mixing angle in degrees [NO/IO:33.82] theta12_degrees = 1e-6 diff --git a/sample_inputs/inputs_fast_flavor_nonzerok b/sample_inputs/inputs_fast_flavor_nonzerok index 49c82619..32ccb2ae 100644 --- a/sample_inputs/inputs_fast_flavor_nonzerok +++ b/sample_inputs/inputs_fast_flavor_nonzerok @@ -4,6 +4,7 @@ st3_wavelength_fraction_of_domain = 1 cfl_factor = 0.5 flavor_cfl_factor = 0.5 +max_adaptive_speedup = 0 maxError = 1e-6 integration.type = 1 diff --git a/sample_inputs/inputs_msw_test b/sample_inputs/inputs_msw_test index a8231fb5..e674a58c 100644 --- a/sample_inputs/inputs_msw_test +++ b/sample_inputs/inputs_msw_test @@ -1,6 +1,7 @@ simulation_type = 0 cfl_factor = 0.5 flavor_cfl_factor = 0.1 +max_adaptive_speedup = 0 maxError = 1e-6 integration.type = 1 @@ -55,7 +56,8 @@ mass1_eV = -0.008596511 mass2_eV = 0 # mass state 3 mass in eV [NO:sqrt(2.449e-3) IO:-sqrt(2.509e-3)] -mass3_eV = 0.049487372 +# mass3_eV = 0.049487372 +mass3_eV = 0 # 1-2 mixing angle in degrees [NO/IO:33.82] theta12_degrees = 33.82