Skip to content

Commit

Permalink
Updated init functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@RemiLehe
RemiLehe committed Jan 10, 2025
1 parent 8516d1a commit 3c0867e
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 26 deletions.
39 changes: 20 additions & 19 deletions Source/EmbeddedBoundary/WarpXInitEB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,9 +292,10 @@ WarpX::ScaleAreas (ablastr::fields::VectorField& face_areas,
}

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

using ablastr::fields::Direction;
Expand All @@ -305,63 +306,63 @@ WarpX::MarkUpdateECells (

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

for (amrex::MFIter mfi(*m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)); mfi.isValid(); ++mfi) {
for (amrex::MFIter mfi(*field[idim]); mfi.isValid(); ++mfi) {

const amrex::Box& box = mfi.tilebox();
auto const & update_E_arr = m_eb_update_E[lev][idim]->array(mfi);
auto 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 E in every cell
// Every cell in box is all regular: update field in every cell
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
update_E_arr(i, j, k) = 1;
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 E
// Every cell in box is all covered: do not update field
amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
update_E_arr(i, j, k) = 0;
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 = m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)->ixType();
auto index_type = field[idim]->ixType();

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

// Stair-case approximation: If neighboring cells are either partially
// or fully covered (i.e. if they are not regular cells): do not update E
// or fully covered (i.e. if they are not regular cells): do not update field

int update_E = 1;
int eb_update = 1;
// Check neighbors in the first spatial direction
if ( index_type.nodeCentered(0) ) {
if ( !flag(i, j, k).isRegular() || !flag(i-1, j, k).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() || !flag(i-1, j, k).isRegular() ) { eb_update = 0; }
} else {
if ( !flag(i, j, k).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() ) { eb_update = 0; }
}
#if AMREX_SPACEDIM > 1
// Check neighbors in the second spatial direction
if ( index_type.nodeCentered(1) ) {
if ( !flag(i, j, k).isRegular() || !flag(i, j-1, k).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() || !flag(i, j-1, k).isRegular() ) { eb_update = 0; }
} else {
if ( !flag(i, j, k).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() ) { eb_update = 0; }
}
#endif
#if AMREX_SPACEDIM > 2
// Check neighbors in the third spatial direction
if ( index_type.nodeCentered(2) ) {
if ( !flag(i, j, k).isRegular() || !flag(i, j, k-1).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() || !flag(i, j, k-1).isRegular() ) { eb_update = 0; }
} else {
if ( !flag(i, j, k).isRegular() ) { update_E = 0; }
if ( !flag(i, j, k).isRegular() ) { eb_update = 0; }
}
#endif
update_E_arr(i, j, k) = update_E;
eb_update_arr(i, j, k) = eb_update;
});

// TODO: Handle the case of the ECT solver
Expand Down
6 changes: 5 additions & 1 deletion Source/Initialization/WarpXInitData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1300,7 +1300,11 @@ void WarpX::InitializeEBGridData (int lev)

auto const eb_fact = fieldEBFactory(lev);

MarkUpdateECells( eb_fact, lev );
// Mark on which grid points E should be updated
MarkUpdateCells(
m_eb_update_E[lev],
m_fields.get_alldirs(FieldType::Efield_fp, lev),
eb_fact );

// TODO: move inside if condition for ECT
auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev);
Expand Down
24 changes: 18 additions & 6 deletions Source/WarpX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1013,13 +1013,25 @@ public:
void InitEB ();

#ifdef AMREX_USE_EB
/** Set a flag to indicate which cells are in the valid domain and should be updated
* by the field solve, and which cells are in the EB and should not be updated.
* TODO: Add detail about criterion for Yee and ECT solver
/** Set a flag to indicate on which grid points the field `field`
* should be updated.
*
* More specifically, this function fill the iMultiFabs in `eb_update`
* (which have the same indexType as the MultiFabs in `field`) with 1
* or 0, depending on whether the grid point is in the valid region
* or in the embedded boundary.
*
* This uses a stair-case approximation of the embedded boundary
* (unless the ECT solver is used): If a grid point touches cells
* that are either partially or fully covered by the embedded
* boundary: do not update field.
*
* TODO Discuss ECT solver
*/
void MarkUpdateECells (
amrex::EBFArrayBoxFactory const & eb_fact,
int const lev );
void MarkUpdateCells (
std::array< std::unique_ptr<amrex::iMultiFab>,3> & eb_update,
ablastr::fields::VectorField const & field,
amrex::EBFArrayBoxFactory const & eb_fact );

/**
* \brief Compute the length of the mesh edges. Here the length is a value in [0, 1].
Expand Down

0 comments on commit 3c0867e

Please sign in to comment.