Flux injection from EB: Pick a random point instead of the centroid (#…

…5493)

Co-authored-by: Remi Lehe <[email protected]>

Examples/Tests/flux_injection/CMakeLists.txt

@@ -26,8 +26,8 @@ add_warpx_test(
     3  # dims
     2  # nprocs
     inputs_test_3d_flux_injection_from_eb  # inputs
-    "analysis_flux_injection_from_eb.py diags/diag1000010"  # analysis
-    "analysis_default_regression.py --path diags/diag1000010"  # checksum
+    "analysis_flux_injection_from_eb.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
     OFF  # dependency
 )
 
@@ -36,8 +36,8 @@ add_warpx_test(
     RZ  # dims
     2  # nprocs
     inputs_test_rz_flux_injection_from_eb  # inputs
-    "analysis_flux_injection_from_eb.py diags/diag1000010"  # analysis
-    "analysis_default_regression.py --path diags/diag1000010"  # checksum
+    "analysis_flux_injection_from_eb.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
     OFF  # dependency
 )
 
@@ -46,7 +46,7 @@ add_warpx_test(
     2  # dims
     2  # nprocs
     inputs_test_2d_flux_injection_from_eb  # inputs
-    "analysis_flux_injection_from_eb.py diags/diag1000010"  # analysis
-    "analysis_default_regression.py --path diags/diag1000010"  # checksum
+    "analysis_flux_injection_from_eb.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
     OFF  # dependency
 )

Examples/Tests/flux_injection/analysis_flux_injection_from_eb.py

@@ -147,7 +147,8 @@ def compare_gaussian_flux(u, w, u_th, u_m, label=""):
 wy = nz * vx - nx * vz
 wz = nx * vy - ny * vx
 u_perp2 = ux * wx + uy * wy + uz * wz
-compare_gaussian(u_perp2, w, u_th=0.01, label="u_perp")
+compare_gaussian(u_perp2, w, u_th=0.01, label="u_perp2")
 
+plt.legend()
 plt.tight_layout()
 plt.savefig("Distribution.png")

Examples/Tests/flux_injection/inputs_base_from_eb

@@ -1,5 +1,5 @@
 # Maximum number of time steps
-max_step = 10
+max_step = 20
 
 # The lo and hi ends of grids are multipliers of blocking factor
 amr.blocking_factor = 8
@@ -13,7 +13,7 @@ amr.max_level = 0
 
 # Deactivate Maxwell solver
 algo.maxwell_solver = none
-warpx.const_dt = 1e-9
+warpx.const_dt = 0.5e-9
 
 # Embedded boundary
 warpx.eb_implicit_function = "-(x**2+y**2+z**2-2**2)"

Examples/Tests/flux_injection/inputs_test_2d_flux_injection_from_eb

@@ -1,7 +1,7 @@
 FILE = inputs_base_from_eb
 
 # number of grid points
-amr.n_cell = 16 16
+amr.n_cell = 32 32
 
 # Geometry
 geometry.dims = 2

Examples/Tests/flux_injection/inputs_test_3d_flux_injection_from_eb

@@ -1,7 +1,7 @@
 FILE = inputs_base_from_eb
 
 # number of grid points
-amr.n_cell = 16 16 16
+amr.n_cell = 32 32 32
 
 # Geometry
 geometry.dims = 3

Examples/Tests/flux_injection/inputs_test_rz_flux_injection_from_eb

@@ -1,7 +1,7 @@
 FILE = inputs_base_from_eb
 
 # number of grid points
-amr.n_cell = 8 16
+amr.n_cell = 16 32
 
 # Geometry
 geometry.dims = RZ

Regression/Checksum/benchmarks_json/test_2d_flux_injection_from_eb.json

@@ -1,11 +1,11 @@
 {
   "lev=0": {},
   "electron": {
-    "particle_momentum_x": 3.4911323396038835e-19,
-    "particle_momentum_y": 2.680312173420972e-20,
-    "particle_momentum_z": 3.4918430443688734e-19,
-    "particle_position_x": 17950.08139982036,
-    "particle_position_y": 17949.47183079554,
-    "particle_weight": 6.269e-08
+    "particle_momentum_x": 1.4013860393698154e-18,
+    "particle_momentum_y": 1.0934049057929508e-19,
+    "particle_momentum_z": 1.4066623146535866e-18,
+    "particle_position_x": 72129.9049362857,
+    "particle_position_y": 72178.76505490103,
+    "particle_weight": 6.279375e-08
   }
 }

