Skip to content

Commit

Permalink
remove buffer each dir and initializeexternal field on Bfield_fp
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Aug 13, 2024
1 parent eea5eb8 commit 0c42c38
Show file tree
Hide file tree
Showing 5 changed files with 100 additions and 74 deletions.
120 changes: 74 additions & 46 deletions Source/Initialization/WarpXInitData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -866,12 +866,12 @@ WarpX::InitLevelData (int lev, Real /*time*/)
// Externally imposed fields are only initialized until the user-defined maxlevel_extEMfield_init.
// The default maxlevel_extEMfield_init value is the total number of levels in the simulation
if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
&& (lev > 0) && (lev <= maxlevel_extEMfield_init)) {
&& (lev <= maxlevel_extEMfield_init)) {

InitializeExternalFieldsOnGridUsingParser(
Bfield_aux[lev][0].get(),
Bfield_aux[lev][1].get(),
Bfield_aux[lev][2].get(),
Bfield_fp[lev][0].get(),
Bfield_fp[lev][1].get(),
Bfield_fp[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
Expand All @@ -880,17 +880,31 @@ WarpX::InitLevelData (int lev, Real /*time*/)
'B',
lev, PatchType::fine);

InitializeExternalFieldsOnGridUsingParser(
Bfield_cp[lev][0].get(),
Bfield_cp[lev][1].get(),
Bfield_cp[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::coarse);
if (lev > 0) {
InitializeExternalFieldsOnGridUsingParser(
Bfield_aux[lev][0].get(),
Bfield_aux[lev][1].get(),
Bfield_aux[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::fine);

InitializeExternalFieldsOnGridUsingParser(
Bfield_cp[lev][0].get(),
Bfield_cp[lev][1].get(),
Bfield_cp[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::coarse);
}
}
ParmParse pp_warpx("warpx");
int IncludeBfieldPerturbation = 0;
Expand Down Expand Up @@ -939,6 +953,18 @@ WarpX::InitLevelData (int lev, Real /*time*/)
if ((m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function)
&& (lev <= maxlevel_extEMfield_init)) {

InitializeExternalFieldsOnGridUsingParser(
Efield_fp[lev][0].get(),
Efield_fp[lev][1].get(),
Efield_fp[lev][2].get(),
m_p_ext_field_params->Exfield_parser->compile<3>(),
m_p_ext_field_params->Eyfield_parser->compile<3>(),
m_p_ext_field_params->Ezfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'E',
lev, PatchType::fine);

#ifdef AMREX_USE_EB
// We initialize ECTRhofield consistently with the Efield
if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) {
Expand Down Expand Up @@ -1392,22 +1418,23 @@ WarpX::LoadExternalFields (int const lev)
"External fields from file are not compatible with the moving window." );
}

// External grid fields
if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) {
// Initialize Bfield_fp_external with external function
InitializeExternalFieldsOnGridUsingParser(
Bfield_fp_external[lev][0].get(),
Bfield_fp_external[lev][1].get(),
Bfield_fp_external[lev][2].get(),
m_p_ext_field_params->Bxfield_parser->compile<3>(),
m_p_ext_field_params->Byfield_parser->compile<3>(),
m_p_ext_field_params->Bzfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'B',
lev, PatchType::fine);
}
else if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) {
//// External grid fields
//if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) {
// // Initialize Bfield_fp_external with external function
// InitializeExternalFieldsOnGridUsingParser(
// Bfield_fp_external[lev][0].get(),
// Bfield_fp_external[lev][1].get(),
// Bfield_fp_external[lev][2].get(),
// m_p_ext_field_params->Bxfield_parser->compile<3>(),
// m_p_ext_field_params->Byfield_parser->compile<3>(),
// m_p_ext_field_params->Bzfield_parser->compile<3>(),
// m_edge_lengths[lev],
// m_face_areas[lev],
// 'B',
// lev, PatchType::fine);
//}
//else
if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) {
#if defined(WARPX_DIM_RZ)
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1,
"External field reading is not implemented for more than one RZ mode (see #3829)");
Expand All @@ -1421,21 +1448,22 @@ WarpX::LoadExternalFields (int const lev)
#endif
}

