From 38d967d83b3f1dbc99114d1a329f6a8658334078 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Jan 2025 16:28:05 -0800 Subject: [PATCH] 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 ); }