Regression/Checksum/benchmarks_json/test_3d_flux_injection_from_eb.json

@@ -1,12 +1,12 @@
 {
   "lev=0": {},
   "electron": {
-    "particle_momentum_x": 2.1855512033870577e-18,
-    "particle_momentum_y": 2.1826030840183147e-18,
-    "particle_momentum_z": 2.181852403122796e-18,
-    "particle_position_x": 111042.81925863726,
-    "particle_position_y": 111012.52928910403,
-    "particle_position_z": 111015.90903542604,
-    "particle_weight": 2.4775750000000003e-07
+    "particle_momentum_x": 1.7587772989573373e-17,
+    "particle_momentum_y": 1.7608560965806728e-17,
+    "particle_momentum_z": 1.7596701993624562e-17,
+    "particle_position_x": 902783.9285213391,
+    "particle_position_y": 902981.7980528818,
+    "particle_position_z": 902777.1246066706,
+    "particle_weight": 2.503818749999996e-07
   }
-}
+}

Regression/Checksum/benchmarks_json/test_rz_flux_injection_from_eb.json

@@ -1,12 +1,12 @@
 {
   "lev=0": {},
   "electron": {
-    "particle_momentum_x": 3.3665608248716305e-19,
-    "particle_momentum_y": 3.392690322852239e-19,
-    "particle_momentum_z": 5.254577143779578e-19,
-    "particle_position_x": 26933.772112044953,
-    "particle_position_y": 26926.994273876346,
-    "particle_theta": 29492.77423173835,
-    "particle_weight": 2.4953304765944705e-07
+    "particle_momentum_x": 1.3547613622259754e-18,
+    "particle_momentum_y": 1.3539614160696825e-18,
+    "particle_momentum_z": 2.102305484242805e-18,
+    "particle_position_x": 108281.74349700565,
+    "particle_position_y": 108222.91506078152,
+    "particle_theta": 118597.06004310239,
+    "particle_weight": 2.5087578786544294e-07
   }
-}
+}

Source/Particles/AddPlasmaUtilities.H

@@ -111,28 +111,24 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 amrex::Real compute_scale_fac_area_eb (
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
     const amrex::Real num_ppc_real,
-    amrex::Array4<const amrex::Real> const& eb_bnd_normal_arr,
-    int i, int j, int k ) {
+    AMREX_D_DECL(const amrex::Real n0,
+                 const amrex::Real n1,
+                 const amrex::Real n2)) {
     using namespace amrex::literals;
     // Scale particle weight by the area of the emitting surface, within one cell
     // By definition, eb_bnd_area_arr is normalized (unitless).
     // Here we undo the normalization (i.e. multiply by the surface used for normalization in amrex:
     // see https://amrex-codes.github.io/amrex/docs_html/EB.html#embedded-boundary-data-structures)
 #if defined(WARPX_DIM_3D)
-    const amrex::Real nx = eb_bnd_normal_arr(i,j,k,0);
-    const amrex::Real ny = eb_bnd_normal_arr(i,j,k,1);
-    const amrex::Real nz = eb_bnd_normal_arr(i,j,k,2);
-    amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(nx*dx[1]*dx[2]) +
-                  amrex::Math::powi<2>(ny*dx[0]*dx[2]) +
-                  amrex::Math::powi<2>(nz*dx[0]*dx[1]));
+    amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(n0*dx[1]*dx[2]) +
+                                      amrex::Math::powi<2>(n1*dx[0]*dx[2]) +
+                                      amrex::Math::powi<2>(n2*dx[0]*dx[1]));
 
 #elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
-    const amrex::Real nx = eb_bnd_normal_arr(i,j,k,0);
-    const amrex::Real nz = eb_bnd_normal_arr(i,j,k,1);
-    amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(nx*dx[1]) +
-                  amrex::Math::powi<2>(nz*dx[0]));
+    amrex::Real scale_fac = std::sqrt(amrex::Math::powi<2>(n0*dx[1]) +
+                                      amrex::Math::powi<2>(n1*dx[0]));
 #else
