From 3fa996fc8835fcadcd6c523ca56b0ba5f4f8cd03 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 2 Jan 2025 16:54:39 -0800 Subject: [PATCH 01/32] Modify stair-case approximation --- .../FiniteDifferenceSolver/EvolveE.cpp | 77 +++++++++++++++---- 1 file changed, 63 insertions(+), 14 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 03a9866fb98..6d30ad8bdca 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -155,11 +155,24 @@ void FiniteDifferenceSolver::EvolveECartesian ( Array4 const& jy = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); - amrex::Array4 lx, ly, lz; + // Extract structures for embedded boundaries + bool eb_fully_covered_box = false; + amrex::Array4::value_type> eb_flag_arr; if (EB::enabled()) { - lx = edge_lengths[0]->array(mfi); - ly = edge_lengths[1]->array(mfi); - lz = edge_lengths[2]->array(mfi); +#ifdef AMREX_USE_EB + auto & warpx = WarpX::GetInstance(); + amrex::EBFArrayBoxFactory const& eb_box_factory = warpx.fieldEBFactory(lev); + amrex::FabArray const& eb_flag = eb_box_factory.getMultiEBCellFlagFab(); + amrex::Box const& box = mfi.tilebox( amrex::IntVect::TheCellVector() ); + amrex::FabType const fab_type = eb_flag[mfi].getType(box); + if (fab_type == amrex::FabType::covered) { + eb_fully_covered_box = true; + } else if (!(fab_type == amrex::FabType::regular)) { + // For cells that are not fully covered or regular, + // we need to check each cell is covered or regular + eb_flag_arr = eb_flag.array(mfi); + } +#endif } // Extract stencil coefficients @@ -179,8 +192,21 @@ void FiniteDifferenceSolver::EvolveECartesian ( amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip field push if this cell is fully covered by embedded boundaries - if (lx && lx(i, j, k) <= 0) { return; } + + // Stair-case approximation to the embedded boundary: + // Skip field push if this edge touches a partially or fully covered cell + if (eb_fully_covered_box) { return; } + else if (eb_flag_arr) { +#ifdef WARPX_DIM_3D + if ( !eb_flag_arr(i, j-1, k-1).isRegular() || + !eb_flag_arr(i, j-1, k ).isRegular() || + !eb_flag_arr(i, j , k-1).isRegular() || + !eb_flag_arr(i, j , k ).isRegular() ) { return; } +#elif (defined WARPX_DIM_XZ) + if ( !eb_flag_arr(i, j-1, k).isRegular() || + !eb_flag_arr(i, j , k).isRegular() ) { return; } +#endif + } Ex(i, j, k) += c2 * dt * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k) @@ -189,14 +215,23 @@ void FiniteDifferenceSolver::EvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip field push if this cell is fully covered by embedded boundaries + + // Stair-case approximation to the embedded boundary: + // Skip field push if this edge touches a partially or fully covered cell + if (eb_fully_covered_box) { return; } + else if (eb_flag_arr) { #ifdef WARPX_DIM_3D - if (ly && ly(i,j,k) <= 0) { return; } -#elif defined(WARPX_DIM_XZ) - //In XZ Ey is associated with a mesh node, so we need to check if the mesh node is covered - amrex::ignore_unused(ly); - if (lx && (lx(i, j, k)<=0 || lx(i-1, j, k)<=0 || lz(i, j-1, k)<=0 || lz(i, j, k)<=0)) { return; } + if ( !eb_flag_arr(i-1, j, k-1).isRegular() || + !eb_flag_arr(i-1, j, k ).isRegular() || + !eb_flag_arr(i , j, k-1).isRegular() || + !eb_flag_arr(i , j, k ).isRegular() ) { return; } +#elif (defined WARPX_DIM_XZ) + if ( !eb_flag_arr(i-1, j-1, k).isRegular() || + !eb_flag_arr(i-1, j , k).isRegular() || + !eb_flag_arr(i , j-1, k).isRegular() || + !eb_flag_arr(i , j , k).isRegular() ) { return; } #endif + } Ey(i, j, k) += c2 * dt * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) @@ -205,8 +240,22 @@ void FiniteDifferenceSolver::EvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip field push if this cell is fully covered by embedded boundaries - if (lz && lz(i,j,k) <= 0) { return; } + + // Stair-case approximation to the embedded boundary: + // Skip field push if this edge touches a partially or fully covered cell + if (eb_fully_covered_box) { return; } + else if (eb_flag_arr) { +#ifdef WARPX_DIM_3D + if ( !eb_flag_arr(i-1, j-1, k).isRegular() || + !eb_flag_arr(i-1, j , k).isRegular() || + !eb_flag_arr(i , j-1, k).isRegular() || + !eb_flag_arr(i , j , k).isRegular() ) { return; } +#elif (defined WARPX_DIM_XZ) + if ( !eb_flag_arr(i-1, j, k).isRegular() || + !eb_flag_arr(i , j, k).isRegular() ) { return; } +#endif + } + Ez(i, j, k) += c2 * dt * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k) + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k) From 112cce1952f1b1f5753dfcf876351c6a33cfb0af Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 5 Jan 2025 08:10:14 -0800 Subject: [PATCH 02/32] Fix compilation without EB --- Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 6d30ad8bdca..bddd69b2393 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -155,11 +155,11 @@ void FiniteDifferenceSolver::EvolveECartesian ( Array4 const& jy = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); +#ifdef AMREX_USE_EB // Extract structures for embedded boundaries bool eb_fully_covered_box = false; amrex::Array4::value_type> eb_flag_arr; if (EB::enabled()) { -#ifdef AMREX_USE_EB auto & warpx = WarpX::GetInstance(); amrex::EBFArrayBoxFactory const& eb_box_factory = warpx.fieldEBFactory(lev); amrex::FabArray const& eb_flag = eb_box_factory.getMultiEBCellFlagFab(); @@ -172,8 +172,8 @@ void FiniteDifferenceSolver::EvolveECartesian ( // we need to check each cell is covered or regular eb_flag_arr = eb_flag.array(mfi); } -#endif } +#endif // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); @@ -193,6 +193,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB // Stair-case approximation to the embedded boundary: // Skip field push if this edge touches a partially or fully covered cell if (eb_fully_covered_box) { return; } @@ -207,6 +208,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( !eb_flag_arr(i, j , k).isRegular() ) { return; } #endif } +#endif Ex(i, j, k) += c2 * dt * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k) @@ -216,6 +218,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB // Stair-case approximation to the embedded boundary: // Skip field push if this edge touches a partially or fully covered cell if (eb_fully_covered_box) { return; } @@ -232,6 +235,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( !eb_flag_arr(i , j , k).isRegular() ) { return; } #endif } +#endif Ey(i, j, k) += c2 * dt * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) @@ -241,6 +245,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB // Stair-case approximation to the embedded boundary: // Skip field push if this edge touches a partially or fully covered cell if (eb_fully_covered_box) { return; } @@ -255,6 +260,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( !eb_flag_arr(i , j, k).isRegular() ) { return; } #endif } +#endif Ez(i, j, k) += c2 * dt * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k) From b0ccce457a74054099c93454f71f7de04b2b5dd9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 6 Jan 2025 11:28:46 -0800 Subject: [PATCH 03/32] Slightly modify boundaries --- .../inputs_test_2d_embedded_boundary_cube | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube index 684325dc030..4fbf13e7f5c 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube +++ b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube @@ -12,10 +12,10 @@ warpx.abort_on_warning_threshold = medium boundary.field_lo = pec pec boundary.field_hi = pec pec -my_constants.xmin = -0.5 -my_constants.zmin = -0.5 -my_constants.xmax = 0.5 -my_constants.zmax = 0.5 +my_constants.xmin = -0.501 +my_constants.zmin = -0.501 +my_constants.xmax = 0.501 +my_constants.zmax = 0.501 # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular. warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))" From 76fa82e783c2420ca9b589db2b48da838532d1ea Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 06:46:16 -0800 Subject: [PATCH 04/32] Added new arrays with flags --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 55 ++++++++++++++++++++++++- Source/Initialization/WarpXInitData.cpp | 5 ++- Source/WarpX.H | 19 +++++++-- Source/WarpX.cpp | 10 +++++ 4 files changed, 83 insertions(+), 6 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index edbc97a8efe..d87b891dd99 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -291,9 +291,60 @@ WarpX::ScaleAreas (ablastr::fields::VectorField& face_areas, } } +void +WarpX::MarkUpdateECells ( + amrex::EBFArrayBoxFactory const & eb_fact, + int const lev ) +{ + + using ablastr::fields::Direction; + using warpx::fields::FieldType; + + // Extract structures for embedded boundaries + amrex::FabArray const& eb_flag = eb_fact.getMultiEBCellFlagFab(); + + for (int idim = 0; idim < 3; ++idim) { + + // TODO: Handle guard cells + + for (amrex::MFIter mfi(*m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)); mfi.isValid(); ++mfi) { + + const amrex::Box& box = mfi.tilebox(); + auto const & update_E_arr = m_eb_update_E[lev][idim]->array(mfi); + + amrex::FabType const fab_type = eb_flag[mfi].getType(mfi.tilebox(amrex::IntVect::TheCellVector())); + + if (fab_type == amrex::FabType::regular) { + + // every cell in box is all regular + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + update_E_arr(i, j, k) = 1; + }); + + } else if (fab_type == amrex::FabType::covered) { + + // every cell in box is all covered + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + update_E_arr(i, j, k) = 0; + }); + + } else { + + auto const & flag = eb_flag[mfi].array(); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + update_E_arr(i, j, k) = flag(i, j, k).isRegular(); + }); + + } + + } + + } + +} void -WarpX::MarkCells () +WarpX::MarkExtensionCells () { using ablastr::fields::Direction; using warpx::fields::FieldType; @@ -302,7 +353,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) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index daecfac8bed..b0f429b3e74 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1300,6 +1300,9 @@ void WarpX::InitializeEBGridData (int lev) auto const eb_fact = fieldEBFactory(lev); + MarkUpdateECells( eb_fact, lev ); + + // TODO: move inside if condition for ECT auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev); ComputeEdgeLengths(edge_lengths_lev, eb_fact); ScaleEdges(edge_lengths_lev, CellSize(lev)); @@ -1309,7 +1312,7 @@ void WarpX::InitializeEBGridData (int lev) ScaleAreas(face_areas_lev, CellSize(lev)); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { - MarkCells(); + MarkExtensionCells(); ComputeFaceExtensions(); } } diff --git a/Source/WarpX.H b/Source/WarpX.H index 73998d6faf2..bb771c18ece 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1013,6 +1013,14 @@ 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 + */ + void MarkUpdateECells ( + amrex::EBFArrayBoxFactory const & eb_fact, + int const lev ); + /** * \brief Compute the length of the mesh edges. Here the length is a value in [0, 1]. * An edge of length 0 is fully covered. @@ -1044,7 +1052,7 @@ public: * - 2 for stable cells which have been intruded * Here we cannot know if a cell is intruded or not so we initialize all stable cells with 1 */ - void MarkCells(); + void MarkExtensionCells(); #endif /** @@ -1393,17 +1401,22 @@ private: mutable amrex::Vector,3 > > Afield_dotMask; mutable amrex::Vector< std::unique_ptr > phi_dotMask; + /** EB: Flag to indicate whether the E location is inside the EB and therefore E should not be updated. + * (One array per level and per direction, due to staggering) + */ + amrex::Vector,3 > > m_eb_update_E; + /** EB: for every mesh face flag_info_face contains a: * * 0 if the face needs to be extended * * 1 if the face is large enough to lend area to other faces * * 2 if the face is actually intruded by other face - * It is initialized in WarpX::MarkCells + * It is initialized in WarpX::MarkExtensionCells * This is only used for the ECT solver.*/ amrex::Vector, 3 > > m_flag_info_face; /** EB: for every mesh face face flag_ext_face contains a: * * 1 if the face needs to be extended * * 0 otherwise - * It is initialized in WarpX::MarkCells and then modified in WarpX::ComputeOneWayExtensions + * It is initialized in WarpX::MarkExtensionCells and then modified in WarpX::ComputeOneWayExtensions * and in WarpX::ComputeEightWaysExtensions * This is only used for the ECT solver.*/ amrex::Vector, 3 > > m_flag_ext_face; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 965235e1078..6b2c025c89f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -323,6 +323,7 @@ WarpX::WarpX () Afield_dotMask.resize(nlevs_max); phi_dotMask.resize(nlevs_max); + m_eb_update_E.resize(nlevs_max); m_flag_info_face.resize(nlevs_max); m_flag_ext_face.resize(nlevs_max); m_borrowing.resize(nlevs_max); @@ -2301,6 +2302,15 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // EB info are needed only at the finest level if (lev == maxLevel()) { if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { + + AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]"); + AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]"); + AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); + + // TODO: do not allocate these arrays anymore with the Yee solver //! EB: Lengths of the mesh edges m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); From 3624a7232e5cb9551f0d9b55dc65888f1ac6e4ca Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 07:01:48 -0800 Subject: [PATCH 05/32] Finalize initialization of the array --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 39 +++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index d87b891dd99..676c28a1e9a 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -316,14 +316,15 @@ WarpX::MarkUpdateECells ( if (fab_type == amrex::FabType::regular) { - // every cell in box is all regular + // Every cell in box is all regular: update E in every cell + // TODO: We actually need to check guard cells amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { update_E_arr(i, j, k) = 1; }); } else if (fab_type == amrex::FabType::covered) { - // every cell in box is all covered + // Every cell in box is all covered: do not update E amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { update_E_arr(i, j, k) = 0; }); @@ -331,10 +332,42 @@ WarpX::MarkUpdateECells ( } else { auto const & flag = eb_flag[mfi].array(); + auto index_type = m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)->ixType(); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - update_E_arr(i, j, k) = flag(i, j, k).isRegular(); + + // If neighboring cells are partially covered: do not update E + int update_E = 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; } + } else { + if ( !flag(i, j, k).isRegular() ) { update_E = 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; } + } else { + if ( !flag(i, j, k).isRegular() ) { update_E = 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; } + } else { + if ( !flag(i, j, k).isRegular() ) { update_E = 0; } + } +#endif + update_E_arr(i, j, k) = update_E; }); + // TODO: Handle the case of the ECT solver + } } From 3181d45e840eccc8222ccd90e5de5b3d4cabfbe3 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 11:38:32 -0800 Subject: [PATCH 06/32] Use guard cell information --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 676c28a1e9a..1259c9f8835 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -305,38 +305,38 @@ WarpX::MarkUpdateECells ( for (int idim = 0; idim < 3; ++idim) { - // TODO: Handle guard cells - for (amrex::MFIter mfi(*m_fields.get(FieldType::Efield_fp, Direction{idim}, lev)); mfi.isValid(); ++mfi) { const amrex::Box& box = mfi.tilebox(); auto const & update_E_arr = m_eb_update_E[lev][idim]->array(mfi); - amrex::FabType const fab_type = eb_flag[mfi].getType(mfi.tilebox(amrex::IntVect::TheCellVector())); + // 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) { + 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 - // TODO: We actually need to check guard cells amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { update_E_arr(i, j, k) = 1; }); - } else if (fab_type == amrex::FabType::covered) { + } 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 amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { update_E_arr(i, j, k) = 0; }); - } else { + } 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(); amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - // If neighboring cells are partially covered: do not update E + // Stair-case approximation: If neighboring cells are either partially + // or fully covered (i.e. if they are not regular cells): do not update E int update_E = 1; // Check neighbors in the first spatial direction From 38d967d83b3f1dbc99114d1a329f6a8658334078 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 16:28:05 -0800 Subject: [PATCH 07/32] Remove edge_length in field update --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 4 +- .../FiniteDifferenceSolver/EvolveE.cpp | 136 +++++------------- .../FiniteDifferenceSolver.H | 5 +- .../ImplicitSolvers/WarpXImplicitOps.cpp | 2 + Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 + 5 files changed, 40 insertions(+), 109 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 1259c9f8835..29dd1da9279 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -337,15 +337,14 @@ WarpX::MarkUpdateECells ( // Stair-case approximation: If neighboring cells are either partially // or fully covered (i.e. if they are not regular cells): do not update E - int update_E = 1; + int update_E = 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; } } else { if ( !flag(i, j, k).isRegular() ) { update_E = 0; } } - #if AMREX_SPACEDIM > 1 // Check neighbors in the second spatial direction if ( index_type.nodeCentered(1) ) { @@ -354,7 +353,6 @@ WarpX::MarkUpdateECells ( if ( !flag(i, j, k).isRegular() ) { update_E = 0; } } #endif - #if AMREX_SPACEDIM > 2 // Check neighbors in the third spatial direction if ( index_type.nodeCentered(2) ) { diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index bddd69b2393..8f9d13dcb45 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -55,6 +55,7 @@ void FiniteDifferenceSolver::EvolveE ( int lev, PatchType patch_type, ablastr::fields::VectorField const& Efield, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::Real const dt ) { @@ -72,40 +73,23 @@ void FiniteDifferenceSolver::EvolveE ( fields.get(FieldType::F_fp, lev) : fields.get(FieldType::F_cp, lev); } - ablastr::fields::VectorField edge_lengths; - if (fields.has_vector(FieldType::edge_lengths, lev)) { - edge_lengths = fields.get_alldirs(FieldType::edge_lengths, lev); - } - ablastr::fields::VectorField face_areas; - if (fields.has_vector(FieldType::face_areas, lev)) { - face_areas = fields.get_alldirs(FieldType::face_areas, lev); - } - ablastr::fields::VectorField area_mod; - if (fields.has_vector(FieldType::area_mod, lev)) { - area_mod = fields.get_alldirs(FieldType::area_mod, lev); - } - ablastr::fields::VectorField ECTRhofield; - if (fields.has_vector(FieldType::ECTRhofield, lev)) { - ECTRhofield = fields.get_alldirs(FieldType::ECTRhofield, lev); - } - // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ if (m_fdtd_algo == ElectromagneticSolverAlgo::Yee){ - EvolveECylindrical ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); + EvolveECylindrical ( Efield, Bfield, Jfield, eb_update_E, Ffield, lev, dt ); #else if (m_grid_type == GridType::Collocated) { - EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, eb_update_E, Ffield, lev, dt ); } else if (m_fdtd_algo == ElectromagneticSolverAlgo::Yee || m_fdtd_algo == ElectromagneticSolverAlgo::ECT) { - EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, eb_update_E, Ffield, lev, dt ); } else if (m_fdtd_algo == ElectromagneticSolverAlgo::CKC) { - EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, eb_update_E, Ffield, lev, dt ); #endif } else { @@ -122,7 +106,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - VectorField const& edge_lengths, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real const dt ) { @@ -155,25 +139,13 @@ void FiniteDifferenceSolver::EvolveECartesian ( Array4 const& jy = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); -#ifdef AMREX_USE_EB - // Extract structures for embedded boundaries - bool eb_fully_covered_box = false; - amrex::Array4::value_type> eb_flag_arr; + // Extract structures indicating whether the E field should be updated + amrex::Array4 update_Ex_arr, update_Ey_arr, update_Ez_arr; if (EB::enabled()) { - auto & warpx = WarpX::GetInstance(); - amrex::EBFArrayBoxFactory const& eb_box_factory = warpx.fieldEBFactory(lev); - amrex::FabArray const& eb_flag = eb_box_factory.getMultiEBCellFlagFab(); - amrex::Box const& box = mfi.tilebox( amrex::IntVect::TheCellVector() ); - amrex::FabType const fab_type = eb_flag[mfi].getType(box); - if (fab_type == amrex::FabType::covered) { - eb_fully_covered_box = true; - } else if (!(fab_type == amrex::FabType::regular)) { - // For cells that are not fully covered or regular, - // we need to check each cell is covered or regular - eb_flag_arr = eb_flag.array(mfi); - } + update_Ex_arr = eb_update_E[0]->array(mfi); + update_Ey_arr = eb_update_E[1]->array(mfi); + update_Ez_arr = eb_update_E[2]->array(mfi); } -#endif // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); @@ -193,22 +165,8 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ -#ifdef AMREX_USE_EB - // Stair-case approximation to the embedded boundary: - // Skip field push if this edge touches a partially or fully covered cell - if (eb_fully_covered_box) { return; } - else if (eb_flag_arr) { -#ifdef WARPX_DIM_3D - if ( !eb_flag_arr(i, j-1, k-1).isRegular() || - !eb_flag_arr(i, j-1, k ).isRegular() || - !eb_flag_arr(i, j , k-1).isRegular() || - !eb_flag_arr(i, j , k ).isRegular() ) { return; } -#elif (defined WARPX_DIM_XZ) - if ( !eb_flag_arr(i, j-1, k).isRegular() || - !eb_flag_arr(i, j , k).isRegular() ) { return; } -#endif - } -#endif + // Skip field push in the embedded boundaries + if (update_Ex_arr && update_Ex_arr(i, j, k) == 0) { return; } Ex(i, j, k) += c2 * dt * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k) @@ -218,24 +176,8 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ -#ifdef AMREX_USE_EB - // Stair-case approximation to the embedded boundary: - // Skip field push if this edge touches a partially or fully covered cell - if (eb_fully_covered_box) { return; } - else if (eb_flag_arr) { -#ifdef WARPX_DIM_3D - if ( !eb_flag_arr(i-1, j, k-1).isRegular() || - !eb_flag_arr(i-1, j, k ).isRegular() || - !eb_flag_arr(i , j, k-1).isRegular() || - !eb_flag_arr(i , j, k ).isRegular() ) { return; } -#elif (defined WARPX_DIM_XZ) - if ( !eb_flag_arr(i-1, j-1, k).isRegular() || - !eb_flag_arr(i-1, j , k).isRegular() || - !eb_flag_arr(i , j-1, k).isRegular() || - !eb_flag_arr(i , j , k).isRegular() ) { return; } -#endif - } -#endif + // Skip field push in the embedded boundaries + if (update_Ey_arr && update_Ey_arr(i, j, k) == 0) { return; } Ey(i, j, k) += c2 * dt * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) @@ -245,22 +187,8 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] AMREX_GPU_DEVICE (int i, int j, int k){ -#ifdef AMREX_USE_EB - // Stair-case approximation to the embedded boundary: - // Skip field push if this edge touches a partially or fully covered cell - if (eb_fully_covered_box) { return; } - else if (eb_flag_arr) { -#ifdef WARPX_DIM_3D - if ( !eb_flag_arr(i-1, j-1, k).isRegular() || - !eb_flag_arr(i-1, j , k).isRegular() || - !eb_flag_arr(i , j-1, k).isRegular() || - !eb_flag_arr(i , j , k).isRegular() ) { return; } -#elif (defined WARPX_DIM_XZ) - if ( !eb_flag_arr(i-1, j, k).isRegular() || - !eb_flag_arr(i , j, k).isRegular() ) { return; } -#endif - } -#endif + // Skip field push in the embedded boundaries + if (update_Ez_arr && update_Ez_arr(i, j, k) == 0) { return; } Ez(i, j, k) += c2 * dt * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k) @@ -311,14 +239,10 @@ void FiniteDifferenceSolver::EvolveECylindrical ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real const dt ) { -#ifndef AMREX_USE_EB - amrex::ignore_unused(edge_lengths); -#endif - amrex::LayoutData* cost = WarpX::getCosts(lev); // Loop through the grids, and over the tiles within each grid @@ -343,10 +267,12 @@ void FiniteDifferenceSolver::EvolveECylindrical ( Array4 const& jt = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); - amrex::Array4 lr, lz; + // Extract structures indicating whether the E field should be updated + amrex::Array4 update_Er_arr, update_Et_arr, update_Ez_arr; if (EB::enabled()) { - lr = edge_lengths[0]->array(mfi); - lz = edge_lengths[2]->array(mfi); + update_Er_arr = eb_update_E[0]->array(mfi); + update_Et_arr = eb_update_E[1]->array(mfi); + update_Ez_arr = eb_update_E[2]->array(mfi); } // Extract stencil coefficients @@ -371,8 +297,9 @@ void FiniteDifferenceSolver::EvolveECylindrical ( amrex::ParallelFor(ter, tet, tez, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip field push if this cell is fully covered by embedded boundaries - if (lr && lr(i, j, 0) <= 0) { return; } + + // Skip field push in the embedded boundaries + if (update_Er_arr && update_Er_arr(i, j, k) == 0) { return; } Real const r = rmin + (i + 0.5_rt)*dr; // r on cell-centered point (Er is cell-centered in r) Er(i, j, 0, 0) += c2 * dt*( @@ -391,9 +318,9 @@ void FiniteDifferenceSolver::EvolveECylindrical ( }, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip field push if this cell is fully covered by embedded boundaries - // The Et field is at a node, so we need to check if the node is covered - if (lr && (lr(i, j, 0)<=0 || lr(i-1, j, 0)<=0 || lz(i, j-1, 0)<=0 || lz(i, j, 0)<=0)) { return; } + + // Skip field push in the embedded boundaries + if (update_Et_arr && update_Et_arr(i, j, k) == 0) { return; } Real const r = rmin + i*dr; // r on a nodal grid (Et is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations @@ -436,8 +363,9 @@ void FiniteDifferenceSolver::EvolveECylindrical ( }, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip field push if this cell is fully covered by embedded boundaries - if (lz && lz(i, j, 0) <= 0) { return; } + + // Skip field push in the embedded boundaries + if (update_Ez_arr && update_Ez_arr(i, j, k) == 0) { return; } Real const r = rmin + i*dr; // r on a nodal grid (Ez is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 45c06584fda..64d93f33935 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -63,6 +63,7 @@ class FiniteDifferenceSolver int lev, PatchType patch_type, ablastr::fields::VectorField const& Efield, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::Real dt ); void EvolveF ( amrex::MultiFab* Ffield, @@ -215,7 +216,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real dt ); @@ -267,7 +268,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 >& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real dt ); diff --git a/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp index eaf96cf77ec..9b62bd91b0c 100644 --- a/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp +++ b/Source/FieldSolver/ImplicitSolvers/WarpXImplicitOps.cpp @@ -385,12 +385,14 @@ WarpX::ImplicitComputeRHSE (int lev, PatchType patch_type, amrex::Real a_dt, War lev, patch_type, a_Erhs_vec.getArrayVec()[lev], + m_eb_update_E[lev], a_dt ); } else { m_fdtd_solver_cp[lev]->EvolveE( m_fields, lev, patch_type, a_Erhs_vec.getArrayVec()[lev], + m_eb_update_E[lev], a_dt ); } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 24640fc63c7..a77fb390fff 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -963,12 +963,14 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt, amrex::Real sta lev, patch_type, m_fields.get_alldirs(FieldType::Efield_fp, lev), + m_eb_update_E[lev], a_dt ); } else { m_fdtd_solver_cp[lev]->EvolveE( m_fields, lev, patch_type, m_fields.get_alldirs(FieldType::Efield_cp, lev), + m_eb_update_E[lev], a_dt ); } From 79ed11409f3f1fb521dd4dade7e333ef731c5bff Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 16:38:13 -0800 Subject: [PATCH 08/32] Fix const-ness --- Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp | 6 +++--- .../FiniteDifferenceSolver/FiniteDifferenceSolver.H | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 8f9d13dcb45..f8bc4ade0b3 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -55,7 +55,7 @@ void FiniteDifferenceSolver::EvolveE ( int lev, PatchType patch_type, ablastr::fields::VectorField const& Efield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real const dt ) { @@ -106,7 +106,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3> const& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real const dt ) { @@ -239,7 +239,7 @@ void FiniteDifferenceSolver::EvolveECylindrical ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real const dt ) { diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 64d93f33935..07596dc8562 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -63,7 +63,7 @@ class FiniteDifferenceSolver int lev, PatchType patch_type, ablastr::fields::VectorField const& Efield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real dt ); void EvolveF ( amrex::MultiFab* Ffield, @@ -216,7 +216,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real dt ); @@ -268,7 +268,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - std::array< std::unique_ptr,3 >& eb_update_E, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::MultiFab const* Ffield, int lev, amrex::Real dt ); From 8516d1ad8defd31f2dc958386d2de49ba6f21174 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 16:55:51 -0800 Subject: [PATCH 09/32] Fix compilation warning --- Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index f8bc4ade0b3..509027c5569 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -110,10 +110,6 @@ void FiniteDifferenceSolver::EvolveECartesian ( amrex::MultiFab const* Ffield, int lev, amrex::Real const dt ) { -#ifndef AMREX_USE_EB - amrex::ignore_unused(edge_lengths); -#endif - amrex::LayoutData* cost = WarpX::getCosts(lev); Real constexpr c2 = PhysConst::c * PhysConst::c; From 3c0867eff6e622dc5aa0a96e60f28abafb264f2d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2025 06:12:49 -0800 Subject: [PATCH 10/32] 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]. From 28d43a13825aea0642596b14226339e4fc164213 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2025 06:21:16 -0800 Subject: [PATCH 11/32] Fix initialization criterion --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 44 +++++++++++++------------ 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 059667d605c..d2a2c991165 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -337,31 +337,33 @@ WarpX::MarkUpdateCells ( 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 field - - 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() ) { eb_update = 0; } - } else { - if ( !flag(i, j, k).isRegular() ) { eb_update = 0; } - } + // or fully covered: do not update field + + // Check neighbors in the each spatial direction + int i_start = 0; + int j_start = 0; + int k_start = 0; + // If nodal in a given direction, we need to start from -1 (left-neighboring cell) + if ( index_type.nodeCentered(0) ) { i_start = -1; }; #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() ) { eb_update = 0; } - } else { - if ( !flag(i, j, k).isRegular() ) { eb_update = 0; } - } + if ( index_type.nodeCentered(1) ) { j_start = -1; }; #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() ) { eb_update = 0; } - } else { - if ( !flag(i, j, k).isRegular() ) { eb_update = 0; } - } + if ( index_type.nodeCentered(2) ) { k_start = -1; }; #endif + // Loop over neighboring cells + int eb_update = 1; + for (int ii = i_start; ii <= 1; ++ii) { + for (int jj = j_start; jj <= 1; ++jj) { + for (int kk = k_start; kk <= 1; ++kk) { + // If one of the neighboring is either partially or fully covered + // (i.e. if they are not regular cells), do not update field + if ( !flag(i+ii, j+jj, k+kk).isRegular() ) { + eb_update = 0; + } + } + } + } eb_update_arr(i, j, k) = eb_update; }); From 02147690ebfda487ba0f8709d8b0f5727f1d0061 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2025 06:29:04 -0800 Subject: [PATCH 12/32] Add flag for B --- Source/Initialization/WarpXInitData.cpp | 5 +++++ Source/WarpX.H | 5 +++-- Source/WarpX.cpp | 8 ++++++++ 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 945376a8eae..1b541f85d8e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1305,6 +1305,11 @@ void WarpX::InitializeEBGridData (int lev) m_eb_update_E[lev], m_fields.get_alldirs(FieldType::Efield_fp, lev), eb_fact ); + // Mark on which grid points B should be updated + MarkUpdateCells( + m_eb_update_B[lev], + m_fields.get_alldirs(FieldType::Bfield_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 063555ea7a7..5bd823ad8ba 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1413,10 +1413,11 @@ private: mutable amrex::Vector,3 > > Afield_dotMask; mutable amrex::Vector< std::unique_ptr > phi_dotMask; - /** EB: Flag to indicate whether the E location is inside the EB and therefore E should not be updated. - * (One array per level and per direction, due to staggering) + /** EB: Flag to indicate whether a gridpoint is inside the embedded boundary and therefore + * whether the E or B should not be updated. (One array per level and per direction, due to staggering) */ amrex::Vector,3 > > m_eb_update_E; + amrex::Vector,3 > > m_eb_update_B; /** EB: for every mesh face flag_info_face contains a: * * 0 if the face needs to be extended diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6b2c025c89f..cf709ca304c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -324,6 +324,7 @@ WarpX::WarpX () phi_dotMask.resize(nlevs_max); m_eb_update_E.resize(nlevs_max); + m_eb_update_B.resize(nlevs_max); m_flag_info_face.resize(nlevs_max); m_flag_ext_face.resize(nlevs_max); m_borrowing.resize(nlevs_max); @@ -2310,6 +2311,13 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); + AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[x]"); + AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[y]"); + AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[z]"); + // TODO: do not allocate these arrays anymore with the Yee solver //! EB: Lengths of the mesh edges m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), From ab87809509adca02439740e7f8c7d72a03406b73 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2025 09:12:40 -0800 Subject: [PATCH 13/32] Fix RZ compilation error --- Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 509027c5569..926f52aa8ee 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -295,7 +295,7 @@ void FiniteDifferenceSolver::EvolveECylindrical ( [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ // Skip field push in the embedded boundaries - if (update_Er_arr && update_Er_arr(i, j, k) == 0) { return; } + if (update_Er_arr && update_Er_arr(i, j, 0) == 0) { return; } Real const r = rmin + (i + 0.5_rt)*dr; // r on cell-centered point (Er is cell-centered in r) Er(i, j, 0, 0) += c2 * dt*( @@ -316,7 +316,7 @@ void FiniteDifferenceSolver::EvolveECylindrical ( [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ // Skip field push in the embedded boundaries - if (update_Et_arr && update_Et_arr(i, j, k) == 0) { return; } + if (update_Et_arr && update_Et_arr(i, j, 0) == 0) { return; } Real const r = rmin + i*dr; // r on a nodal grid (Et is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations @@ -361,7 +361,7 @@ void FiniteDifferenceSolver::EvolveECylindrical ( [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ // Skip field push in the embedded boundaries - if (update_Ez_arr && update_Ez_arr(i, j, k) == 0) { return; } + if (update_Ez_arr && update_Ez_arr(i, j, 0) == 0) { return; } Real const r = rmin + i*dr; // r on a nodal grid (Ez is nodal in r) if (r != 0) { // Off-axis, regular Maxwell equations From 78bb8c1bb20ad12c0e5dd8316af4448a9edf1f8e Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 09:52:12 -0800 Subject: [PATCH 14/32] Update initialization of external fields --- .../HybridPICModel/HybridPICModel.cpp | 10 +- Source/Initialization/WarpXInitData.cpp | 105 +++++------------- Source/WarpX.H | 12 +- 3 files changed, 38 insertions(+), 89 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index 20989cbeca9..c4d48fc0c26 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -226,9 +226,8 @@ void HybridPICModel::InitData () m_J_external[0], m_J_external[1], m_J_external[2], - lev, PatchType::fine, 'e', - warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev), - warpx.m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, + warpx.m_eb_update_E); } } @@ -244,9 +243,8 @@ void HybridPICModel::GetCurrentExternal () m_J_external[0], m_J_external[1], m_J_external[2], - lev, PatchType::fine, 'e', - warpx.m_fields.get_alldirs(FieldType::edge_lengths, lev), - warpx.m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, + warpx.m_eb_update_E); } } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 1b541f85d8e..57f1802e939 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -976,23 +976,21 @@ WarpX::InitLevelData (int lev, Real /*time*/) // 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)) { + + // TODO: raise error when EB is activated ComputeExternalFieldOnGridUsingParser( FieldType::Bfield_aux, m_p_ext_field_params->Bxfield_parser->compile<4>(), m_p_ext_field_params->Byfield_parser->compile<4>(), m_p_ext_field_params->Bzfield_parser->compile<4>(), - lev, PatchType::fine, 'f', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, m_eb_update_B); ComputeExternalFieldOnGridUsingParser( FieldType::Bfield_cp, m_p_ext_field_params->Bxfield_parser->compile<4>(), m_p_ext_field_params->Byfield_parser->compile<4>(), m_p_ext_field_params->Bzfield_parser->compile<4>(), - lev, PatchType::coarse, 'f', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_mr_levels_alldirs(FieldType::face_areas, max_level)[lev]); + lev, PatchType::coarse, m_eb_update_E); } // if the input string for the E-field is "parse_e_ext_grid_function", @@ -1023,18 +1021,14 @@ WarpX::InitLevelData (int lev, Real /*time*/) m_p_ext_field_params->Exfield_parser->compile<4>(), m_p_ext_field_params->Eyfield_parser->compile<4>(), m_p_ext_field_params->Ezfield_parser->compile<4>(), - lev, PatchType::fine, 'e', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, m_eb_update_E); ComputeExternalFieldOnGridUsingParser( FieldType::Efield_cp, m_p_ext_field_params->Exfield_parser->compile<4>(), m_p_ext_field_params->Eyfield_parser->compile<4>(), m_p_ext_field_params->Ezfield_parser->compile<4>(), - lev, PatchType::coarse, 'e', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::coarse, m_eb_update_E); #ifdef AMREX_USE_EB if (eb_enabled) { if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { @@ -1068,9 +1062,8 @@ void WarpX::ComputeExternalFieldOnGridUsingParser ( amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, [[maybe_unused]] const char topology, - std::optional const& edge_lengths, - std::optional const& face_areas) + int lev, PatchType patch_type, + amrex::Vector,3 > > const& eb_update_field) { auto t = gett_new(lev); @@ -1093,8 +1086,6 @@ void WarpX::ComputeExternalFieldOnGridUsingParser ( const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect(); const amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect(); - const bool eb_enabled = EB::enabled(); - for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const amrex::Box& tbx = mfi.tilebox( x_nodal_flag, mfx->nGrowVect() ); const amrex::Box& tby = mfi.tilebox( y_nodal_flag, mfy->nGrowVect() ); @@ -1104,44 +1095,19 @@ void WarpX::ComputeExternalFieldOnGridUsingParser ( auto const& mfyfab = mfy->array(mfi); auto const& mfzfab = mfz->array(mfi); - amrex::Array4 lx, ly, lz, Sx, Sy, Sz; - if (eb_enabled) { - if (edge_lengths.has_value()) { - const auto& edge_lengths_array = edge_lengths.value(); - lx = edge_lengths_array[0]->array(mfi); - ly = edge_lengths_array[1]->array(mfi); - lz = edge_lengths_array[2]->array(mfi); - } - if (face_areas.has_value()) { - const auto& face_areas_array = face_areas.value(); - Sx = face_areas_array[0]->array(mfi); - Sy = face_areas_array[1]->array(mfi); - Sz = face_areas_array[2]->array(mfi); - } - } - -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Dim3 lx_lo, lx_hi, lz_lo, lz_hi; -#endif - if (eb_enabled) { -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - lx_lo = amrex::lbound(lx); - lx_hi = amrex::ubound(lx); - lz_lo = amrex::lbound(lz); - lz_hi = amrex::ubound(lz); -#endif + amrex::Array4 update_fx_arr, update_fy_arr, update_fz_arr; + if (EB::enabled()) { + update_fx_arr = eb_update_field[lev][0]->array(mfi); + update_fy_arr = eb_update_field[lev][1]->array(mfi); + update_fz_arr = eb_update_field[lev][2]->array(mfi); } amrex::ParallelFor (tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { -#ifdef AMREX_USE_EB -#ifdef WARPX_DIM_3D - if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and Sx(i, j, k)<=0))) { return; } -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - //In XZ and RZ Ex is associated with a x-edge, while Bx is associated with a z-edge - if(lx && ((topology=='e' and lx(i, j, k)<=0) or (topology=='f' and lz(i, j, k)<=0))) { return; } -#endif -#endif + + // Do not set fields inside the embedded boundary + if (update_fx_arr && update_fx_arr(i,j,k) == 0) { return; } + // Shift required in the x-, y-, or z- position // depending on the index type of the multifab #if defined(WARPX_DIM_1D_Z) @@ -1167,19 +1133,10 @@ void WarpX::ComputeExternalFieldOnGridUsingParser ( mfxfab(i,j,k) = fx_parser(x,y,z,t); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { -#ifdef AMREX_USE_EB -#ifdef WARPX_DIM_3D - if(ly && ((topology=='e' and ly(i, j, k)<=0) or (topology=='f' and Sy(i, j, k)<=0))) { return; } -#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 && - ((topology=='e' and (lx(std::min(i , lx_hi.x), std::min(j , lx_hi.y), k)<=0 - || lx(std::max(i-1, lx_lo.x), std::min(j , lx_hi.y), k)<=0 - || lz(std::min(i , lz_hi.x), std::min(j , lz_hi.y), k)<=0 - || lz(std::min(i , lz_hi.x), std::max(j-1, lz_lo.y), k)<=0)) or - (topology=='f' and Sy(i,j,k)<=0))) { return; } -#endif -#endif + + // Do not set fields inside the embedded boundary + if (update_fy_arr && update_fy_arr(i,j,k) == 0) { return; } + #if defined(WARPX_DIM_1D_Z) const amrex::Real x = 0._rt; const amrex::Real y = 0._rt; @@ -1203,14 +1160,10 @@ void WarpX::ComputeExternalFieldOnGridUsingParser ( mfyfab(i,j,k) = fy_parser(x,y,z,t); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { -#ifdef AMREX_USE_EB -#ifdef WARPX_DIM_3D - if(lz && ((topology=='e' and lz(i, j, k)<=0) or (topology=='f' and Sz(i, j, k)<=0))) { return; } -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - //In XZ and RZ Ez is associated with a z-edge, while Bz is associated with a x-edge - if(lz && ((topology=='e' and lz(i, j, k)<=0) or (topology=='f' and lx(i, j, k)<=0))) { return; } -#endif -#endif + + // Do not set fields inside the embedded boundary + if (update_fz_arr && update_fz_arr(i,j,k) == 0) { return; } + #if defined(WARPX_DIM_1D_Z) const amrex::Real x = 0._rt; const amrex::Real y = 0._rt; @@ -1407,9 +1360,7 @@ WarpX::LoadExternalFields (int const lev) m_p_ext_field_params->Bxfield_parser->compile<4>(), m_p_ext_field_params->Byfield_parser->compile<4>(), m_p_ext_field_params->Bzfield_parser->compile<4>(), - lev, PatchType::fine, 'f', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, m_eb_update_B); } else if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { #if defined(WARPX_DIM_RZ) @@ -1432,9 +1383,7 @@ WarpX::LoadExternalFields (int const lev) m_p_ext_field_params->Exfield_parser->compile<4>(), m_p_ext_field_params->Eyfield_parser->compile<4>(), m_p_ext_field_params->Ezfield_parser->compile<4>(), - lev, PatchType::fine, 'e', - m_fields.get_alldirs(FieldType::edge_lengths, lev), - m_fields.get_alldirs(FieldType::face_areas, lev)); + lev, PatchType::fine, m_eb_update_E ); } else if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { #if defined(WARPX_DIM_RZ) diff --git a/Source/WarpX.H b/Source/WarpX.H index 5bd823ad8ba..8fb1eaf4408 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -937,9 +937,8 @@ public: amrex::ParserExecutor<4> const& fx_parser, amrex::ParserExecutor<4> const& fy_parser, amrex::ParserExecutor<4> const& fz_parser, - int lev, PatchType patch_type, [[maybe_unused]] char topology, - std::optional const& edge_lengths = std::nullopt, - std::optional const& face_areas = std::nullopt); + int lev, PatchType patch_type, + amrex::Vector,3 > > const& eb_update_field); /** * \brief Load field values from a user-specified openPMD file, @@ -1416,8 +1415,11 @@ private: /** EB: Flag to indicate whether a gridpoint is inside the embedded boundary and therefore * whether the E or B should not be updated. (One array per level and per direction, due to staggering) */ - amrex::Vector,3 > > m_eb_update_E; - amrex::Vector,3 > > m_eb_update_B; + // TODO: leave as private, add setters, getters + public: + amrex::Vector,3 > > m_eb_update_E; + amrex::Vector,3 > > m_eb_update_B; + private: /** EB: for every mesh face flag_info_face contains a: * * 0 if the face needs to be extended From df9ae82e0e0ce90360ec3c2e0aa6246e64eab3a9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 11:12:01 -0800 Subject: [PATCH 15/32] Fix setting of eb_update --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index d2a2c991165..c6cb4f7bceb 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -353,12 +353,12 @@ WarpX::MarkUpdateCells ( #endif // Loop over neighboring cells int eb_update = 1; - for (int ii = i_start; ii <= 1; ++ii) { - for (int jj = j_start; jj <= 1; ++jj) { - for (int kk = k_start; kk <= 1; ++kk) { + for (int i_shift = i_start; i_shift <= 0; ++i_shift) { + for (int j_shift = j_start; j_shift <= 0; ++j_shift) { + for (int k_shift = k_start; k_shift <= 0; ++k_shift) { // If one of the neighboring is either partially or fully covered // (i.e. if they are not regular cells), do not update field - if ( !flag(i+ii, j+jj, k+kk).isRegular() ) { + if ( !flag(i+i_shift, j+j_shift, k+k_shift).isRegular() ) { eb_update = 0; } } From 59579b4da85f58e7196fe61eb28135c628665722 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 12:26:28 -0800 Subject: [PATCH 16/32] Implement ECT condition for B --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 71 +++++++++++++++++++++++-- Source/Initialization/WarpXInitData.cpp | 27 ++++++---- Source/WarpX.H | 11 ++++ 3 files changed, 95 insertions(+), 14 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index c6cb4f7bceb..8fe0a9d2fa8 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -306,10 +306,13 @@ WarpX::MarkUpdateCells ( 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(); - auto const & eb_update_arr = eb_update[idim]->array(mfi); + amrex::Array4 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); @@ -367,8 +370,6 @@ WarpX::MarkUpdateCells ( eb_update_arr(i, j, k) = eb_update; }); - // TODO: Handle the case of the ECT solver - } } @@ -377,6 +378,70 @@ WarpX::MarkUpdateCells ( } +void +WarpX::MarkECTUpdateECells ( + std::array< std::unique_ptr,3> & eb_update_E, + ablastr::fields::VectorField const& edge_lengths ) +{ +} + +void +WarpX::MarkECTUpdateBCells ( + std::array< std::unique_ptr,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 const & eb_update_Bx_arr = eb_update_B[0]->array(mfi); + amrex::Array4 const & eb_update_By_arr = eb_update_B[1]->array(mfi); + amrex::Array4 const & eb_update_Bz_arr = eb_update_B[2]->array(mfi); + +#ifdef WARPX_DIM_3D + amrex::Array4 const & Sx_arr = face_areas[0]->array(mfi); + amrex::Array4 const & Sy_arr = face_areas[1]->array(mfi); + amrex::Array4 const & Sz_arr = face_areas[2]->array(mfi); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Array4 const & Sy_arr = face_areas[1]->array(mfi); + amrex::Array4 const & lx_arr = edge_lengths[0]->array(mfi); + amrex::Array4 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::MarkExtensionCells () { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 57f1802e939..64d675878c7 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1253,17 +1253,6 @@ void WarpX::InitializeEBGridData (int lev) auto const eb_fact = fieldEBFactory(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 ); - // Mark on which grid points B should be updated - MarkUpdateCells( - m_eb_update_B[lev], - m_fields.get_alldirs(FieldType::Bfield_fp, lev), - eb_fact ); - // TODO: move inside if condition for ECT auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev); ComputeEdgeLengths(edge_lengths_lev, eb_fact); @@ -1274,8 +1263,24 @@ void WarpX::InitializeEBGridData (int lev) ScaleAreas(face_areas_lev, CellSize(lev)); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + // Mark on which grid points E should be updated + MarkECTUpdateECells( m_eb_update_E[lev], edge_lengths_lev ); + // Mark on which grid points B should be updated + MarkECTUpdateBCells( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); + // Compute additional quantities required for the ECT solver MarkExtensionCells(); ComputeFaceExtensions(); + } else { + // Mark on which grid points E should be updated (stair-case approximation) + MarkUpdateCells( + m_eb_update_E[lev], + m_fields.get_alldirs(FieldType::Efield_fp, lev), + eb_fact ); + // Mark on which grid points B should be updated (stair-case approximation) + MarkUpdateCells( + m_eb_update_B[lev], + m_fields.get_alldirs(FieldType::Bfield_fp, lev), + eb_fact ); } } diff --git a/Source/WarpX.H b/Source/WarpX.H index 8fb1eaf4408..0dce1b75905 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1032,6 +1032,17 @@ public: ablastr::fields::VectorField const & field, amrex::EBFArrayBoxFactory const & eb_fact ); + // TODO documentation + void MarkECTUpdateECells ( + std::array< std::unique_ptr,3> & eb_update_B, + ablastr::fields::VectorField const& edge_lengths ); + + // TODO documentation + void MarkECTUpdateBCells ( + std::array< std::unique_ptr,3> & eb_update_B, + ablastr::fields::VectorField const& face_areas, + ablastr::fields::VectorField const& edge_lengths ); + /** * \brief Compute the length of the mesh edges. Here the length is a value in [0, 1]. * An edge of length 0 is fully covered. From 4fd4dc7a0eddcd7f5811d539c2ad9475e93c2f81 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 12:36:32 -0800 Subject: [PATCH 17/32] Do not use any guard cells --- Source/WarpX.cpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index cf709ca304c..1df4bcde2e2 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2304,19 +2304,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (lev == maxLevel()) { if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { - AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]"); - AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]"); - AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); - - AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[x]"); - AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[y]"); - AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, - guard_cells.ng_FieldSolver, lev, "m_eb_update_B[z]"); + AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[x]"); + AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[y]"); + AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[z]"); + + AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[x]"); + AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[y]"); + AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, 0, + amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[z]"); // TODO: do not allocate these arrays anymore with the Yee solver //! EB: Lengths of the mesh edges From 4232db7040633d3a3bd3d671b2845c1bfc48fa82 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 15:53:11 -0800 Subject: [PATCH 18/32] Fixed ECT solver --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 54 +++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 8fe0a9d2fa8..559da0472a2 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -383,6 +383,60 @@ WarpX::MarkECTUpdateECells ( std::array< std::unique_ptr,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 const & eb_update_Ex_arr = eb_update_E[0]->array(mfi); + amrex::Array4 const & eb_update_Ey_arr = eb_update_E[1]->array(mfi); + amrex::Array4 const & eb_update_Ez_arr = eb_update_E[2]->array(mfi); + + amrex::Array4 const & lx_arr = edge_lengths[0]->array(mfi); + amrex::Array4 const & lz_arr = edge_lengths[2]->array(mfi); +#if defined(WARPX_DIM_3D) + amrex::Array4 const & ly_arr = edge_lengths[1]->array(mfi); +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + amrex::Dim3 lx_lo = amrex::lbound(lx_arr); + amrex::Dim3 lx_hi = amrex::ubound(lx_arr); + amrex::Dim3 lz_lo = amrex::lbound(lz_arr); + amrex::Dim3 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 From f45b257f2c6f5bf8d6ca9c15c1bcd971cceb4563 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 16:00:04 -0800 Subject: [PATCH 19/32] Revert "Do not use any guard cells" This reverts commit 4fd4dc7a0eddcd7f5811d539c2ad9475e93c2f81. --- Source/WarpX.cpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1df4bcde2e2..cf709ca304c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2304,19 +2304,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (lev == maxLevel()) { if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { - AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[x]"); - AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[y]"); - AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_E[z]"); - - AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[x]"); - AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[y]"); - AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, 0, - amrex::IntVect::TheZeroVector(), lev, "m_eb_update_B[z]"); + AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]"); + AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]"); + AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); + + AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[x]"); + AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[y]"); + AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, + guard_cells.ng_FieldSolver, lev, "m_eb_update_B[z]"); // TODO: do not allocate these arrays anymore with the Yee solver //! EB: Lengths of the mesh edges From 2003d8593385b504f76ac5c6e619ccef9b4b5edf Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 18:49:37 -0800 Subject: [PATCH 20/32] Fix compilation errors --- .../embedded_boundary_cube/inputs_base_3d | 4 +- .../inputs_test_2d_embedded_boundary_cube | 8 ++-- Source/EmbeddedBoundary/WarpXInitEB.cpp | 1 + .../FiniteDifferenceSolver.H | 4 +- .../MacroscopicEvolveE.cpp | 44 +++++++++---------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- 6 files changed, 30 insertions(+), 33 deletions(-) diff --git a/Examples/Tests/embedded_boundary_cube/inputs_base_3d b/Examples/Tests/embedded_boundary_cube/inputs_base_3d index 9710701d871..70ddd8f8f64 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_base_3d +++ b/Examples/Tests/embedded_boundary_cube/inputs_base_3d @@ -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.501 -0.501 -0.501 # Ensures that the stair-case EB is exactly at -0.5 +eb2.box_hi = 0.501 0.501 0.501 # 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. diff --git a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube index 4fbf13e7f5c..46272052c2c 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube +++ b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube @@ -12,10 +12,10 @@ warpx.abort_on_warning_threshold = medium boundary.field_lo = pec pec boundary.field_hi = pec pec -my_constants.xmin = -0.501 -my_constants.zmin = -0.501 -my_constants.xmax = 0.501 -my_constants.zmax = 0.501 +my_constants.xmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5 +my_constants.zmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5 +my_constants.xmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5 +my_constants.zmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5 # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular. warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))" diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 559da0472a2..fca6cc3114b 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -463,6 +463,7 @@ WarpX::MarkECTUpdateBCells ( amrex::Array4 const & Sx_arr = face_areas[0]->array(mfi); amrex::Array4 const & Sy_arr = face_areas[1]->array(mfi); amrex::Array4 const & Sz_arr = face_areas[2]->array(mfi); + amrex::ignore_unused(edge_lengths); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::Array4 const & Sy_arr = face_areas[1]->array(mfi); amrex::Array4 const & lx_arr = edge_lengths[0]->array(mfi); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 07596dc8562..7726a2ed5bd 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -111,7 +111,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real dt, std::unique_ptr const& macroscopic_properties); @@ -313,7 +313,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real dt, std::unique_ptr const& macroscopic_properties); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp index 708728c4e5b..20dab1b2d12 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -40,7 +40,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveE ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real const dt, std::unique_ptr const& macroscopic_properties) { @@ -48,7 +48,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveE ( // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ - amrex::ignore_unused(Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties); + amrex::ignore_unused(Efield, Bfield, Jfield, eb_update_E, dt, macroscopic_properties); WARPX_ABORT_WITH_MESSAGE("currently macro E-push does not work for RZ"); #else @@ -61,13 +61,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveE ( if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) { MacroscopicEvolveECartesian - ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties); + ( Efield, Bfield, Jfield, eb_update_E, dt, macroscopic_properties); } if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) { MacroscopicEvolveECartesian - ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties); + ( Efield, Bfield, Jfield, eb_update_E, dt, macroscopic_properties); } @@ -78,12 +78,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveE ( if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) { MacroscopicEvolveECartesian - ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties); + ( Efield, Bfield, Jfield, eb_update_E, dt, macroscopic_properties); } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) { MacroscopicEvolveECartesian - ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties); + ( Efield, Bfield, Jfield, eb_update_E, dt, macroscopic_properties); } @@ -103,7 +103,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Bfield, ablastr::fields::VectorField const& Jfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, amrex::Real const dt, std::unique_ptr const& macroscopic_properties) { @@ -141,15 +141,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( Array4 const& jy = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); - amrex::Array4 eb_lx, eb_ly, eb_lz; + amrex::Array4 update_Ex_arr, update_Ey_arr, update_Ez_arr; if (EB::enabled()) { - eb_lx = edge_lengths[0]->array(mfi); - eb_ly = edge_lengths[1]->array(mfi); - eb_lz = edge_lengths[2]->array(mfi); + update_Ex_arr = eb_update_E[0]->array(mfi); + update_Ey_arr = eb_update_E[1]->array(mfi); + update_Ez_arr = eb_update_E[2]->array(mfi); } -#ifdef WARPX_DIM_XZ - amrex::ignore_unused(eb_ly); -#endif // material prop // amrex::Array4 const& sigma_arr = sigma_mf.array(mfi); @@ -180,8 +177,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( // Loop over the cells and update the fields amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip field push if this cell is fully covered by embedded boundaries - if (eb_lx && eb_lx(i, j, k) <= 0) { return; } + + // Skip field push in the embedded boundaries + if (update_Ex_arr && update_Ex_arr(i, j, k) == 0) { return; } // Interpolate conductivity, sigma, to Ex position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, @@ -198,12 +196,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ -#ifdef WARPX_DIM_3D - if (eb_ly && eb_ly(i,j,k) <= 0) { return; } -#elif defined(WARPX_DIM_XZ) - //In XZ Ey is associated with a mesh node, so we need to check if the mesh node is covered - if (eb_lx && (eb_lx(i, j, k)<=0 || eb_lx(i-1, j, k)<=0 || eb_lz(i, j, k)<=0 || eb_lz(i, j-1, k)<=0)) { return; } -#endif + + // Skip field push in the embedded boundaries + if (update_Ey_arr && update_Ey_arr(i, j, k) == 0) { return; } // Interpolate conductivity, sigma, to Ey position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, @@ -221,8 +216,9 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip field push if this cell is fully covered by embedded boundaries - if (eb_lz && eb_lz(i, j, k) <= 0) { return; } + + // Skip field push in the embedded boundaries + if (update_Ez_arr && update_Ez_arr(i, j, k) == 0) { return; } // Interpolate conductivity, sigma, to Ez position on the grid amrex::Real const sigma_interp = ablastr::coarsen::sample::Interp(sigma_arr, sigma_stag, diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index a77fb390fff..139988b69a2 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -1157,7 +1157,7 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt, amre m_fields.get_alldirs(FieldType::Efield_fp, lev), m_fields.get_alldirs(FieldType::Bfield_fp, lev), m_fields.get_alldirs(FieldType::current_fp, lev), - m_fields.get_alldirs(FieldType::edge_lengths, lev), + m_eb_update_E[lev], a_dt, m_macroscopic_properties); if (do_pml && pml[lev]->ok()) { From d6e64af21c125a31c8b65e850a95c6ea71cbc742 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 11 Jan 2025 20:34:59 -0800 Subject: [PATCH 21/32] Fix compilation error --- .../FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp index 20dab1b2d12..33d368925f7 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -107,10 +107,6 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian ( amrex::Real const dt, std::unique_ptr const& macroscopic_properties) { -#ifndef AMREX_USE_EB - amrex::ignore_unused(edge_lengths); -#endif - amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf(); amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf(); amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf(); From e6d93f9f9be813e2a743819145a1f0c195a1aebd Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 12 Jan 2025 07:10:05 -0800 Subject: [PATCH 22/32] Fix const-ness --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index fca6cc3114b..3cc29d17b39 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -343,16 +343,17 @@ WarpX::MarkUpdateCells ( // or fully covered: do not update field // Check neighbors in the each spatial direction - int i_start = 0; - int j_start = 0; - int k_start = 0; // If nodal in a given direction, we need to start from -1 (left-neighboring cell) - if ( index_type.nodeCentered(0) ) { i_start = -1; }; + int const i_start = ( index_type.nodeCentered(0) )? -1 : 0; #if AMREX_SPACEDIM > 1 - if ( index_type.nodeCentered(1) ) { j_start = -1; }; + int const j_start = ( index_type.nodeCentered(1) )? -1 : 0; +#else + int const j_start = 0; #endif #if AMREX_SPACEDIM > 2 - if ( index_type.nodeCentered(2) ) { k_start = -1; }; + int const k_start = ( index_type.nodeCentered(2) )? -1 : 0; +#else + int const k_start = 0; #endif // Loop over neighboring cells int eb_update = 1; @@ -402,10 +403,10 @@ WarpX::MarkECTUpdateECells ( #if defined(WARPX_DIM_3D) amrex::Array4 const & ly_arr = edge_lengths[1]->array(mfi); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Dim3 lx_lo = amrex::lbound(lx_arr); - amrex::Dim3 lx_hi = amrex::ubound(lx_arr); - amrex::Dim3 lz_lo = amrex::lbound(lz_arr); - amrex::Dim3 lz_hi = amrex::ubound(lz_arr); + 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, From eb4308775f1e8dd66ea5d9e2db993dcd21aa0ef8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 12 Jan 2025 08:28:15 -0800 Subject: [PATCH 23/32] Activate load-balancing --- Source/Parallelization/WarpXRegrid.cpp | 10 ++++++++++ Source/WarpX.H | 4 ++-- Source/WarpX.cpp | 4 ++-- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index a0a2d4929df..c12b9f7e0f7 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -174,6 +174,14 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi using ablastr::fields::Direction; using warpx::fields::FieldType; + const auto RemakeMultiFab = [&](auto& mf){ + if (mf == nullptr) { return; } + const IntVect& ng = mf->nGrowVect(); + auto pmf = std::remove_reference_t{}; + AllocInitMultiFab(pmf, mf->boxArray(), dm, mf->nComp(), ng, lev, mf->tags()[0]); + mf = std::move(pmf); + }; + bool const eb_enabled = EB::enabled(); if (ba == boxArray(lev)) { @@ -187,6 +195,8 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi { if (eb_enabled) { if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { + RemakeMultiFab( m_eb_update_E[lev][idim] ); + RemakeMultiFab( m_eb_update_B[lev][idim] ); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { m_borrowing[lev][idim] = std::make_unique>(amrex::convert(ba, Bfield_fp[lev][idim]->ixType().toIntVect()), dm); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 0dce1b75905..5e4650e3745 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1015,7 +1015,7 @@ public: /** 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` + * More specifically, this function fills 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. @@ -1034,7 +1034,7 @@ public: // TODO documentation void MarkECTUpdateECells ( - std::array< std::unique_ptr,3> & eb_update_B, + std::array< std::unique_ptr,3> & eb_update_E, ablastr::fields::VectorField const& edge_lengths ); // TODO documentation diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index cf709ca304c..265603eae24 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2304,14 +2304,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (lev == maxLevel()) { if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { - AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, + AllocInitMultiFab(m_eb_update_E[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_E[x]"); AllocInitMultiFab(m_eb_update_E[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_E[y]"); AllocInitMultiFab(m_eb_update_E[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_E[z]"); - AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, + AllocInitMultiFab(m_eb_update_B[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_B[x]"); AllocInitMultiFab(m_eb_update_B[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_B[y]"); From 2eb3383a1d26dbb6dd7a2a756e4227dd3e5248fa Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 12 Jan 2025 11:10:17 -0800 Subject: [PATCH 24/32] Update checksums --- .../test_2d_embedded_boundary_cube.json | 4 ++-- .../benchmarks_json/test_2d_field_probe.json | 8 ++++---- .../benchmarks_json/test_3d_eb_picmi.json | 14 +++++++------- .../test_3d_embedded_boundary_cube.json | 10 +++++----- ...est_3d_embedded_boundary_cube_macroscopic.json | 10 +++++----- .../test_3d_embedded_boundary_rotated_cube.json | 15 +++++++-------- .../test_3d_particle_absorption.json | 14 +++++++------- .../benchmarks_json/test_3d_particle_scrape.json | 14 +++++++------- .../test_3d_particle_scrape_picmi.json | 14 +++++++------- .../test_rz_embedded_boundary_diffraction.json | 12 ++++++------ 10 files changed, 57 insertions(+), 58 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json index a3e609bd9a9..dbb5ffa39ae 100644 --- a/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json @@ -3,8 +3,8 @@ "Bx": 9.263694545408503e-05, "By": 0.00031905198933489145, "Bz": 7.328424783762594e-05, - "Ex": 8553.906698053046, + "Ex": 8553.90669811286, "Ey": 60867.04830538045, - "Ez": 8.439422682267567e-07 + "Ez": 4.223902107031194e-06 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_2d_field_probe.json b/Regression/Checksum/benchmarks_json/test_2d_field_probe.json index cb82acfc067..8aabe6c8301 100644 --- a/Regression/Checksum/benchmarks_json/test_2d_field_probe.json +++ b/Regression/Checksum/benchmarks_json/test_2d_field_probe.json @@ -1,10 +1,10 @@ { "lev=0": { "Bx": 0.0, - "By": 126826.78487921853, + "By": 123510.69657444415, "Bz": 0.0, - "Ex": 32517064310550.266, + "Ex": 31206368949280.34, "Ey": 0.0, - "Ez": 17321323003697.61 + "Ez": 16921005306450.537 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json index ad0d2cee5a3..1f9f0a77b5a 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 148673.005859208, - "By": 148673.00585920806, - "Bz": 3371.758117878558, - "Ex": 55378581103426.71, - "Ey": 55378581103426.72, - "Ez": 68412803445328.25 + "Bx": 144495.08082507108, + "By": 144495.08082507114, + "Bz": 8481.958724628861, + "Ex": 54500496182517.92, + "Ey": 54500496182517.91, + "Ez": 70231240245509.39 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json index 58ee8806540..9563c52adbe 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 4.060477854092961e-18, + "Bx": 4.166971025838921e-18, "By": 0.006628374119786834, "Bz": 0.006628374119786834, - "Ex": 5102618.4711524295, - "Ey": 6.323754160591239e-05, - "Ez": 6.323754160591239e-05 + "Ex": 5102618.471153786, + "Ey": 1.4283859321773714e-05, + "Ez": 1.4283859321773714e-05 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json index 8cc6af7cb93..67bdbea18ca 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 4.20930075273562e-18, + "Bx": 4.228863291892693e-18, "By": 0.005101824310293573, "Bz": 0.005101824310293573, - "Ex": 4414725.184731115, - "Ey": 6.32375413967707e-05, - "Ez": 6.32375413967707e-05 + "Ex": 4414725.184732471, + "Ey": 1.4283895626502055e-05, + "Ez": 1.4283895626502055e-05 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json index b2b4aa569c1..d45887ca932 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json @@ -1,11 +1,10 @@ { "lev=0": { - "Bx": 1.280747509243305e-05, - "By": 2.473900144296397e-02, - "Bz": 2.473890786894079e-02, - "Ex": 1.025322901921306e+07, - "Ey": 1.042254197269831e+04, - "Ez": 1.040011664019071e+04 + "Bx": 1.252616939910365e-05, + "By": 0.02473895628331097, + "Bz": 0.024738956316621142, + "Ex": 10253221.850298548, + "Ey": 10387.334582977643, + "Ez": 10387.532806510022 } -} - +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json b/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json index ce6e2fcf79b..3dc9d956b79 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 202106.71291347666, - "By": 202106.71291347663, - "Bz": 3371.897999274175, - "Ex": 38304043178806.11, - "Ey": 38304043178806.11, - "Ez": 83057027925874.84 + "Bx": 198610.0530604908, + "By": 198610.0530604909, + "Bz": 8482.656173586969, + "Ex": 37232105734622.53, + "Ey": 37232105734622.54, + "Ez": 85094015810307.19 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json index b03a954397a..9437ebed275 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 148673.005859208, - "By": 148673.00585920803, - "Bz": 3371.758117878557, - "Ex": 55378581103426.695, - "Ey": 55378581103426.7, - "Ez": 68412803445328.25 + "Bx": 144495.08082507108, + "By": 144495.0808250711, + "Bz": 8481.95872462886, + "Ex": 54500496182517.914, + "Ey": 54500496182517.914, + "Ez": 70231240245509.4 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json index b03a954397a..1f9f0a77b5a 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 148673.005859208, - "By": 148673.00585920803, - "Bz": 3371.758117878557, - "Ex": 55378581103426.695, - "Ey": 55378581103426.7, - "Ez": 68412803445328.25 + "Bx": 144495.08082507108, + "By": 144495.08082507114, + "Bz": 8481.958724628861, + "Ex": 54500496182517.92, + "Ey": 54500496182517.91, + "Ez": 70231240245509.39 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json b/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json index 0e5fad8db8a..e4b9d9c07ff 100644 --- a/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json +++ b/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json @@ -1,10 +1,10 @@ { "lev=0": { - "Br": 6.821267675779345e-19, - "Bt": 5.564905732478707e-05, - "Bz": 2.368259586613272e-19, - "Er": 16503.98082446463, - "Et": 1.5299584682447838e-10, - "Ez": 1466.854467399728 + "Br": 6.7914286131989935e-19, + "Bt": 5.4557350206853276e-05, + "Bz": 2.357229221622199e-19, + "Er": 16481.39008058988, + "Et": 1.5258937379236053e-10, + "Ez": 1508.1064116028576 } } \ No newline at end of file From 8981179a16ccd0c80eb61bbd7ba107cd9f7364e2 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 12 Jan 2025 14:26:49 -0800 Subject: [PATCH 25/32] Add documentation --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 6 +- .../HybridPICModel/HybridPICModel.cpp | 4 +- Source/Initialization/WarpXInitData.cpp | 8 +-- Source/WarpX.H | 58 ++++++++++++------- 4 files changed, 45 insertions(+), 31 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 3cc29d17b39..d7eb98489dc 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -292,7 +292,7 @@ WarpX::ScaleAreas (ablastr::fields::VectorField& face_areas, } void -WarpX::MarkUpdateCells ( +WarpX::MarkUpdateCellsStairCase ( std::array< std::unique_ptr,3> & eb_update, ablastr::fields::VectorField const& field, amrex::EBFArrayBoxFactory const & eb_fact ) @@ -380,7 +380,7 @@ WarpX::MarkUpdateCells ( } void -WarpX::MarkECTUpdateECells ( +WarpX::MarkUpdateECellsECT ( std::array< std::unique_ptr,3> & eb_update_E, ablastr::fields::VectorField const& edge_lengths ) { @@ -441,7 +441,7 @@ WarpX::MarkECTUpdateECells ( } void -WarpX::MarkECTUpdateBCells ( +WarpX::MarkUpdateBCellsECT ( std::array< std::unique_ptr,3> & eb_update_B, ablastr::fields::VectorField const& face_areas, ablastr::fields::VectorField const& edge_lengths ) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index c4d48fc0c26..abda59e40ba 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -227,7 +227,7 @@ void HybridPICModel::InitData () m_J_external[1], m_J_external[2], lev, PatchType::fine, - warpx.m_eb_update_E); + warpx.GetEBUpdateEFlag()); } } @@ -244,7 +244,7 @@ void HybridPICModel::GetCurrentExternal () m_J_external[1], m_J_external[2], lev, PatchType::fine, - warpx.m_eb_update_E); + warpx.GetEBUpdateEFlag()); } } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 64d675878c7..c8c7d1f0337 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1264,20 +1264,20 @@ void WarpX::InitializeEBGridData (int lev) if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { // Mark on which grid points E should be updated - MarkECTUpdateECells( m_eb_update_E[lev], edge_lengths_lev ); + MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev ); // Mark on which grid points B should be updated - MarkECTUpdateBCells( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); + MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); // Compute additional quantities required for the ECT solver MarkExtensionCells(); ComputeFaceExtensions(); } else { // Mark on which grid points E should be updated (stair-case approximation) - MarkUpdateCells( + MarkUpdateCellsStairCase( m_eb_update_E[lev], m_fields.get_alldirs(FieldType::Efield_fp, lev), eb_fact ); // Mark on which grid points B should be updated (stair-case approximation) - MarkUpdateCells( + MarkUpdateCellsStairCase( m_eb_update_B[lev], m_fields.get_alldirs(FieldType::Bfield_fp, lev), eb_fact ); diff --git a/Source/WarpX.H b/Source/WarpX.H index 5e4650e3745..ad28c672ad5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -169,6 +169,7 @@ public: } #endif ParticleBoundaryBuffer& GetParticleBoundaryBuffer () { return *m_particle_boundary_buffer; } + amrex::Vector,3 > >& GetEBUpdateEFlag() { return m_eb_update_E; } static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, int dir, int lev, bool update_cost_flag, @@ -1012,33 +1013,49 @@ public: void InitEB (); #ifdef AMREX_USE_EB - /** Set a flag to indicate on which grid points the field `field` - * should be updated. + /** \brief Set a flag to indicate on which grid points the field `field` + * should be updated, depending on their position relative to the embedded boundary. + * + * This function is used by all finite-difference solvers, except the + * ECT solver, which instead uses `MarkUpdateECellsECT` and `MarkUpdateBCellsECT`. + * It uses a stair-case approximation of the embedded boundary: + * If a grid point touches cells that are either partially or fully covered + * by the embedded boundary: the corresponding field is not updated. * * More specifically, this function fills 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 + * or 0, depending on whether the grid point should be updated or not. */ - void MarkUpdateCells ( + void MarkUpdateCellsStairCase ( std::array< std::unique_ptr,3> & eb_update, ablastr::fields::VectorField const & field, amrex::EBFArrayBoxFactory const & eb_fact ); - // TODO documentation - void MarkECTUpdateECells ( + /** \brief Set a flag to indicate on which grid points the E field + * should be updated, depending on their position relative to the embedded boundary. + * + * This function is used by ECT solver. The E field is not updated if + * the edge on which it is defined is fully covered by the embedded boundary. + * + * More specifically, this function fills the iMultiFabs in `eb_update_E` + * (which have the same indexType as the E field) with 1 or 0, depending + * on whether the grid point should be updated or not. + */ + void MarkUpdateECellsECT ( std::array< std::unique_ptr,3> & eb_update_E, ablastr::fields::VectorField const& edge_lengths ); - // TODO documentation - void MarkECTUpdateBCells ( + /** \brief Set a flag to indicate on which grid points the B field + * should be updated, depending on their position relative to the embedded boundary. + * + * This function is used by ECT solver. The B field is not updated if + * the face on which it is defined is fully covered by the embedded boundary. + * + * More specifically, this function fills the iMultiFabs in `eb_update_B` + * (which have the same indexType as the B field) with 1 or 0, depending + * on whether the grid point should be updated or not. + */ + void MarkUpdateBCellsECT ( std::array< std::unique_ptr,3> & eb_update_B, ablastr::fields::VectorField const& face_areas, ablastr::fields::VectorField const& edge_lengths ); @@ -1425,12 +1442,9 @@ private: /** EB: Flag to indicate whether a gridpoint is inside the embedded boundary and therefore * whether the E or B should not be updated. (One array per level and per direction, due to staggering) - */ - // TODO: leave as private, add setters, getters - public: - amrex::Vector,3 > > m_eb_update_E; - amrex::Vector,3 > > m_eb_update_B; - private: + */ + amrex::Vector,3 > > m_eb_update_E; + amrex::Vector,3 > > m_eb_update_B; /** EB: for every mesh face flag_info_face contains a: * * 0 if the face needs to be extended From 8182df74a86d528604ed3c2b461bacedbfdb25a5 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 14 Jan 2025 12:52:35 -0800 Subject: [PATCH 26/32] Apply suggestions from code review --- Source/Initialization/WarpXInitData.cpp | 1 - Source/WarpX.cpp | 1 - 2 files changed, 2 deletions(-) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index c8c7d1f0337..5575e6211f9 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1253,7 +1253,6 @@ void WarpX::InitializeEBGridData (int lev) auto const eb_fact = fieldEBFactory(lev); - // TODO: move inside if condition for ECT auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev); ComputeEdgeLengths(edge_lengths_lev, eb_fact); ScaleEdges(edge_lengths_lev, CellSize(lev)); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 265603eae24..db9abb851c1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2318,7 +2318,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(m_eb_update_B[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_eb_update_B[z]"); - // TODO: do not allocate these arrays anymore with the Yee solver //! EB: Lengths of the mesh edges m_fields.alloc_init(FieldType::edge_lengths, Direction{0}, lev, amrex::convert(ba, Ex_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); From 0980516a26ab1ccdcf89fd3a539ce8eb0452ca8a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 17 Jan 2025 08:19:50 -0800 Subject: [PATCH 27/32] Minor cleanup --- Source/Initialization/WarpXInitData.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 5575e6211f9..4db7096eb05 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -977,7 +977,6 @@ WarpX::InitLevelData (int lev, Real /*time*/) if ((m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::parse_ext_grid_function) && (lev > 0) && (lev <= maxlevel_extEMfield_init)) { - // TODO: raise error when EB is activated ComputeExternalFieldOnGridUsingParser( FieldType::Bfield_aux, m_p_ext_field_params->Bxfield_parser->compile<4>(), @@ -990,7 +989,7 @@ WarpX::InitLevelData (int lev, Real /*time*/) m_p_ext_field_params->Bxfield_parser->compile<4>(), m_p_ext_field_params->Byfield_parser->compile<4>(), m_p_ext_field_params->Bzfield_parser->compile<4>(), - lev, PatchType::coarse, m_eb_update_E); + lev, PatchType::coarse, m_eb_update_B); } // if the input string for the E-field is "parse_e_ext_grid_function", From ace399b6a3d044436563165a27bc93cb25b190b3 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 18 Jan 2025 08:34:54 -0800 Subject: [PATCH 28/32] Modify initialization in ECT solver --- .../test_3d_embedded_boundary_rotated_cube.json | 12 ++++++------ Source/Initialization/WarpXInitData.cpp | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json index d45887ca932..118214948a5 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 1.252616939910365e-05, - "By": 0.02473895628331097, - "Bz": 0.024738956316621142, - "Ex": 10253221.850298548, - "Ey": 10387.334582977643, - "Ez": 10387.532806510022 + "Bx": 1.280747509243305e-05, + "By": 2.473900144296397e-02, + "Bz": 2.473890786894079e-02, + "Ex": 1.025322901921306e+07, + "Ey": 1.042254197269831e+04, + "Ez": 1.040011664019071e+04 } } \ No newline at end of file diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 4db7096eb05..248448e817e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1261,13 +1261,13 @@ void WarpX::InitializeEBGridData (int lev) ScaleAreas(face_areas_lev, CellSize(lev)); if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + // Compute additional quantities required for the ECT solver + MarkExtensionCells(); + ComputeFaceExtensions(); // Mark on which grid points E should be updated MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev ); // Mark on which grid points B should be updated MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); - // Compute additional quantities required for the ECT solver - MarkExtensionCells(); - ComputeFaceExtensions(); } else { // Mark on which grid points E should be updated (stair-case approximation) MarkUpdateCellsStairCase( From 16ec9b2867d754e292ba5cf3e8e36d7da2782496 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 18 Jan 2025 08:53:19 -0800 Subject: [PATCH 29/32] Update comments --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 37 +++++++++++++++---------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index d7eb98489dc..0ddd8c22d1f 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -339,30 +339,39 @@ WarpX::MarkUpdateCellsStairCase ( amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - // Stair-case approximation: If neighboring cells are either partially - // or fully covered: do not update field - - // Check neighbors in the each spatial direction - // If nodal in a given direction, we need to start from -1 (left-neighboring cell) - int const i_start = ( index_type.nodeCentered(0) )? -1 : 0; + // 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) )? -1 : 0; + int const j_start = ( index_type.nodeCentered(1) )? j-1 : j; #else - int const j_start = 0; + int const j_start = j; #endif #if AMREX_SPACEDIM > 2 - int const k_start = ( index_type.nodeCentered(2) )? -1 : 0; + int const k_start = ( index_type.nodeCentered(2) )? k-1 : k; #else - int const k_start = 0; + int const k_start = k; #endif // Loop over neighboring cells int eb_update = 1; - for (int i_shift = i_start; i_shift <= 0; ++i_shift) { - for (int j_shift = j_start; j_shift <= 0; ++j_shift) { - for (int k_shift = k_start; k_shift <= 0; ++k_shift) { + 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 - if ( !flag(i+i_shift, j+j_shift, k+k_shift).isRegular() ) { + // (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; } } From f2b0580026588d827f1fd20d8c00c2402df81687 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 17 Jan 2025 06:15:00 -0800 Subject: [PATCH 30/32] Revert "Update checksums" This reverts commit 2eb3383a1d26dbb6dd7a2a756e4227dd3e5248fa. --- .../test_2d_embedded_boundary_cube.json | 4 ++-- .../benchmarks_json/test_2d_field_probe.json | 8 ++++---- .../Checksum/benchmarks_json/test_3d_eb_picmi.json | 14 +++++++------- .../test_3d_embedded_boundary_cube.json | 10 +++++----- ...test_3d_embedded_boundary_cube_macroscopic.json | 10 +++++----- .../test_3d_embedded_boundary_rotated_cube.json | 3 ++- .../test_3d_particle_absorption.json | 14 +++++++------- .../benchmarks_json/test_3d_particle_scrape.json | 14 +++++++------- .../test_3d_particle_scrape_picmi.json | 14 +++++++------- .../test_rz_embedded_boundary_diffraction.json | 12 ++++++------ 10 files changed, 52 insertions(+), 51 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json index dbb5ffa39ae..a3e609bd9a9 100644 --- a/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json @@ -3,8 +3,8 @@ "Bx": 9.263694545408503e-05, "By": 0.00031905198933489145, "Bz": 7.328424783762594e-05, - "Ex": 8553.90669811286, + "Ex": 8553.906698053046, "Ey": 60867.04830538045, - "Ez": 4.223902107031194e-06 + "Ez": 8.439422682267567e-07 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_2d_field_probe.json b/Regression/Checksum/benchmarks_json/test_2d_field_probe.json index 8aabe6c8301..cb82acfc067 100644 --- a/Regression/Checksum/benchmarks_json/test_2d_field_probe.json +++ b/Regression/Checksum/benchmarks_json/test_2d_field_probe.json @@ -1,10 +1,10 @@ { "lev=0": { "Bx": 0.0, - "By": 123510.69657444415, + "By": 126826.78487921853, "Bz": 0.0, - "Ex": 31206368949280.34, + "Ex": 32517064310550.266, "Ey": 0.0, - "Ez": 16921005306450.537 + "Ez": 17321323003697.61 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json index 1f9f0a77b5a..ad0d2cee5a3 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 144495.08082507108, - "By": 144495.08082507114, - "Bz": 8481.958724628861, - "Ex": 54500496182517.92, - "Ey": 54500496182517.91, - "Ez": 70231240245509.39 + "Bx": 148673.005859208, + "By": 148673.00585920806, + "Bz": 3371.758117878558, + "Ex": 55378581103426.71, + "Ey": 55378581103426.72, + "Ez": 68412803445328.25 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json index 9563c52adbe..58ee8806540 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 4.166971025838921e-18, + "Bx": 4.060477854092961e-18, "By": 0.006628374119786834, "Bz": 0.006628374119786834, - "Ex": 5102618.471153786, - "Ey": 1.4283859321773714e-05, - "Ez": 1.4283859321773714e-05 + "Ex": 5102618.4711524295, + "Ey": 6.323754160591239e-05, + "Ez": 6.323754160591239e-05 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json index 67bdbea18ca..8cc6af7cb93 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 4.228863291892693e-18, + "Bx": 4.20930075273562e-18, "By": 0.005101824310293573, "Bz": 0.005101824310293573, - "Ex": 4414725.184732471, - "Ey": 1.4283895626502055e-05, - "Ez": 1.4283895626502055e-05 + "Ex": 4414725.184731115, + "Ey": 6.32375413967707e-05, + "Ez": 6.32375413967707e-05 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json index 118214948a5..b2b4aa569c1 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json +++ b/Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json @@ -7,4 +7,5 @@ "Ey": 1.042254197269831e+04, "Ez": 1.040011664019071e+04 } -} \ No newline at end of file +} + diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json b/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json index 3dc9d956b79..ce6e2fcf79b 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 198610.0530604908, - "By": 198610.0530604909, - "Bz": 8482.656173586969, - "Ex": 37232105734622.53, - "Ey": 37232105734622.54, - "Ez": 85094015810307.19 + "Bx": 202106.71291347666, + "By": 202106.71291347663, + "Bz": 3371.897999274175, + "Ex": 38304043178806.11, + "Ey": 38304043178806.11, + "Ez": 83057027925874.84 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json index 9437ebed275..b03a954397a 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 144495.08082507108, - "By": 144495.0808250711, - "Bz": 8481.95872462886, - "Ex": 54500496182517.914, - "Ey": 54500496182517.914, - "Ez": 70231240245509.4 + "Bx": 148673.005859208, + "By": 148673.00585920803, + "Bz": 3371.758117878557, + "Ex": 55378581103426.695, + "Ey": 55378581103426.7, + "Ez": 68412803445328.25 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json index 1f9f0a77b5a..b03a954397a 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_3d_particle_scrape_picmi.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 144495.08082507108, - "By": 144495.08082507114, - "Bz": 8481.958724628861, - "Ex": 54500496182517.92, - "Ey": 54500496182517.91, - "Ez": 70231240245509.39 + "Bx": 148673.005859208, + "By": 148673.00585920803, + "Bz": 3371.758117878557, + "Ex": 55378581103426.695, + "Ey": 55378581103426.7, + "Ez": 68412803445328.25 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json b/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json index e4b9d9c07ff..0e5fad8db8a 100644 --- a/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json +++ b/Regression/Checksum/benchmarks_json/test_rz_embedded_boundary_diffraction.json @@ -1,10 +1,10 @@ { "lev=0": { - "Br": 6.7914286131989935e-19, - "Bt": 5.4557350206853276e-05, - "Bz": 2.357229221622199e-19, - "Er": 16481.39008058988, - "Et": 1.5258937379236053e-10, - "Ez": 1508.1064116028576 + "Br": 6.821267675779345e-19, + "Bt": 5.564905732478707e-05, + "Bz": 2.368259586613272e-19, + "Er": 16503.98082446463, + "Et": 1.5299584682447838e-10, + "Ez": 1466.854467399728 } } \ No newline at end of file From 3a3f36efe252fe60645dafe89146f56748f41892 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 17 Jan 2025 06:16:41 -0800 Subject: [PATCH 31/32] Use previous criterion to update EB --- Source/Initialization/WarpXInitData.cpp | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 248448e817e..147d2bfd0de 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1264,22 +1264,13 @@ void WarpX::InitializeEBGridData (int lev) // Compute additional quantities required for the ECT solver MarkExtensionCells(); ComputeFaceExtensions(); - // Mark on which grid points E should be updated - MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev ); - // Mark on which grid points B should be updated - MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); - } else { - // Mark on which grid points E should be updated (stair-case approximation) - MarkUpdateCellsStairCase( - m_eb_update_E[lev], - m_fields.get_alldirs(FieldType::Efield_fp, lev), - eb_fact ); - // Mark on which grid points B should be updated (stair-case approximation) - MarkUpdateCellsStairCase( - m_eb_update_B[lev], - m_fields.get_alldirs(FieldType::Bfield_fp, lev), - eb_fact ); } + + // Mark on which grid points E should be updated + MarkUpdateECellsECT( m_eb_update_E[lev], edge_lengths_lev ); + // Mark on which grid points B should be updated + MarkUpdateBCellsECT( m_eb_update_B[lev], face_areas_lev, edge_lengths_lev); + } ComputeDistanceToEB(); From 4b4e761ede75fe113a8b4bec6a07b712f4bc29e8 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 17 Jan 2025 06:18:32 -0800 Subject: [PATCH 32/32] Use previous input scripts --- Examples/Tests/embedded_boundary_cube/inputs_base_3d | 4 ++-- .../inputs_test_2d_embedded_boundary_cube | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Examples/Tests/embedded_boundary_cube/inputs_base_3d b/Examples/Tests/embedded_boundary_cube/inputs_base_3d index 70ddd8f8f64..90ae2996635 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_base_3d +++ b/Examples/Tests/embedded_boundary_cube/inputs_base_3d @@ -13,8 +13,8 @@ boundary.field_lo = pec pec pec boundary.field_hi = pec pec pec eb2.geom_type = box -eb2.box_lo = -0.501 -0.501 -0.501 # Ensures that the stair-case EB is exactly at -0.5 -eb2.box_hi = 0.501 0.501 0.501 # Ensures that the stair-case EB is exactly at 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. diff --git a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube index 46272052c2c..684325dc030 100644 --- a/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube +++ b/Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube @@ -12,10 +12,10 @@ warpx.abort_on_warning_threshold = medium boundary.field_lo = pec pec boundary.field_hi = pec pec -my_constants.xmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5 -my_constants.zmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5 -my_constants.xmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5 -my_constants.zmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5 +my_constants.xmin = -0.5 +my_constants.zmin = -0.5 +my_constants.xmax = 0.5 +my_constants.zmax = 0.5 # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular. warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))"