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.

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

[Do not merge] Test checksums with new EB structures #5574

Open
wants to merge 33 commits into
base: development
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
3fa996f
Modify stair-case approximation
RemiLehe Jan 3, 2025
112cce1
Fix compilation without EB
RemiLehe Jan 5, 2025
b0ccce4
Slightly modify boundaries
RemiLehe Jan 6, 2025
76fa82e
Added new arrays with flags
RemiLehe Jan 8, 2025
3624a72
Finalize initialization of the array
RemiLehe Jan 8, 2025
3181d45
Use guard cell information
RemiLehe Jan 8, 2025
38d967d
Remove edge_length in field update
RemiLehe Jan 9, 2025
79ed114
Fix const-ness
RemiLehe Jan 9, 2025
8516d1a
Fix compilation warning
RemiLehe Jan 9, 2025
3c0867e
Updated init functions
RemiLehe Jan 10, 2025
28d43a1
Fix initialization criterion
RemiLehe Jan 10, 2025
0214769
Add flag for B
RemiLehe Jan 10, 2025
ab87809
Fix RZ compilation error
RemiLehe Jan 10, 2025
78bb8c1
Update initialization of external fields
RemiLehe Jan 11, 2025
df9ae82
Fix setting of eb_update
RemiLehe Jan 11, 2025
59579b4
Implement ECT condition for B
RemiLehe Jan 11, 2025
4fd4dc7
Do not use any guard cells
RemiLehe Jan 11, 2025
4232db7
Fixed ECT solver
RemiLehe Jan 11, 2025
4b8079f
Merge branch 'move_stair_case_approx' of github.com:RemiLehe/WarpX in…
RemiLehe Jan 11, 2025
f45b257
Revert "Do not use any guard cells"
RemiLehe Jan 12, 2025
2003d85
Fix compilation errors
RemiLehe Jan 12, 2025
d6e64af
Fix compilation error
RemiLehe Jan 12, 2025
e6d93f9
Fix const-ness
RemiLehe Jan 12, 2025
eb43087
Activate load-balancing
RemiLehe Jan 12, 2025
2eb3383
Update checksums
RemiLehe Jan 12, 2025
8981179
Add documentation
RemiLehe Jan 12, 2025
8182df7
Apply suggestions from code review
RemiLehe Jan 14, 2025
0980516
Minor cleanup
RemiLehe Jan 17, 2025
ace399b
Modify initialization in ECT solver
RemiLehe Jan 18, 2025
16ec9b2
Update comments
RemiLehe Jan 18, 2025
f2b0580
Revert "Update checksums"
RemiLehe Jan 17, 2025
3a3f36e
Use previous criterion to update EB
RemiLehe Jan 17, 2025
4b4e761
Use previous input scripts
RemiLehe Jan 17, 2025
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
4 changes: 2 additions & 2 deletions Examples/Tests/embedded_boundary_cube/inputs_base_3d
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
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
// 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
Loading