-    amrex::ignore_unused(dx, eb_bnd_normal_arr, i, j, k);
+    amrex::ignore_unused(dx, AMREX_D_DECL(n0,n1,n2));
     amrex::Real scale_fac = 0.0_rt;
 #endif
     // Do not multiply by eb_bnd_area_arr(i,j,k) here because this
@@ -159,18 +155,19 @@ amrex::Real compute_scale_fac_area_eb (
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void rotate_momentum_eb (
     PDim3 & pu,
-    amrex::Array4<const amrex::Real> const& eb_bnd_normal_arr,
-    int i, int j, int k )
+    AMREX_D_DECL(const amrex::Real n0,
+                 const amrex::Real n1,
+                 const amrex::Real n2))
 {
     using namespace amrex::literals;
 
 #if defined(WARPX_DIM_3D)
     // The minus sign below takes into account the fact that eb_bnd_normal_arr
     // points towards the covered region, while particles are to be emitted
     // *away* from the covered region.
-    amrex::Real const nx = -eb_bnd_normal_arr(i,j,k,0);
-    amrex::Real const ny = -eb_bnd_normal_arr(i,j,k,1);
-    amrex::Real const nz = -eb_bnd_normal_arr(i,j,k,2);
+    amrex::Real const nx = -n0;
+    amrex::Real const ny = -n1;
+    amrex::Real const nz = -n2;
 
     // Rotate the momentum in theta and phi
     amrex::Real const cos_theta = nz;
@@ -194,14 +191,14 @@ void rotate_momentum_eb (
     // The minus sign below takes into account the fact that eb_bnd_normal_arr
     // points towards the covered region, while particles are to be emitted
     // *away* from the covered region.
-    amrex::Real const sin_theta = -eb_bnd_normal_arr(i,j,k,0);
-    amrex::Real const cos_theta = -eb_bnd_normal_arr(i,j,k,1);
+    amrex::Real const sin_theta = -n0;
+    amrex::Real const cos_theta = -n1;
     amrex::Real const uz = pu.z*cos_theta - pu.x*sin_theta;
     amrex::Real const ux = pu.x*cos_theta + pu.z*sin_theta;
     pu.x = ux;
     pu.z = uz;
 #else
-    amrex::ignore_unused(pu, eb_bnd_normal_arr, i, j, k);
+    amrex::ignore_unused(pu, AMREX_D_DECL(n0,n1,n2));
 #endif
 }
 #endif //AMREX_USE_EB

Source/Particles/PhysicalParticleContainer.cpp

@@ -1351,16 +1351,11 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
 #ifdef AMREX_USE_EB
     bool const inject_from_eb = plasma_injector.m_inject_from_eb; // whether to inject from EB or from a plane
     // Extract data structures for embedded boundaries
+    amrex::EBFArrayBoxFactory const* eb_factory = nullptr;
     amrex::FabArray<amrex::EBCellFlagFab> const* eb_flag = nullptr;
-    amrex::MultiCutFab const* eb_bnd_area = nullptr;
-    amrex::MultiCutFab const* eb_bnd_normal = nullptr;
-    amrex::MultiCutFab const* eb_bnd_cent = nullptr;
     if (inject_from_eb) {
-        amrex::EBFArrayBoxFactory const& eb_box_factory = WarpX::GetInstance().fieldEBFactory(0);
-        eb_flag = &eb_box_factory.getMultiEBCellFlagFab();
-        eb_bnd_area = &eb_box_factory.getBndryArea();
-        eb_bnd_normal = &eb_box_factory.getBndryNormal();
-        eb_bnd_cent = &eb_box_factory.getBndryCent();
+        eb_factory = &(WarpX::GetInstance().fieldEBFactory(0));
+        eb_flag = &(eb_factory->getMultiEBCellFlagFab());
     }
 #endif
 
@@ -1456,17 +1451,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
         }
 
 #ifdef AMREX_USE_EB
