From 3c0867eff6e622dc5aa0a96e60f28abafb264f2d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2025 06:12:49 -0800 Subject: [PATCH] Updated init functions --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 39 +++++++++++++------------ Source/Initialization/WarpXInitData.cpp | 6 +++- Source/WarpX.H | 24 +++++++++++---- 3 files changed, 43 insertions(+), 26 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 29dd1da9279..059667d605c 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -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,3> & eb_update, + ablastr::fields::VectorField const& field, + amrex::EBFArrayBoxFactory const & eb_fact ) { using ablastr::fields::Direction; @@ -305,10 +306,10 @@ 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); @@ -316,52 +317,52 @@ WarpX::MarkUpdateECells ( 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 diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index b0f429b3e74..945376a8eae 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -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); diff --git a/Source/WarpX.H b/Source/WarpX.H index bb771c18ece..063555ea7a7 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -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,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].