if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function) {
// Initialize Efield_fp_external with external function
InitializeExternalFieldsOnGridUsingParser(
Efield_fp_external[lev][0].get(),
Efield_fp_external[lev][1].get(),
Efield_fp_external[lev][2].get(),
m_p_ext_field_params->Exfield_parser->compile<3>(),
m_p_ext_field_params->Eyfield_parser->compile<3>(),
m_p_ext_field_params->Ezfield_parser->compile<3>(),
m_edge_lengths[lev],
m_face_areas[lev],
'E',
lev, PatchType::fine);
}
else if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) {
//if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::parse_ext_grid_function) {
// // Initialize Efield_fp_external with external function
// InitializeExternalFieldsOnGridUsingParser(
// Efield_fp_external[lev][0].get(),
// Efield_fp_external[lev][1].get(),
// Efield_fp_external[lev][2].get(),
// m_p_ext_field_params->Exfield_parser->compile<3>(),
// m_p_ext_field_params->Eyfield_parser->compile<3>(),
// m_p_ext_field_params->Ezfield_parser->compile<3>(),
// m_edge_lengths[lev],
// m_face_areas[lev],
// 'E',
// lev, PatchType::fine);
//}
//else
if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) {
#if defined(WARPX_DIM_RZ)
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1,
"External field reading is not implemented for more than one RZ mode (see #3829)");
Expand Down
2 changes: 1 addition & 1 deletion Source/Parallelization/WarpXRegrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -465,7 +465,7 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi
#endif
}

