Skip to content

Commit

Permalink
smooth tanh weighting for fp and cp interpolation in buffer region
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Mar 18, 2024
1 parent 6be87d6 commit 3059022
Show file tree
Hide file tree
Showing 5 changed files with 275 additions and 14 deletions.
1 change: 1 addition & 0 deletions Source/Parallelization/WarpXRegrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,6 +334,7 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi
// we can avoid redistributing these since we immediately re-build the values via BuildBufferMasks()
RemakeMultiFab(current_buffer_masks[lev], dm, false ,lev);
RemakeMultiFab(gather_buffer_masks[lev], dm, false ,lev);
RemakeMultiFab(interp_weight_gbuffer[lev], dm, false, lev);

if (current_buffer_masks[lev] || gather_buffer_masks[lev]) {
BuildBufferMasks();
Expand Down
22 changes: 22 additions & 0 deletions Source/Particles/PhysicalParticleContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,6 +326,28 @@ public:
amrex::FArrayBox const * & ez_ptr, amrex::FArrayBox const * & bx_ptr,
amrex::FArrayBox const * & by_ptr, amrex::FArrayBox const * & bz_ptr);

void InterpolateFieldsInGatherBuffer (const amrex::Box& box,
amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz,
amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz,
amrex::Elixir& buf_exeli, amrex::Elixir& buf_eyeli, amrex::Elixir& buf_ezeli,
amrex::Elixir& buf_bxeli, amrex::Elixir& buf_byeli, amrex::Elixir& buf_bzeli,
amrex::FArrayBox const *& ex_ptr, amrex::FArrayBox const *& ey_ptr, amrex::FArrayBox const* & ez_ptr,
amrex::FArrayBox const *& bx_ptr, amrex::FArrayBox const *& by_ptr, amrex::FArrayBox const* & bz_ptr,
amrex::FArrayBox const *& cex_ptr, amrex::FArrayBox const *& cey_ptr, amrex::FArrayBox const* & cez_ptr,
amrex::FArrayBox const *& cbx_ptr, amrex::FArrayBox const *& cby_ptr, amrex::FArrayBox const* & cbz_ptr,
const amrex::IArrayBox& gather_mask, const amrex::FArrayBox& weight_gbuffer,
const amrex::IntVect ref_ratio);

void WeightedSumInGatherBuffer (const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox,
amrex::Array4<amrex::Real> const & xbuf_field_arr, amrex::Array4<amrex::Real> const& ybuf_field_arr,
amrex::Array4<amrex::Real> const & zbuf_field_arr,
amrex::Array4<const amrex::Real> cfx_arr, amrex::Array4<const amrex::Real> cfy_arr,
amrex::Array4<const amrex::Real> cfz_arr,
amrex::Array4<const amrex::Real> fx_arr, amrex::Array4<const amrex::Real> fy_arr,
amrex::Array4<const amrex::Real> fz_arr,
amrex::Array4<const int > gm_arr, amrex::Array4<const amrex::Real> wt_arr,
const amrex::IntVect ref_ratio);

