Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 2 commits into
base: development
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3c6a88d
Flux injection from EB: Pick a random point instead of the centroid
WeiqunZhang Dec 2, 2024
59ef806
Fix non-EB
WeiqunZhang Dec 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 18 additions & 21 deletions Source/Particles/AddPlasmaUtilities.H
View file
Open in desktop
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
Open in desktop
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 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
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
Loading