Skip to content

Commit

Permalink
New embedded boundary data structures in preparation for #5534 (#5574)
Browse files Browse the repository at this point in the history 
When updating `E` (and when initializing the field values for `E` and
`B` with a parser), we used to check somewhat complicated conditions
that relied on the MultiFabs `edge_lengths` and `face_areas`.

This PR introduces separate arrays (`m_eb_update_E` and `m_eb_update_B`,
which are `iMultiFab` instead of `MultiFab` and thus take up much less
memory). These arrays contain flags that keep track of where to udpate
the fields (i.e. the black crosses in the above figure). These arrays
are initialized in separate functions which keep the same EB definition.
(PR #5534 will then change this definition.) The code for the field
pusher and field initialization uses these arrays (instead of
edge_lengths and face_areas. It is thus significantly simpler and avoids
duplicating complicated if conditions.

The above changes have not yet been implemented for the hybrid solver.
This will instead be done in a separate PR:
#5558. Once this is complete, the
MultiFab `edge_lengths` and `face_areas` will not be needed anymore
(except for the ECT solver), and we will thus only allocate them for the
ECT solver. This should save a significant amount of memory.

---------

Co-authored-by: Roelof Groenewald <[email protected]>
  • Loading branch information
@RemiLehe @roelof-groenewald
RemiLehe and roelof-groenewald authored Jan 21, 2025
1 parent 54333c7 commit e3378b0
Show file tree
Hide file tree
Showing 12 changed files with 409 additions and 185 deletions.
4 changes: 2 additions & 2 deletions Examples/Tests/embedded_boundary_cube/inputs_base_3d
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ boundary.field_lo = pec pec pec
boundary.field_hi = pec pec pec

eb2.geom_type = box
eb2.box_lo = -0.5 -0.5 -0.5
eb2.box_hi = 0.5 0.5 0.5
eb2.box_lo = -0.5 -0.5 -0.5 # Ensures that the stair-case EB is exactly at -0.5
eb2.box_hi = 0.5 0.5 0.5 # Ensures that the stair-case EB is exactly at 0.5
eb2.box_has_fluid_inside = true
# Alternatively one could use parser to build EB
# Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.
Expand Down
219 changes: 217 additions & 2 deletions Source/EmbeddedBoundary/WarpXInitEB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,9 +291,224 @@ WarpX::ScaleAreas (ablastr::fields::VectorField& face_areas,
}
}

void
WarpX::MarkUpdateCellsStairCase (
std::array< std::unique_ptr<amrex::iMultiFab>,3> & eb_update,
ablastr::fields::VectorField const& field,
amrex::EBFArrayBoxFactory const & eb_fact )
{

using ablastr::fields::Direction;
using warpx::fields::FieldType;

// Extract structures for embedded boundaries
amrex::FabArray<amrex::EBCellFlagFab> const& eb_flag = eb_fact.getMultiEBCellFlagFab();

for (int idim = 0; idim < 3; ++idim) {

#ifdef AMREX_USE_OMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for (amrex::MFIter mfi(*field[idim]); mfi.isValid(); ++mfi) {

const amrex::Box& box = mfi.tilebox();
amrex::Array4<int> const & eb_update_arr = eb_update[idim]->array(mfi);

// Check if the box (including one layer of guard cells) contains a mix of covered and regular cells
const amrex::Box& eb_info_box = mfi.tilebox(amrex::IntVect::TheCellVector()).grow(1);
amrex::FabType const fab_type = eb_flag[mfi].getType( eb_info_box );

if (fab_type == amrex::FabType::regular) { // All cells in the box are regular

// Every cell in box is all regular: update field in every cell
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
eb_update_arr(i, j, k) = 1;
});

} else if (fab_type == amrex::FabType::covered) { // All cells in the box are covered

// Every cell in box is all covered: do not update field
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
eb_update_arr(i, j, k) = 0;
});

} else { // The box contains a mix of covered and regular cells

auto const & flag = eb_flag[mfi].array();
auto index_type = field[idim]->ixType();

amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {

// Stair-case approximation: If neighboring cells of this gridpoint
// are either partially or fully covered: do not update field

// The number of cells that we need to check depend on the index type
// of the `eb_update_arr` in each direction.
// If `eb_update_arr` is nodal in a given direction, we need to check the cells
// to the left and right of this nodal gridpoint.
// For instance, if `eb_update_arr` is nodal in the first dimension, we need
// to check the cells at index i-1 and at index i, since, with AMReX indexing conventions,
// these are the neighboring cells for the nodal gripoint at index i.
// If `eb_update_arr` is cell-centerd in a given direction, we only need to check
// the cell at the same position (e.g., in the first dimension: the cell at index i).
int const i_start = ( index_type.nodeCentered(0) )? i-1 : i;
#if AMREX_SPACEDIM > 1
int const j_start = ( index_type.nodeCentered(1) )? j-1 : j;
#else
int const j_start = j;
#endif
#if AMREX_SPACEDIM > 2
int const k_start = ( index_type.nodeCentered(2) )? k-1 : k;
#else
int const k_start = k;
#endif
// Loop over neighboring cells
int eb_update = 1;
for (int i_cell = i_start; i_cell <= i; ++i_cell) {
for (int j_cell = j_start; j_cell <= j; ++j_cell) {
for (int k_cell = k_start; k_cell <= k; ++k_cell) {
// If one of the neighboring is either partially or fully covered
// (i.e. if they are not regular cells), do not update field
// (Note that `flag` is a cell-centered object, and `isRegular`
// returns `false` if the cell is either partially or fully covered.)
if ( !flag(i_cell, j_cell, k_cell).isRegular() ) {
eb_update = 0;
}
}
}
}
eb_update_arr(i, j, k) = eb_update;
});

}

}

}

}

