Merge pull request #52 from AMReX-Astro/development

Optimizations for particle communication and bugfixes

.gitmodules

@@ -1,3 +1,3 @@
 [submodule "amrex"]
 	path = amrex
-	url = https://github.com/srichers/amrex.git
+	url = https://github.com/dwillcox/amrex.git

CHANGES.md

@@ -0,0 +1,29 @@
+
+# 1.2
+
+   * Using the new particle shape factor implementation, we now redistribute
+     particles only at the end of each timestep instead of at every runge-kutta
+     stage (#50).
+
+   * Removed WarpX bilinear filter, since we don't need it to get stable
+     evolution. Emu now does not include any WarpX code (#49).
+
+   * Removed WarpX shape factors. Wrote new shape factor implementation in
+     `ParticleInterpolator` so we now can evolve particles within ghost cells (#48).
+
+   * We now determine the shape factor order taking into account both the
+     user-specified order and the dimensionality of the problem so we use order
+     0 in the directions where the domain contains only one cell (#47).
+
+   * Fixed issue with flavor vector length so we can now run simulations with
+     nonzero initial amounts of non-electron flavor neutrinos (#44).
+
+# 1.1.1
+
+   * Emu has been through the LBNL licensing process for open source code. We
+     release Emu under the BSD-3-Clause-LBNL license, the code is identical to v1.0.
+
+# 1.0
+
+   * Version of Emu used for paper 1 _Particle-in-cell Simulation of the Neutrino Fast Flavor Instability_ (https://arxiv.org/abs/2101.02745)
+

Exec/GNUmakefile

@@ -1,8 +1,7 @@
 EMU_HOME ?= ../
 AMREX_HOME ?= ../amrex
 
-DIM = 3
-
+SHAPE_FACTOR_ORDER ?= 2
 NUM_FLAVORS ?= 3
 
 COMP = gnu

Make.Emu

@@ -1,5 +1,6 @@
 NUM_FLAVORS ?= 2
 SHAPE_FACTOR_ORDER ?= 2
+DIM = 3
 
 TOP := $(EMU_HOME)
 

Scripts/symbolic_hermitians/HermitianUtils.py

@@ -93,6 +93,12 @@ def SU_vector_magnitude(self):
 
         return sympy.sqrt(mag2)
 
+    def trace(self):
+        result = 0
+        for i in range(self.size):
+            result += self.H[i,i]
+        return result
+
     def times(self, x):
         # Apply self.H = self.H * x
         # where x is a Sympy expression

Scripts/symbolic_hermitians/generate_code.py

@@ -150,16 +150,16 @@ def delete_generated_files():
     # Evolve.cpp_deposit_to_mesh_fill #
     #=================================#
     tails = ["","bar"]
-    string1 = "amrex::Gpu::Atomic::Add(&sarr(i+ii-1, j+jj-1, k+kk-1, GIdx::"
-    string2 = "-start_comp), sx[ii]*sy[jj]*sz[kk] * p.rdata(PIdx::"
+    string1 = "amrex::Gpu::Atomic::AddNoRet(&sarr(i, j, k, GIdx::"
+    string2 = "-start_comp), sx(i) * sy(j) * sz(k) * p.rdata(PIdx::"
     string4 = [");",
                "*p.rdata(PIdx::pupx)/p.rdata(PIdx::pupt));",
                "*p.rdata(PIdx::pupy)/p.rdata(PIdx::pupt));",
                "*p.rdata(PIdx::pupz)/p.rdata(PIdx::pupt));"]
     deposit_vars = ["N","Fx","Fy","Fz"]
     code = []
     for t in tails:
-        string3 = ")*2.*p.rdata(PIdx::L"+t+")"
+        string3 = ")*p.rdata(PIdx::N"+t+")"
         flist = HermitianMatrix(args.N, "f{}{}_{}"+t).header()
         for ivar in range(len(deposit_vars)):
             deplist = HermitianMatrix(args.N, deposit_vars[ivar]+"{}{}_{}"+t).header()
@@ -260,8 +260,8 @@ def delete_generated_files():
     # matter and SI potentials require interpolating from grid
     tails = ["","bar"]
     string1 = "p.rdata(PIdx::"
-    string2 = ") +=  sqrt(2.) * PhysConst::GF * sx[ii]*sy[jj]*sz[kk] * ("
-    string_interp = "sarr(i+ii-1,j+jj-1,k+kk-1,GIdx::"
+    string2 = ") +=  sqrt(2.) * PhysConst::GF * sx(i) * sy(j) * sz(k) * ("
+    string_interp = "sarr(i, j, k, GIdx::"
     direction = ["x","y","z"]
     string3 = ["*p.rdata(PIdx::pupx)"]
     string4 = "/p.rdata(PIdx::pupt)"
@@ -312,7 +312,7 @@ def sgn(t,var):
             else:
                 line += " -= "
 
-            line += "sqrt(2.) * PhysConst::GF * inv_cell_volume * sx[ii]*sy[jj]*sz[kk] * (inside_parentheses);"
+            line += "sqrt(2.) * PhysConst::GF * inv_cell_volume * sx(i) * sy(j) * sz(k) * (inside_parentheses);"
             code.append(line)
             code.append("")
     write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_interpolate_from_mesh_fill"))
@@ -378,7 +378,7 @@ def sgn(t,var):
         # make sure the flavor vector length is what it would be with a 1 in only one diagonal
         length = sympy.symbols("length",real=True)
         length = f.SU_vector_magnitude()
-        target_length = SU_vector_ideal_magnitude(args.N)
+        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();")
@@ -392,3 +392,12 @@ def sgn(t,var):
     # Write code to output file, using a template if one is provided
     # write_code(code, "code.cpp", args.output_template)
 
+
+    #====================================================#
+    # FlavoredNeutrinoContainerInit.cpp_set_trace_length #
+    #====================================================#
+    code = []
+    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"))

Source/Evolve.H

@@ -7,7 +7,6 @@
 #include <AMReX_MultiFab.H>
 #include <AMReX_ParticleMesh.H>
 #include <FlavoredNeutrinoContainer.H>
-#include <BilinearFilter.H>
 
 namespace GIdx
 {
@@ -29,7 +28,7 @@ namespace GIdx
 
 amrex::Real compute_dt(const amrex::Geometry& geom, const amrex::Real cfl_factor, const MultiFab& state, const Real flavor_cfl_factor);
 
-void deposit_to_mesh(const FlavoredNeutrinoContainer& neutrinos, amrex::MultiFab& state, const amrex::Geometry& geom, BilinearFilter& bilinear_filter, bool use_filter);
+void deposit_to_mesh(const FlavoredNeutrinoContainer& neutrinos, amrex::MultiFab& state, const amrex::Geometry& geom);
 
 void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const amrex::MultiFab& state, const amrex::Geometry& geom, const TestParams* parms);