-        // Extract data structures for embedded boundaries
-        amrex::Array4<const typename FabArray<EBCellFlagFab>::value_type> eb_flag_arr;
-        amrex::Array4<const amrex::Real> eb_bnd_area_arr;
-        amrex::Array4<const amrex::Real> eb_bnd_normal_arr;
-        amrex::Array4<const amrex::Real> eb_bnd_cent_arr;
-        if (inject_from_eb) {
-            eb_flag_arr = eb_flag->array(mfi);
-            eb_bnd_area_arr = eb_bnd_area->array(mfi);
-            eb_bnd_normal_arr = eb_bnd_normal->array(mfi);
-            eb_bnd_cent_arr = eb_bnd_cent->array(mfi);
-        }
+        auto eb_flag_arr = eb_flag ? eb_flag->const_array(mfi) : Array4<EBCellFlag const>{};
+        auto eb_data = eb_factory ? eb_factory->getEBData(mfi) : EBData{};
 #endif
 
         amrex::ParallelForRNG(overlap_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept
@@ -1482,7 +1468,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 // Skip cells that are not partially covered by the EB
                 if (eb_flag_arr(i,j,k).isRegular() || eb_flag_arr(i,j,k).isCovered()) { return; }
                 // Scale by the (normalized) area of the EB surface in this cell
-                num_ppc_real_in_this_cell *= eb_bnd_area_arr(i,j,k);
+                num_ppc_real_in_this_cell *= eb_data.get<amrex::EBData_t::bndryarea>(i,j,k);
             }
 #else
             amrex::Real const num_ppc_real_in_this_cell = num_ppc_real; // user input: number of macroparticles per cell
@@ -1574,7 +1560,10 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
             Real scale_fac;
 #ifdef AMREX_USE_EB
             if (inject_from_eb) {
-                scale_fac = compute_scale_fac_area_eb(dx, num_ppc_real, eb_bnd_normal_arr, i, j, k );
+                scale_fac = compute_scale_fac_area_eb(dx, num_ppc_real,
+                                                      AMREX_D_DECL(eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,0),
+                                                                   eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,1),
+                                                                   eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,2)));
             } else
 #endif
             {
@@ -1595,14 +1584,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 XDim3 r;
 #ifdef AMREX_USE_EB
                 if (inject_from_eb) {
+                    auto const& pt = eb_data.randomPointOnEB(i,j,k,engine);
 #if defined(WARPX_DIM_3D)
-                    pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + eb_bnd_cent_arr(i,j,k,0))*dx[0];
-                    pos.y = overlap_corner[1] + (iv[1] + 0.5_rt + eb_bnd_cent_arr(i,j,k,1))*dx[1];
-                    pos.z = overlap_corner[2] + (iv[2] + 0.5_rt + eb_bnd_cent_arr(i,j,k,2))*dx[2];
+                    pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + pt[0])*dx[0];
+                    pos.y = overlap_corner[1] + (iv[1] + 0.5_rt + pt[1])*dx[1];
+                    pos.z = overlap_corner[2] + (iv[2] + 0.5_rt + pt[2])*dx[2];
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                    pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + eb_bnd_cent_arr(i,j,k,0))*dx[0];
+                    pos.x = overlap_corner[0] + (iv[0] + 0.5_rt + pt[0])*dx[0];
                     pos.y = 0.0_rt;
-                    pos.z = overlap_corner[1] + (iv[1] + 0.5_rt + eb_bnd_cent_arr(i,j,k,1))*dx[1];
+                    pos.z = overlap_corner[1] + (iv[1] + 0.5_rt + pt[1])*dx[1];
 #endif
                 } else
 #endif
@@ -1661,7 +1651,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                     // Injection from EB: rotate momentum according to the normal of the EB surface
                     // (The above code initialized the momentum by assuming that z is the direction
                     // normal to the EB surface. Thus we need to rotate from z to the normal.)
-                    rotate_momentum_eb(pu, eb_bnd_normal_arr, i, j , k);
+                    rotate_momentum_eb(pu, AMREX_D_DECL(eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,0),
+                                                        eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,1),
+                                                        eb_data.get<amrex::EBData_t::bndrynorm>(i,j,k,2)));
                 }
 #endif