void
WarpX::MarkUpdateECellsECT (
std::array< std::unique_ptr<amrex::iMultiFab>,3> & eb_update_E,
ablastr::fields::VectorField const& edge_lengths )
{

#ifdef AMREX_USE_OMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( amrex::MFIter mfi(*eb_update_E[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {

const amrex::Box& tbx = mfi.tilebox( eb_update_E[0]->ixType().toIntVect(), eb_update_E[0]->nGrowVect() );
const amrex::Box& tby = mfi.tilebox( eb_update_E[1]->ixType().toIntVect(), eb_update_E[1]->nGrowVect() );
const amrex::Box& tbz = mfi.tilebox( eb_update_E[2]->ixType().toIntVect(), eb_update_E[2]->nGrowVect() );

amrex::Array4<int> const & eb_update_Ex_arr = eb_update_E[0]->array(mfi);
amrex::Array4<int> const & eb_update_Ey_arr = eb_update_E[1]->array(mfi);
amrex::Array4<int> const & eb_update_Ez_arr = eb_update_E[2]->array(mfi);

amrex::Array4<amrex::Real> const & lx_arr = edge_lengths[0]->array(mfi);
amrex::Array4<amrex::Real> const & lz_arr = edge_lengths[2]->array(mfi);
#if defined(WARPX_DIM_3D)
amrex::Array4<amrex::Real> const & ly_arr = edge_lengths[1]->array(mfi);
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
amrex::Dim3 const lx_lo = amrex::lbound(lx_arr);
amrex::Dim3 const lx_hi = amrex::ubound(lx_arr);
amrex::Dim3 const lz_lo = amrex::lbound(lz_arr);
amrex::Dim3 const lz_hi = amrex::ubound(lz_arr);
#endif

amrex::ParallelFor (tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// Do not update Ex if the edge on which it lives is fully covered
eb_update_Ex_arr(i, j, k) = (lx_arr(i, j, k) == 0)? 0 : 1;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef WARPX_DIM_3D
// In 3D: Do not update Ey if the edge on which it lives is fully covered
eb_update_Ey_arr(i, j, k) = (ly_arr(i, j, k) == 0)? 0 : 1;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
// In XZ and RZ: Ey is associated with a mesh node,
// so we need to check if the mesh node is covered
if((lx_arr(std::min(i , lx_hi.x), std::min(j , lx_hi.y), k)==0)
||(lx_arr(std::max(i-1, lx_lo.x), std::min(j , lx_hi.y), k)==0)
||(lz_arr(std::min(i , lz_hi.x), std::min(j , lz_hi.y), k)==0)
||(lz_arr(std::min(i , lz_hi.x), std::max(j-1, lz_lo.y), k)==0)) {
eb_update_Ey_arr(i, j, k) = 0;
} else {
eb_update_Ey_arr(i, j, k) = 1;
}
#endif
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// Do not update Ez if the edge on which it lives is fully covered
eb_update_Ez_arr(i, j, k) = (lz_arr(i, j, k) == 0)? 0 : 1;
}
);

}
}

