Skip to content

Commit

Permalink
Flux injection from EB: Pick a random point instead of the centroid
Browse files Browse the repository at this point in the history
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Dec 3, 2024
1 parent 780a584 commit 40e38e7
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 48 deletions.
39 changes: 18 additions & 21 deletions Source/Particles/AddPlasmaUtilities.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
Expand All @@ -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
Expand Down
46 changes: 19 additions & 27 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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 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
Expand All @@ -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
#else
amrex::Real const num_ppc_real_in_this_cell = num_ppc_real; // user input: number of macroparticles per cell
Expand Down Expand Up @@ -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<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
{
Expand All @@ -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
Expand Down Expand Up @@ -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<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

Expand Down

0 comments on commit 40e38e7

Please sign in to comment.