/**
* \brief This function determines if resampling should be done for the current species, and
* if so, performs the resampling.
Expand Down
142 changes: 131 additions & 11 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2016,12 +2016,11 @@ PhysicalParticleContainer::Evolve (int lev,
WARPX_PROFILE_VAR_NS("PhysicalParticleContainer::Evolve::GatherAndPush", blp_fg);

BL_ASSERT(OnSameGrids(lev,jx));

amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);

const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev);
const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev);

const amrex::MultiFab* weight_gbuffer = WarpX::GetInstance().getInterpWeightGBuffer(lev);
const bool has_buffer = cEx || cjx;

if (m_do_back_transformed_particles)
Expand Down Expand Up @@ -2050,7 +2049,8 @@ PhysicalParticleContainer::Evolve (int lev,

FArrayBox filtered_Ex, filtered_Ey, filtered_Ez;
FArrayBox filtered_Bx, filtered_By, filtered_Bz;

FArrayBox bufferEx, bufferEy, bufferEz;
FArrayBox bufferBx, bufferBy, bufferBz;
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
Expand Down Expand Up @@ -2079,6 +2079,7 @@ PhysicalParticleContainer::Evolve (int lev,
FArrayBox const* bzfab = &Bz[pti];

Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli;
Elixir buf_exeli, buf_eyeli, buf_ezeli, buf_bxeli, buf_byeli, buf_bzeli;

if (WarpX::use_fdtd_nci_corr)
{
Expand Down Expand Up @@ -2120,7 +2121,7 @@ PhysicalParticleContainer::Evolve (int lev,

DepositCharge(pti, wp, ion_lev, rho, 0, 0,
np_current, thread_num, lev, lev);
if (has_buffer){
if ((np-np_current)> 0 ){
DepositCharge(pti, wp, ion_lev, crho, 0, np_current,
np-np_current, thread_num, lev, lev-1);
}
Expand Down Expand Up @@ -2163,8 +2164,19 @@ PhysicalParticleContainer::Evolve (int lev,
FArrayBox const* cbyfab = &(*cBy)[pti];
FArrayBox const* cbzfab = &(*cBz)[pti];

//// buffer box (Both 2D and 3D)
InterpolateFieldsInGatherBuffer(box, bufferEx, bufferEy, bufferEz,
bufferBx, bufferBy, bufferBz,
buf_exeli, buf_eyeli, buf_ezeli,
buf_bxeli, buf_byeli, buf_bzeli,
exfab, eyfab, ezfab, bxfab, byfab, bzfab,
cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab,
(*gather_masks)[pti], (*weight_gbuffer)[pti],
ref_ratio);

if (WarpX::use_fdtd_nci_corr)
{
// should this be bufferEFields*
// Filter arrays (*cEx)[pti], store the result in
// filtered_Ex and update pointer cexfab so that it
// points to filtered_Ex (and do the same for all
Expand All @@ -2176,15 +2188,17 @@ PhysicalParticleContainer::Evolve (int lev,
(*cBx)[pti], (*cBy)[pti], (*cBz)[pti],
cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab);
}

// Field gather and push for particles in gather buffers
e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
if (push_type == PushType::Explicit) {
PushPX(pti, cexfab, ceyfab, cezfab,
cbxfab, cbyfab, cbzfab,
//PushPX(pti, cexfab, ceyfab, cezfab,
// cbxfab, cbyfab, cbzfab,
PushPX(pti, &bufferEx, &bufferEy, &bufferEz,
&bufferBx, &bufferBy, &bufferBz,
cEx->nGrowVect(), e_is_nodal,
nfine_gather, np-nfine_gather,
lev, lev-1, dt, ScaleFields(false), a_dt_type);
lev, gather_lev, dt, ScaleFields(false), a_dt_type);
//lev, lev-1, dt, ScaleFields(false), a_dt_type);
} else if (push_type == PushType::Implicit) {
ImplicitPushXP(pti, cexfab, ceyfab, cezfab,
cbxfab, cbyfab, cbzfab,
Expand All @@ -2193,7 +2207,6 @@ PhysicalParticleContainer::Evolve (int lev,
lev, lev-1, dt, ScaleFields(false), a_dt_type);
}
}

WARPX_PROFILE_VAR_STOP(blp_fg);

// Current Deposition
Expand All @@ -2210,7 +2223,7 @@ PhysicalParticleContainer::Evolve (int lev,
0, np_current, thread_num,
lev, lev, dt, relative_time, push_type);

if (has_buffer)
if ((np-np_current)>0)
{
// Deposit in buffers
DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz,
Expand All @@ -2230,7 +2243,7 @@ PhysicalParticleContainer::Evolve (int lev,

DepositCharge(pti, wp, ion_lev, rho, 1, 0,
np_current, thread_num, lev, lev);
if (has_buffer){
if ((np-np_current)>0 ){
DepositCharge(pti, wp, ion_lev, crho, 1, np_current,
np-np_current, thread_num, lev, lev-1);
}
Expand Down Expand Up @@ -2333,6 +2346,113 @@ PhysicalParticleContainer::applyNCIFilter (
#endif
}

void
PhysicalParticleContainer::InterpolateFieldsInGatherBuffer (const Box& box,
amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz,
amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz,
Elixir& buf_exeli, Elixir& buf_eyeli, Elixir& buf_ezeli,
Elixir& buf_bxeli, Elixir& buf_byeli, Elixir& buf_bzeli,
FArrayBox const *& ex_ptr, FArrayBox const *& ey_ptr, FArrayBox const* & ez_ptr,
FArrayBox const *& bx_ptr, FArrayBox const *& by_ptr, FArrayBox const* & bz_ptr,
FArrayBox const *& cex_ptr, FArrayBox const *& cey_ptr, FArrayBox const* & cez_ptr,
FArrayBox const *& cbx_ptr, FArrayBox const *& cby_ptr, FArrayBox const* & cbz_ptr,
const IArrayBox& gather_mask, const FArrayBox& weight_gbuffer,
const amrex::IntVect ref_ratio)
{
amrex::Box tmp_exbox = amrex::convert(box,ex_ptr->box().ixType());
amrex::Box tmp_eybox = amrex::convert(box,ey_ptr->box().ixType());
amrex::Box tmp_ezbox = amrex::convert(box,ez_ptr->box().ixType());

bufferEx.resize(tmp_exbox);
buf_exeli = bufferEx.elixir();
bufferEy.resize(tmp_eybox);
buf_eyeli = bufferEy.elixir();
bufferEz.resize(tmp_ezbox);
buf_ezeli = bufferEz.elixir();
WeightedSumInGatherBuffer(tmp_exbox, tmp_eybox, tmp_ezbox,
bufferEx.array(), bufferEy.array(), bufferEz.array(),
cex_ptr->array(), cey_ptr->array(), cez_ptr->array(),
ex_ptr->array(), ey_ptr->array(), ez_ptr->array(),
gather_mask.array(), weight_gbuffer.array(), ref_ratio);

amrex::Box tmp_bxbox = amrex::convert(box,bx_ptr->box().ixType());
amrex::Box tmp_bybox = amrex::convert(box,by_ptr->box().ixType());
amrex::Box tmp_bzbox = amrex::convert(box,bz_ptr->box().ixType());
bufferBx.resize(tmp_bxbox);
buf_bxeli = bufferBx.elixir();
bufferBy.resize(tmp_bybox);
buf_byeli = bufferBy.elixir();
bufferBz.resize(tmp_bzbox);
buf_bzeli = bufferBz.elixir();

WeightedSumInGatherBuffer(tmp_bxbox, tmp_bybox, tmp_bzbox,
bufferBx.array(), bufferBy.array(), bufferBz.array(),
cbx_ptr->array(), cby_ptr->array(), cbz_ptr->array(),
bx_ptr->array(), by_ptr->array(), bz_ptr->array(),
gather_mask.array(), weight_gbuffer.array(), ref_ratio);
}

void
PhysicalParticleContainer::WeightedSumInGatherBuffer (
const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox,
amrex::Array4<amrex::Real> const & xbuf_field_arr, amrex::Array4<amrex::Real> const& ybuf_field_arr,
amrex::Array4<amrex::Real> const & zbuf_field_arr,
amrex::Array4<const amrex::Real> cfx_arr, amrex::Array4<const amrex::Real> cfy_arr,
amrex::Array4<const amrex::Real> cfz_arr,
amrex::Array4<const amrex::Real> fx_arr, amrex::Array4<const amrex::Real> fy_arr,
amrex::Array4<const amrex::Real> fz_arr,
amrex::Array4<const int > gm_arr, amrex::Array4<const amrex::Real> wt_arr,
const amrex::IntVect ref_ratio)
{
amrex::ParallelFor( tmp_xbox, tmp_ybox, tmp_zbox,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
int ii = amrex::coarsen(i,ref_ratio[0]);
int jj = j;
int kk = k;
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
jj = amrex::coarsen(j,ref_ratio[1]);
#elif defined(WARPX_DIM_3D)
kk = amrex::coarsen(k,ref_ratio[2]);
#endif
if (gm_arr(i,j,k) == 0) {
xbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fx_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfx_arr(ii,jj,kk);
} else {
xbuf_field_arr(i,j,k) = fx_arr(i,j,k);
}
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
int ii = amrex::coarsen(i,ref_ratio[0]);
int jj = j;
int kk = k;
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
jj = amrex::coarsen(j,ref_ratio[1]);
#elif defined(WARPX_DIM_3D)
kk = amrex::coarsen(k,ref_ratio[2]);
#endif
if (gm_arr(i,j,k) == 0) {
ybuf_field_arr(i,j,k) = wt_arr(i,j,k)*fy_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfy_arr(ii,jj,kk);
} else {
ybuf_field_arr(i,j,k) = fy_arr(i,j,k);
}
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
int ii = amrex::coarsen(i,ref_ratio[0]);
int jj = j;
int kk = k;
#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
jj = amrex::coarsen(j,ref_ratio[1]);
#elif defined(WARPX_DIM_3D)
kk = amrex::coarsen(k,ref_ratio[2]);
#endif
if (gm_arr(i,j,k) == 0) {
zbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fz_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfz_arr(ii,jj,kk);
} else {
zbuf_field_arr(i,j,k) = fz_arr(i,j,k);
}
}
);
}

// Loop over all particles in the particle container and
// split particles tagged with p.id()=DoSplitParticleID
void
Expand Down
36 changes: 36 additions & 0 deletions Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,6 +355,32 @@ public:
//! #n_field_gather_buffer cells of the edge of the patch, will gather the fields
//! from the lower refinement level instead of the refinement patch itself
static int n_field_gather_buffer;

/**
* If true, particles located inside refinement patch but within
* #n_field_gather_buffer cells of the patch edge will gather fields
* that are a weighted sum of the fields from lower refinement
* level and the refinement patch. By default, no interpolation will be done,
* and the fields will be gathered only from the lower refinement level.
* The option can be used to allow for a smooth transition from the fine
* and coarse solutions with a tanh profile, such that, the value
* of the gather fields is close to the solution on the refinement patch
* near the gather buffer edge, and close to the solution on the lower
* refinement level close to the edge of the patch.
* #tanh_midpoint_gather_buffer can be used to set the tanh profile
*/
static bool do_fieldinterp_gather_buffer;

