added travis test case for nonzero wavevector fast flavor instability

AMReX-Astro · Sep 2, 2020 · 952ad47 · 952ad47
1 parent dbdb37c
commit 952ad47
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 59 deletions.
diff --git 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
+- 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
@@ -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
@@ -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]);
 		}