diff --git a/.travis.yml b/.travis.yml index 34211cb4..bd44f824 100644 --- a/.travis.yml +++ b/.travis.yml @@ -32,4 +32,5 @@ install: script: - cd Exec; mpirun -np 2 ./main3d.gnu.DEBUG.TPROF.MPI.ex ../sample_inputs/inputs_msw_test; python ../Scripts/tests/msw_test.py - cd Exec; mpirun -np 2 ./main3d.gnu.DEBUG.TPROF.MPI.ex ../sample_inputs/inputs_bipolar_test -- cd Exec; mpirun -np 2 ./main3d.gnu.DEBUG.TPROF.MPI.ex ../sample_inputs/inputs_fast_flavor; python ../Scripts/tests/fast_flavor_test.py \ No newline at end of file +- cd Exec; mpirun -np 2 ./main3d.gnu.DEBUG.TPROF.MPI.ex ../sample_inputs/inputs_fast_flavor; python ../Scripts/tests/fast_flavor_test.py +- cd Exec; mpirun -np 2 ./main3d.gnu.DEBUG.TPROF.MPI.ex ../sample_inputs/inputs_fast_flavor_nonzerok; python ../Scripts/tests/fast_flavor_k_test.py diff --git a/Scripts/tests/fast_flavor_k_test.py b/Scripts/tests/fast_flavor_k_test.py index 7d9b2c55..0d86c893 100644 --- a/Scripts/tests/fast_flavor_k_test.py +++ b/Scripts/tests/fast_flavor_k_test.py @@ -2,21 +2,19 @@ import argparse import glob import EmuReader +import sys +import os +importpath = os.path.dirname(os.path.realpath(__file__))+"/../visualization/" +sys.path.append(importpath) +import amrex_plot_tools as amrex parser = argparse.ArgumentParser() parser.add_argument("-na", "--no_assert", action="store_true", help="If --no_assert is supplied, do not raise assertion errors if the test error > tolerance.") args = parser.parse_args() -# physical constants -clight = 2.99792458e10 # cm/s -hbar = 1.05457266e-27 # erg s -theta12 = 33.82*np.pi/180. # radians -eV = 1.60218e-12 # erg -dm21c4 = 7.39e-5 * eV**2 # erg^2 -mp = 1.6726219e-24 # g -GF = 1.1663787e-5 / (1e9*eV)**2 * (hbar*clight)**3 #erg cm^3 -tolerance = 1e-3 +dm21c4 = 0 # parameter file has both masses set to 0 +tolerance = 1e-2 i0 = 100 i1 = 160 @@ -28,30 +26,8 @@ if __name__ == "__main__": - rkey = { - "x":0, - "y":1, - "z":2, - "time":3, - "pupx":4, - "pupy":5, - "pypz":6, - "pupt":7, - "N":8, - "f00_Re":9, - "f01_Re":10, - "f01_Im":11, - "f11_Re":12, - "Nbar":13, - "f00_Rebar":14, - "f01_Rebar":15, - "f01_Imbar":16, - "f11_Rebar":17, - } - ikey = { - # no ints are stored - } - + rkey, ikey = amrex.get_particle_keys() + t = [] fexR = [] fexI = [] @@ -64,45 +40,42 @@ plotfile = "plt"+str(i).zfill(5) idata, rdata = EmuReader.read_particle_data(plotfile, ptype="neutrinos") - p = rdata[0] - t.append(p[rkey["time"]]) - fexR.append(p[rkey["f01_Re"]]) - fexI.append(p[rkey["f01_Im"]]) - fexRbar.append(p[rkey["f01_Rebar"]]) - fexIbar.append(p[rkey["f01_Imbar"]]) - pupt.append(p[rkey["pupt"]]) + p = rdata + t.append(p[0][rkey["time"]]) + fexR.append(np.max(np.abs(p[:,rkey["f01_Re"]]))) + fexI.append(np.max(np.abs(p[:,rkey["f01_Im"]]))) + fexRbar.append(np.max(np.abs(p[:,rkey["f01_Rebar"]]))) + fexIbar.append(np.max(np.abs(p[:,rkey["f01_Imbar"]]))) + pupt.append(p[0][rkey["pupt"]]) t = np.array(t) fexR = np.array(fexR) fexI = np.array(fexI) fexRbar = np.array(fexRbar) fexIbar = np.array(fexIbar) - + print(fexR) + # The neutrino energy we set - E = 50. * 1e6*eV + E = 50. * 1e6*amrex.eV V = dm21c4 / (2.*E) # wavevector of the perturbation k = 2*np.pi/Lx - - # SI potential - mu = 1.585229243e-24/hbar # 1/s - print(mu) + print("k=",k) # theoretical growth rate according to Chakraborty 2016 Equation 2.7 a=0 mu=0.5(V+k) - ImOmega = (V+k*hbar*clight)/hbar + ImOmega = (V+k*amrex.hbar*amrex.clight)/amrex.hbar print("Theoretical growth rate:",ImOmega," s^-1") - print(V,k*hbar*clight) # get growth rate from each diagonal component dt = t[i1]-t[i0] - fexRomega = np.log(fexR[i1]/fexR[i0]) / dt - fexIomega = np.log(fexI[i1]/fexI[i0]) / dt - fexRbaromega = np.log(fexRbar[i1]/fexRbar[i0]) / dt - fexIbaromega = np.log(fexIbar[i1]/fexIbar[i0]) / dt + fexRomega = np.log(np.abs(fexR[i1]/fexR[i0])) / dt + fexIomega = np.log(np.abs(fexI[i1]/fexI[i0])) / dt + fexRbaromega = np.log(np.abs(fexRbar[i1]/fexRbar[i0])) / dt + fexIbaromega = np.log(np.abs(fexIbar[i1]/fexIbar[i0])) / dt print("growth rates:",fexRomega,fexIomega,fexRbaromega,fexIbaromega) - print("growth rates:",fexRomega/ImOmega,fexIomega/ImOmega,fexRbaromega/ImOmega,fexIbaromega/ImOmega) + print("growth rates / theoretical:",fexRomega/ImOmega,fexIomega/ImOmega,fexRbaromega/ImOmega,fexIbaromega/ImOmega) def myassert(condition): if not args.no_assert: diff --git a/Source/FlavoredNeutrinoContainerInit.cpp b/Source/FlavoredNeutrinoContainerInit.cpp index 8cd98b55..b88d6bc3 100644 --- a/Source/FlavoredNeutrinoContainerInit.cpp +++ b/Source/FlavoredNeutrinoContainerInit.cpp @@ -324,7 +324,6 @@ InitParticles(const TestParams* parms) Real amplitude = 1e-6; Real lambda = Geom(lev).ProbLength(dir); Real k = (2.*M_PI) / lambda; - Print() << "k=" << k*PhysConst::hbarc << " erg" << std::endl; // Set particle flavor p.rdata(PIdx::f00_Re) = 1.0; @@ -347,13 +346,9 @@ InitParticles(const TestParams* parms) // to get maximal growth according to Chakraborty 2016 Equation 2.10 Real dm2 = (parms->mass2-parms->mass1)*(parms->mass2-parms->mass1); //g^2 Real omega = dm2*PhysConst::c4 / (2.* p.rdata(PIdx::pupt)); - Print() << "omega=" << omega << " erg" << std::endl; Real mu_ndens = sqrt(2.) * PhysConst::GF; // SI potential divided by the number density Real ndens = (omega+k*PhysConst::hbarc) / (2.*mu_ndens); // want omega/2mu to be 1 Real mu = mu_ndens*ndens; - Print() << "omega+k=" << omega+PhysConst::hbarc*k << std::endl; - Print() << "mu=" << mu << " erg" << std::endl; - Print() << "mu/omegatilde="<< mu/(omega+PhysConst::hbarc*k) << std::endl; p.rdata(PIdx::N) = ndens * scale_fac * (1. + u[2]); p.rdata(PIdx::Nbar) = ndens * scale_fac * (1. - u[2]); }