/**
* With mesh-refinement, finite number of gather buffer cells, and
* do_fieldinterp_gather_buffer, the particles on the fine patch
* gather a weighted sum of fields from fine patch and coarse patch.
* A tanh profile is used to set the smoothing kernel
* #tanh_midpoint_gather_buffer must be a fraction of the buffer width
* measured from coarsefine interface, where the tanh profile is 0.5
* The default value is 0.5
*/
static amrex::Real tanh_midpoint_gather_buffer;
//! With mesh refinement, particles located inside a refinement patch, but within
//! #n_current_deposition_buffer cells of the edge of the patch, will deposit their charge
//! and current onto the lower refinement level instead of the refinement patch itself
Expand Down Expand Up @@ -878,6 +904,10 @@ public:
static std::array<amrex::Real,3> CellSize (int lev);
static amrex::RealBox getRealBox(const amrex::Box& bx, int lev);

[[nodiscard]] const amrex::MultiFab* getInterpWeightGBuffer (int lev) const
{
return interp_weight_gbuffer[lev].get();
}
/**
* \brief Return the lower corner of the box in real units.
* \param bx The box
Expand Down Expand Up @@ -1391,6 +1421,7 @@ private:
return gather_buffer_masks[lev].get();
}


/**
* \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for
* finite-order centering operations. For example, for finite-order centering of order 6,
Expand Down Expand Up @@ -1564,6 +1595,11 @@ private:
amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cax;
amrex::Vector<std::unique_ptr<amrex::iMultiFab> > current_buffer_masks;
amrex::Vector<std::unique_ptr<amrex::iMultiFab> > gather_buffer_masks;
/** Multifab that stores weights for interpolating fine and coarse solutions
* in the buffer gather region, such that, the weight for fine solution is
* 0 near the coarse-fine interface, and 1 as the buffer gather region terminates in the fine patch
*/
amrex::Vector<std::unique_ptr<amrex::MultiFab> > interp_weight_gbuffer;

// If charge/current deposition buffers are used
amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_buf;
Expand Down
Loading

0 comments on commit 3059022

Please sign in to comment.