From 3c6a88d6e613a5d9b361ad8eb672256d5fb5a324 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 2 Dec 2024 15:29:50 -0800 Subject: [PATCH 1/2] Flux injection from EB: Pick a random point instead of the centroid --- Source/Particles/AddPlasmaUtilities.H | 39 ++++++++-------- .../Particles/PhysicalParticleContainer.cpp | 46 ++++++++----------- 2 files changed, 37 insertions(+), 48 deletions(-) diff --git a/Source/Particles/AddPlasmaUtilities.H b/Source/Particles/AddPlasmaUtilities.H index bb05d7be3c8..6f0be2d62c4 100644 --- a/Source/Particles/AddPlasmaUtilities.H +++ b/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& dx, const amrex::Real num_ppc_real, - amrex::Array4 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,8 +155,9 @@ amrex::Real compute_scale_fac_area_eb ( AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void rotate_momentum_eb ( PDim3 & pu, - amrex::Array4 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; @@ -168,9 +165,9 @@ 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 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 diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c973e9afafa..348e0c95fcf 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/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 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::value_type> eb_flag_arr; - amrex::Array4 eb_bnd_area_arr; - amrex::Array4 eb_bnd_normal_arr; - amrex::Array4 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 const& eb_flag_arr = eb_flag->const_array(mfi); + auto const& eb_data = eb_factory->getEBData(mfi); #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(i,j,k); } else #else amrex::Real const num_ppc_real_in_this_cell = num_ppc_real; // user input: number of macroparticles per cell @@ -1577,7 +1563,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(i,j,k,0), + eb_data.get(i,j,k,1), + eb_data.get(i,j,k,2))); } else #endif { @@ -1598,14 +1587,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 @@ -1664,7 +1654,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(i,j,k,0), + eb_data.get(i,j,k,1), + eb_data.get(i,j,k,2))); } #endif From 59ef806e3c6ad9c4898df6cc0cb0d8752822ca26 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 16 Dec 2024 10:06:03 -0800 Subject: [PATCH 2/2] Fix non-EB --- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 348e0c95fcf..afe625efbe4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1451,8 +1451,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef AMREX_USE_EB - auto const& eb_flag_arr = eb_flag->const_array(mfi); - auto const& eb_data = eb_factory->getEBData(mfi); + auto eb_flag_arr = eb_flag ? eb_flag->const_array(mfi) : Array4{}; + 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