void
WarpX::MarkUpdateBCellsECT (
std::array< std::unique_ptr<amrex::iMultiFab>,3> & eb_update_B,
ablastr::fields::VectorField const& face_areas,
ablastr::fields::VectorField const& edge_lengths )
{

#ifdef AMREX_USE_OMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( amrex::MFIter mfi(*eb_update_B[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {

const amrex::Box& tbx = mfi.tilebox( eb_update_B[0]->ixType().toIntVect(), eb_update_B[0]->nGrowVect() );
const amrex::Box& tby = mfi.tilebox( eb_update_B[1]->ixType().toIntVect(), eb_update_B[1]->nGrowVect() );
const amrex::Box& tbz = mfi.tilebox( eb_update_B[2]->ixType().toIntVect(), eb_update_B[2]->nGrowVect() );

amrex::Array4<int> const & eb_update_Bx_arr = eb_update_B[0]->array(mfi);
amrex::Array4<int> const & eb_update_By_arr = eb_update_B[1]->array(mfi);
amrex::Array4<int> const & eb_update_Bz_arr = eb_update_B[2]->array(mfi);

#ifdef WARPX_DIM_3D
amrex::Array4<amrex::Real> const & Sx_arr = face_areas[0]->array(mfi);
amrex::Array4<amrex::Real> const & Sy_arr = face_areas[1]->array(mfi);
amrex::Array4<amrex::Real> const & Sz_arr = face_areas[2]->array(mfi);
amrex::ignore_unused(edge_lengths);
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
amrex::Array4<amrex::Real> const & Sy_arr = face_areas[1]->array(mfi);
amrex::Array4<amrex::Real> const & lx_arr = edge_lengths[0]->array(mfi);
amrex::Array4<amrex::Real> const & lz_arr = edge_lengths[2]->array(mfi);
#endif
amrex::ParallelFor (tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef WARPX_DIM_3D
// In 3D: do not update Bx if the face on which it lives is fully covered
eb_update_Bx_arr(i, j, k) = (Sx_arr(i, j, k) == 0)? 0 : 1;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
//In XZ and RZ, Bx lives on a z-edge ; do not update if fully covered
eb_update_Bx_arr(i, j, k) = (lz_arr(i, j, k) == 0)? 0 : 1;
#endif
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
// Do not update By if the face on which it lives is fully covered
eb_update_By_arr(i, j, k) = (Sy_arr(i, j, k) == 0)? 0 : 1;
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
#ifdef WARPX_DIM_3D
// In 3D: do not update Bz if the face on which it lives is fully covered
eb_update_Bz_arr(i, j, k) = (Sz_arr(i, j, k) == 0)? 0 : 1;
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
//In XZ and RZ, Bz lives on a x-edge ; do not update if fully covered
eb_update_Bz_arr(i, j, k) = (lx_arr(i, j, k) == 0)? 0 : 1;
#endif
}
);

}
}

void
WarpX::MarkCells ()
WarpX::MarkExtensionCells ()
{
using ablastr::fields::Direction;
using warpx::fields::FieldType;
Expand All @@ -302,7 +517,7 @@ WarpX::MarkCells ()
auto const &cell_size = CellSize(maxLevel());

#if !defined(WARPX_DIM_3D) && !defined(WARPX_DIM_XZ)
WARPX_ABORT_WITH_MESSAGE("MarkCells only implemented in 2D and 3D");
WARPX_ABORT_WITH_MESSAGE("MarkExtensionCells only implemented in 2D and 3D");
#endif

for (int idim = 0; idim < 3; ++idim) {
Expand Down
Loading

0 comments on commit e3378b0

Please sign in to comment.