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