if (lev > 0 && (n_field_gather_buffer_each_dir.max() > 0 || n_current_deposition_buffer_each_dir.max() > 0)) {
if (lev > 0 && (n_field_gather_buffer > 0 || n_current_deposition_buffer > 0)) {
for (int idim=0; idim < 3; ++idim)
{
RemakeMultiFab(Bfield_cax[lev][idim], false);
Expand Down
8 changes: 4 additions & 4 deletions Source/Particles/Sorting/Partition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers(

// - Select the larger buffer
iMultiFab const* bmasks =
(WarpX::n_field_gather_buffer_each_dir.max() >= WarpX::n_current_deposition_buffer_each_dir.max()) ?
(WarpX::n_field_gather_buffer >= WarpX::n_current_deposition_buffer) ?
gather_masks : current_masks;
// - For each particle, find whether it is in the larger buffer,
// by looking up the mask. Store the answer in `inexflag`.
Expand All @@ -86,7 +86,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers(
// Second, among particles that are in the larger buffer, partition
// particles into the smaller buffer

if (WarpX::n_current_deposition_buffer_each_dir.max() == WarpX::n_field_gather_buffer_each_dir.max()) {
if (WarpX::n_current_deposition_buffer == WarpX::n_field_gather_buffer ) {
// No need to do anything if the buffers have the same size
nfine_current = nfine_gather = iteratorDistance(pid.begin(), sep);
} else if (sep == pid.end()) {
Expand All @@ -97,11 +97,11 @@ PhysicalParticleContainer::PartitionParticlesInBuffers(
if (bmasks == gather_masks) {
nfine_gather = n_fine;
bmasks = current_masks;
n_buf = WarpX::n_current_deposition_buffer_each_dir.max();
n_buf = WarpX::n_current_deposition_buffer;
} else {
nfine_current = n_fine;
bmasks = gather_masks;
n_buf = WarpX::n_field_gather_buffer_each_dir.max();
n_buf = WarpX::n_field_gather_buffer;
}
if (n_buf > 0)
{
Expand Down
20 changes: 10 additions & 10 deletions Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -914,10 +914,10 @@ public:
static amrex::XDim3 InvCellSize (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();
}
// [[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 @@ -1126,7 +1126,7 @@ public:

// for cuda
void BuildBufferMasksInBox ( amrex::Box tbx, amrex::IArrayBox &buffer_mask,
const amrex::IArrayBox &guard_mask, const amrex::IntVect ng );
const amrex::IArrayBox &guard_mask, const int ng );

#ifdef AMREX_USE_EB
amrex::EBFArrayBoxFactory const& fieldEBFactory (int lev) const noexcept {
Expand Down Expand Up @@ -1602,11 +1602,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;
// /** 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
24 changes: 11 additions & 13 deletions Source/WarpX.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -402,7 +402,7 @@ WarpX::WarpX ()
gather_buffer_masks.resize(nlevs_max);
current_buf.resize(nlevs_max);
charge_buf.resize(nlevs_max);
interp_weight_gbuffer.resize(nlevs_max);
// interp_weight_gbuffer.resize(nlevs_max);

pml.resize(nlevs_max);
#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT)
Expand Down Expand Up @@ -2160,7 +2160,7 @@ WarpX::ClearLevel (int lev)

current_buffer_masks[lev].reset();
gather_buffer_masks[lev].reset();
interp_weight_gbuffer[lev].reset();
// interp_weight_gbuffer[lev].reset();

F_fp [lev].reset();
G_fp [lev].reset();
Expand Down Expand Up @@ -2253,7 +2253,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
}
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
if (n_field_gather_buffer_each_dir[i] < 0) {
n_field_gather_buffer_each_dir[i] = n_current_deposition_buffer + 1;
n_field_gather_buffer_each_dir[i] = n_current_deposition_buffer_each_dir[i] + 1;
}
}
AllocLevelMFs(lev, ba, dm, guard_cells.ng_alloc_EB, guard_cells.ng_alloc_J,
Expand Down Expand Up @@ -2799,14 +2799,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
//
// Copy of the coarse aux
//
if (lev > 0 && (n_field_gather_buffer_each_dir.max() > 0 ||
n_current_deposition_buffer_each_dir.max() > 0 ||
if (lev > 0 && (n_field_gather_buffer > 0 ||
n_current_deposition_buffer > 0 ||
mypc->nSpeciesGatherFromMainGrid() > 0))
{
BoxArray cba = ba;
cba.coarsen(refRatio(lev-1));

if (n_field_gather_buffer_each_dir.max() > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) {
if (n_field_gather_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) {
if (aux_is_nodal) {
BoxArray const& cnba = amrex::convert(cba,IntVect::TheNodeVector());
AllocInitMultiFab(Bfield_cax[lev][0], cnba,dm,ncomps,ngEB,lev, "Bfield_cax[x]", 0.0_rt);
Expand All @@ -2830,10 +2830,10 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
// Gather buffer masks have 1 ghost cell, because of the fact
// that particles may move by more than one cell when using subcycling.
AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks");
AllocInitMultiFab(interp_weight_gbuffer[lev], amrex::convert(ba, IntVect::TheNodeVector()), dm, ncomps, ngEB, lev, "interp_weight_gbuffer", 0.0_rt);
// AllocInitMultiFab(interp_weight_gbuffer[lev], amrex::convert(ba, IntVect::TheNodeVector()), dm, ncomps, ngEB, lev, "interp_weight_gbuffer", 0.0_rt);
}

if (n_current_deposition_buffer_each_dir.max() > 0) {
if (n_current_deposition_buffer > 0) {
amrex::Print() << " current dep buffer : " << n_current_deposition_buffer << "\n";
AllocInitMultiFab(current_buf[lev][0], amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[x]");
AllocInitMultiFab(current_buf[lev][1], amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[y]");
Expand Down Expand Up @@ -3215,9 +3215,7 @@ WarpX::BuildBufferMasks ()
{
for (int ipass = 0; ipass < 2; ++ipass)
{
//const int ngbuffer = (ipass == 0) ? n_current_deposition_buffer : n_field_gather_buffer;
const amrex::IntVect ngbuffer = (ipass == 0) ? n_current_deposition_buffer_each_dir
: n_field_gather_buffer_each_dir;
const int ngbuffer = (ipass == 0) ? n_current_deposition_buffer : n_field_gather_buffer;
iMultiFab* bmasks = (ipass == 0) ? current_buffer_masks[lev].get() : gather_buffer_masks[lev].get();
if (bmasks)
{
Expand Down Expand Up @@ -3253,11 +3251,11 @@ WarpX::BuildBufferMasks ()
*/
void
WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_mask,
const amrex::IArrayBox &guard_mask, const amrex::IntVect ng )
const amrex::IArrayBox &guard_mask, const int ng )
{
auto const& msk = buffer_mask.array();
auto const& gmsk = guard_mask.const_array();
const amrex::Dim3 ng3 = ng.dim3();
const amrex::Dim3 ng3 = amrex::IntVect(ng).dim3();
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
for (int kk = k-ng3.z; kk <= k+ng3.z; ++kk) {
Expand Down

0 comments on commit 0c42c38

Please sign in to comment.