diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 7726a2ed5bd..bf1b4f102a6 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -103,7 +103,7 @@ class FiniteDifferenceSolver * \param[out] Efield vector of electric field MultiFabs updated at a given level * \param[in] Bfield vector of magnetic field MultiFabs at a given level * \param[in] Jfield vector of current density MultiFabs at a given level - * \param[in] edge_lengths length of edges along embedded boundaries + * \param[in] eb_update_E indicate in which cell E should be updated (related to embedded boundaries) * \param[in] dt timestep of the simulation * \param[in] macroscopic_properties contains user-defined properties of the medium. */ @@ -147,7 +147,7 @@ class FiniteDifferenceSolver * \param[in] Bfield vector of magnetic field MultiFabs at a given level * \param[in] rhofield scalar ion charge density Multifab at a given level * \param[in] Pefield scalar electron pressure MultiFab at a given level - * \param[in] edge_lengths length of edges along embedded boundaries + * \param[in] eb_update_E indicate in which cell E should be udpated (related to embedded boundaries) * \param[in] lev level number for the calculation * \param[in] hybrid_model instance of the hybrid-PIC model * \param[in] solve_for_Faraday boolean flag for whether the E-field is solved to be used in Faraday's equation @@ -158,7 +158,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); @@ -168,13 +168,13 @@ class FiniteDifferenceSolver * * \param[out] Jfield vector of current MultiFabs at a given level * \param[in] Bfield vector of magnetic field MultiFabs at a given level - * \param[in] edge_lengths length of edges along embedded boundaries + * \param[in] eb_update_E indicate in which cell E should be udpated (related to embedded boundaries) * \param[in] lev level number for the calculation */ void CalculateCurrentAmpere ( ablastr::fields::VectorField& Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev ); private: @@ -243,7 +243,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); @@ -251,7 +251,7 @@ class FiniteDifferenceSolver void CalculateCurrentAmpereCylindrical ( ablastr::fields::VectorField& Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev ); @@ -329,7 +329,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField Efield, std::array< amrex::MultiFab*, 3 > Bfield, std::array< amrex::MultiFab*, 3 > Jfield, - std::array< amrex::MultiFab*, 3 > edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, amrex::MultiFab* Ffield, MultiSigmaBox const& sigba, amrex::Real dt, bool pml_has_particles ); @@ -347,7 +347,7 @@ class FiniteDifferenceSolver ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev, HybridPICModel const* hybrid_model, bool solve_for_Faraday ); @@ -355,7 +355,7 @@ class FiniteDifferenceSolver void CalculateCurrentAmpereCartesian ( ablastr::fields::VectorField& Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3> const& eb_update_E, int lev ); #endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H index 7e8dd260a6e..7ecf33b4945 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H @@ -62,15 +62,15 @@ public: * subtracted as well. Used in the Ohm's law solver (kinetic-fluid hybrid model). * * \param[in] Bfield Magnetic field from which the current is calculated. - * \param[in] edge_lengths Length of cell edges taking embedded boundaries into account + * \param[in] eb_update_E Indicate in which cell J should be calculated (related to embedded boundaries). */ void CalculatePlasmaCurrent ( ablastr::fields::MultiLevelVectorField const& Bfield, - ablastr::fields::MultiLevelVectorField const& edge_lengths + amrex::Vector,3 > > eb_update_E ); void CalculatePlasmaCurrent ( ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, int lev ); @@ -83,7 +83,7 @@ public: ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelVectorField const& Bfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, bool solve_for_Faraday) const; void HybridPICSolveE ( @@ -91,15 +91,14 @@ public: ablastr::fields::VectorField const& Jfield, ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, - ablastr::fields::VectorField const& edge_lengths, - int lev, bool solve_for_Faraday) const; + amrex::Vector,3 > > eb_update_E, int lev, bool solve_for_Faraday) const; void HybridPICSolveE ( ablastr::fields::VectorField const& Efield, ablastr::fields::VectorField const& Jfield, ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, - ablastr::fields::VectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, int lev, PatchType patch_type, bool solve_for_Faraday) const; void BfieldEvolveRK ( @@ -107,7 +106,7 @@ public: ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, amrex::Real dt, DtType a_dt_type, amrex::IntVect ng, std::optional nodal_sync); @@ -116,7 +115,7 @@ public: ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, amrex::Real dt, int lev, DtType dt_type, amrex::IntVect ng, std::optional nodal_sync); @@ -125,7 +124,7 @@ public: ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > eb_update_E, amrex::Real dt, DtType dt_type, amrex::IntVect ng, std::optional nodal_sync); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index abda59e40ba..56d3a2e7e2b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -250,18 +250,18 @@ void HybridPICModel::GetCurrentExternal () void HybridPICModel::CalculatePlasmaCurrent ( ablastr::fields::MultiLevelVectorField const& Bfield, - ablastr::fields::MultiLevelVectorField const& edge_lengths) + amrex::Vector,3 > > const& eb_update_E) { auto& warpx = WarpX::GetInstance(); for (int lev = 0; lev <= warpx.finestLevel(); ++lev) { - CalculatePlasmaCurrent(Bfield[lev], edge_lengths[lev], lev); + CalculatePlasmaCurrent(Bfield[lev], eb_update_E[lev], lev); } } void HybridPICModel::CalculatePlasmaCurrent ( ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > > const& eb_update_E, const int lev) { WARPX_PROFILE("HybridPICModel::CalculatePlasmaCurrent()"); @@ -269,7 +269,7 @@ void HybridPICModel::CalculatePlasmaCurrent ( auto& warpx = WarpX::GetInstance(); ablastr::fields::VectorField current_fp_plasma = warpx.m_fields.get_alldirs(FieldType::hybrid_current_fp_plasma, lev); warpx.get_pointer_fdtd_solver_fp(lev)->CalculateCurrentAmpere( - current_fp_plasma, Bfield, edge_lengths, lev + current_fp_plasma, Bfield, eb_update_E, lev ); // we shouldn't apply the boundary condition to J since J = J_i - J_e but @@ -293,7 +293,7 @@ void HybridPICModel::HybridPICSolveE ( ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelVectorField const& Bfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > const& eb_update_E, const bool solve_for_Faraday) const { auto& warpx = WarpX::GetInstance(); @@ -301,7 +301,7 @@ void HybridPICModel::HybridPICSolveE ( { HybridPICSolveE( Efield[lev], Jfield[lev], Bfield[lev], *rhofield[lev], - edge_lengths[lev], lev, solve_for_Faraday + eb_update_E[lev], lev, solve_for_Faraday ); } } @@ -311,13 +311,13 @@ void HybridPICModel::HybridPICSolveE ( ablastr::fields::VectorField const& Jfield, ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, const int lev, const bool solve_for_Faraday) const { WARPX_PROFILE("WarpX::HybridPICSolveE()"); HybridPICSolveE( - Efield, Jfield, Bfield, rhofield, edge_lengths, lev, + Efield, Jfield, Bfield, rhofield, eb_update_E, lev, PatchType::fine, solve_for_Faraday ); if (lev > 0) @@ -332,7 +332,7 @@ void HybridPICModel::HybridPICSolveE ( ablastr::fields::VectorField const& Jfield, ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, const int lev, PatchType patch_type, const bool solve_for_Faraday) const { @@ -344,7 +344,7 @@ void HybridPICModel::HybridPICSolveE ( // Solve E field in regular cells warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE( Efield, current_fp_plasma, Jfield, Bfield, rhofield, - *electron_pressure_fp, edge_lengths, lev, this, solve_for_Faraday + *electron_pressure_fp, eb_update_E, lev, this, solve_for_Faraday ); amrex::Real const time = warpx.gett_old(0) + warpx.getdt(0); warpx.ApplyEfieldBoundary(lev, patch_type, time); @@ -411,7 +411,7 @@ void HybridPICModel::BfieldEvolveRK ( ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > const& eb_update_E, amrex::Real dt, DtType dt_type, IntVect ng, std::optional nodal_sync ) { @@ -419,7 +419,7 @@ void HybridPICModel::BfieldEvolveRK ( for (int lev = 0; lev <= warpx.finestLevel(); ++lev) { BfieldEvolveRK( - Bfield, Efield, Jfield, rhofield, edge_lengths, dt, lev, dt_type, + Bfield, Efield, Jfield, rhofield, eb_update_E, dt, lev, dt_type, ng, nodal_sync ); } @@ -430,7 +430,7 @@ void HybridPICModel::BfieldEvolveRK ( ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > const& eb_update_E, amrex::Real dt, int lev, DtType dt_type, IntVect ng, std::optional nodal_sync ) { @@ -457,7 +457,7 @@ void HybridPICModel::BfieldEvolveRK ( // The Runge-Kutta scheme begins here. // Step 1: FieldPush( - Bfield, Efield, Jfield, rhofield, edge_lengths, + Bfield, Efield, Jfield, rhofield, eb_update_E, 0.5_rt*dt, dt_type, ng, nodal_sync ); @@ -473,7 +473,7 @@ void HybridPICModel::BfieldEvolveRK ( // Step 2: FieldPush( - Bfield, Efield, Jfield, rhofield, edge_lengths, + Bfield, Efield, Jfield, rhofield, eb_update_E, 0.5_rt*dt, dt_type, ng, nodal_sync ); @@ -493,7 +493,7 @@ void HybridPICModel::BfieldEvolveRK ( // Step 3: FieldPush( - Bfield, Efield, Jfield, rhofield, edge_lengths, + Bfield, Efield, Jfield, rhofield, eb_update_E, dt, dt_type, ng, nodal_sync ); @@ -509,7 +509,7 @@ void HybridPICModel::BfieldEvolveRK ( // Step 4: FieldPush( - Bfield, Efield, Jfield, rhofield, edge_lengths, + Bfield, Efield, Jfield, rhofield, eb_update_E, 0.5_rt*dt, dt_type, ng, nodal_sync ); @@ -543,16 +543,16 @@ void HybridPICModel::FieldPush ( ablastr::fields::MultiLevelVectorField const& Efield, ablastr::fields::MultiLevelVectorField const& Jfield, ablastr::fields::MultiLevelScalarField const& rhofield, - ablastr::fields::MultiLevelVectorField const& edge_lengths, + amrex::Vector,3 > > const& eb_update_E, amrex::Real dt, DtType dt_type, IntVect ng, std::optional nodal_sync ) { auto& warpx = WarpX::GetInstance(); // Calculate J = curl x B / mu0 - J_ext - CalculatePlasmaCurrent(Bfield, edge_lengths); + CalculatePlasmaCurrent(Bfield, eb_update_E); // Calculate the E-field from Ohm's law - HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true); + HybridPICSolveE(Efield, Jfield, Bfield, rhofield, eb_update_E, true); warpx.FillBoundaryE(ng, nodal_sync); // Push forward the B-field using Faraday's law amrex::Real const t_old = warpx.gett_old(0); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index 76fedbf4dea..3d3e2c79256 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -26,7 +26,7 @@ using namespace amrex; void FiniteDifferenceSolver::CalculateCurrentAmpere ( ablastr::fields::VectorField & Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev ) { // Select algorithm (The choice of algorithm is a runtime option, @@ -34,12 +34,12 @@ void FiniteDifferenceSolver::CalculateCurrentAmpere ( if (m_fdtd_algo == ElectromagneticSolverAlgo::HybridPIC) { #ifdef WARPX_DIM_RZ CalculateCurrentAmpereCylindrical ( - Jfield, Bfield, edge_lengths, lev + Jfield, Bfield, eb_update_E, lev ); #else CalculateCurrentAmpereCartesian ( - Jfield, Bfield, edge_lengths, lev + Jfield, Bfield, eb_update_E, lev ); #endif @@ -61,7 +61,7 @@ template void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical ( ablastr::fields::VectorField& Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev ) { @@ -92,12 +92,13 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical ( Array4 const& Bt = Bfield[1]->array(mfi); Array4 const& Bz = Bfield[2]->array(mfi); - amrex::Array4 lr, lt, lz; - + // Extract structures indicating where the fields + // should be updated, given the position of the embedded boundaries + amrex::Array4 update_Jr_arr, update_Jt_arr, update_Jz_arr; if (EB::enabled()) { - lr = edge_lengths[0]->array(mfi); - lt = edge_lengths[1]->array(mfi); - lz = edge_lengths[2]->array(mfi); + update_Jr_arr = eb_update_E[0]->array(mfi); + update_Jt_arr = eb_update_E[1]->array(mfi); + update_Jz_arr = eb_update_E[2]->array(mfi); } // Extract stencil coefficients @@ -124,8 +125,10 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical ( // Jr calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip if this cell is fully covered by embedded boundaries - if (lr && lr(i, j, 0) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Jr_arr && update_Jr_arr(i, j, 0) == 0) { return; } + // Mode m=0 Jr(i, j, 0, 0) = one_over_mu0 * ( - T_Algo::DownwardDz(Bt, coefs_z, n_coefs_z, i, j, 0, 0) @@ -148,8 +151,9 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical ( // Jt calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // In RZ Jt is associated with a mesh node, so we need to check if the mesh 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 update in the embedded boundaries + if (update_Jt_arr && update_Jt_arr(i, j, 0) == 0) { return; } // r on a nodal point (Jt is nodal in r) Real const r = rmin + i*dr; @@ -194,8 +198,10 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCylindrical ( // Jz calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip if this cell is fully covered by embedded boundaries - if (lz && lz(i, j, 0) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Jz_arr && update_Jz_arr(i, j, 0) == 0) { return; } + // r on a nodal point (Jz is nodal in r) Real const r = rmin + i*dr; // Off-axis, regular curl @@ -244,7 +250,7 @@ template void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian ( ablastr::fields::VectorField& Jfield, ablastr::fields::VectorField const& Bfield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev ) { @@ -274,11 +280,13 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian ( Array4 const &By = Bfield[1]->const_array(mfi); Array4 const &Bz = Bfield[2]->const_array(mfi); - amrex::Array4 lx, ly, lz; + // Extract structures indicating where the fields + // should be updated, given the position of the embedded boundaries + amrex::Array4 update_Jx_arr, update_Jy_arr, update_Jz_arr; if (EB::enabled()) { - lx = edge_lengths[0]->array(mfi); - ly = edge_lengths[1]->array(mfi); - lz = edge_lengths[2]->array(mfi); + update_Jx_arr = eb_update_E[0]->array(mfi); + update_Jy_arr = eb_update_E[1]->array(mfi); + update_Jz_arr = eb_update_E[2]->array(mfi); } // Extract stencil coefficients @@ -302,8 +310,9 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian ( // Jx calculation [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip if this cell is fully covered by embedded boundaries - if (lx && lx(i, j, k) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Jx_arr && update_Jx_arr(i, j, k) == 0) { return; } Jx(i, j, k) = one_over_mu0 * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k) @@ -313,14 +322,10 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian ( // Jy calculation [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip if this cell is fully covered by embedded boundaries -#ifdef WARPX_DIM_3D - if (ly && ly(i,j,k) <= 0) { return; } -#elif defined(WARPX_DIM_XZ) - // In XZ Jy 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; } -#endif + + // Skip field update in the embedded boundaries + if (update_Jy_arr && update_Jy_arr(i, j, k) == 0) { return; } + Jy(i, j, k) = one_over_mu0 * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k) @@ -329,8 +334,9 @@ void FiniteDifferenceSolver::CalculateCurrentAmpereCartesian ( // Jz calculation [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip if this cell is fully covered by embedded boundaries - if (lz && lz(i,j,k) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Jz_arr && update_Jz_arr(i, j, k) == 0) { return; } Jz(i, j, k) = one_over_mu0 * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k) @@ -357,7 +363,7 @@ void FiniteDifferenceSolver::HybridPICSolveE ( ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday) { @@ -368,14 +374,14 @@ void FiniteDifferenceSolver::HybridPICSolveE ( HybridPICSolveECylindrical ( Efield, Jfield, Jifield, Bfield, rhofield, Pefield, - edge_lengths, lev, hybrid_model, solve_for_Faraday + eb_update_E, lev, hybrid_model, solve_for_Faraday ); #else HybridPICSolveECartesian ( Efield, Jfield, Jifield, Bfield, rhofield, Pefield, - edge_lengths, lev, hybrid_model, solve_for_Faraday + eb_update_E, lev, hybrid_model, solve_for_Faraday ); #endif @@ -394,7 +400,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday ) { @@ -537,11 +543,13 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( Array4 const& rho = rhofield.const_array(mfi); Array4 const& Pe = Pefield.const_array(mfi); - amrex::Array4 lr, lz; + // Extract structures indicating where the fields + // should be updated, given the position of the embedded boundaries + amrex::Array4 update_Er_arr, update_Et_arr, update_Ez_arr; if (EB::enabled()) { - lr = edge_lengths[0]->array(mfi); - // edge_lengths[1] is `lt` and is not needed - 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 @@ -563,8 +571,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Er calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip if this cell is fully covered by embedded boundaries - if (lr && lr(i, j, 0) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Er_arr && update_Er_arr(i, j, 0) == 0) { return; } // Interpolate to get the appropriate charge density in space Real rho_val = Interp(rho, nodal, Er_stag, coarsen, i, j, 0, 0); @@ -605,8 +614,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Et calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // In RZ Et is associated with a mesh node, so we need to check if the mesh 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 update in the embedded boundaries + if (update_Et_arr && update_Et_arr(i, j, 0) == 0) { return; } // r on a nodal grid (Et is nodal in r) Real const r = rmin + i*dr; @@ -648,8 +658,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Ez calculation [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ - // Skip field solve if this cell is fully covered by embedded boundaries - if (lz && lz(i,j,0) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Ez_arr && update_Ez_arr(i, j, 0) == 0) { return; } // Interpolate to get the appropriate charge density in space Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, 0, 0); @@ -705,7 +716,7 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( ablastr::fields::VectorField const& Bfield, amrex::MultiFab const& rhofield, amrex::MultiFab const& Pefield, - ablastr::fields::VectorField const& edge_lengths, + std::array< std::unique_ptr,3 > const& eb_update_E, int lev, HybridPICModel const* hybrid_model, const bool solve_for_Faraday ) { @@ -842,11 +853,13 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Array4 const& rho = rhofield.const_array(mfi); Array4 const& Pe = Pefield.array(mfi); - amrex::Array4 lx, ly, lz; + // Extract structures indicating where the fields + // should be updated, given the position of the embedded boundaries + amrex::Array4 update_Ex_arr, update_Ey_arr, update_Ez_arr; if (EB::enabled()) { - lx = edge_lengths[0]->array(mfi); - ly = edge_lengths[1]->array(mfi); - 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); } // Extract stencil coefficients @@ -866,8 +879,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Ex calculation [=] AMREX_GPU_DEVICE (int i, int j, int k){ - // Skip if this cell is fully covered by embedded boundaries - if (lx && lx(i, j, k) <= 0) { return; } + + // Skip field update in the embedded boundaries + if (update_Ex_arr && update_Ex_arr(i, j, k) == 0) { return; } // Interpolate to get the appropriate charge density in space Real rho_val = Interp(rho, nodal, Ex_stag, coarsen, i, j, k, 0); @@ -905,14 +919,10 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Ey calculation [=] AMREX_GPU_DEVICE (int i, int j, int k) { - // Skip field solve if this cell is fully covered by embedded boundaries -#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; } -#endif + + // Skip field update in the embedded boundaries + if (update_Ey_arr && update_Ey_arr(i, j, k) == 0) { return; } + // Interpolate to get the appropriate charge density in space Real rho_val = Interp(rho, nodal, Ey_stag, coarsen, i, j, k, 0); @@ -949,10 +959,10 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( // Ez calculation [=] AMREX_GPU_DEVICE (int i, int j, int k){ -#ifdef AMREX_USE_EB - // Skip field solve if this cell is fully covered by embedded boundaries - if (lz && lz(i,j,k) <= 0) { return; } -#endif + + // Skip field update in the embedded boundaries + if (update_Ez_arr && update_Ez_arr(i, j, k) == 0) { return; } + // Interpolate to get the appropriate charge density in space Real rho_val = Interp(rho, nodal, Ez_stag, coarsen, i, j, k, 0); diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index 8e9e0daa274..46950030322 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -108,7 +108,7 @@ void WarpX::HybridPICEvolveFields () m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level), m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), current_fp_temp, rho_fp_temp, - m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level), + m_eb_update_E, 0.5_rt/sub_steps*dt[0], DtType::FirstHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points @@ -135,7 +135,7 @@ void WarpX::HybridPICEvolveFields () m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), rho_fp_temp, - m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level), + m_eb_update_E, 0.5_rt/sub_steps*dt[0], DtType::SecondHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points @@ -166,14 +166,13 @@ void WarpX::HybridPICEvolveFields () // Update the E field to t=n+1 using the extrapolated J_i^n+1 value m_hybrid_pic_model->CalculatePlasmaCurrent( m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level), - m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level)); + m_eb_update_E); m_hybrid_pic_model->HybridPICSolveE( m_fields.get_mr_levels_alldirs(FieldType::Efield_fp, finest_level), current_fp_temp, m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level), m_fields.get_mr_levels(FieldType::rho_fp, finest_level), - m_fields.get_mr_levels_alldirs(FieldType::edge_lengths, finest_level), false - ); + m_eb_update_E, false); FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); // Copy the rho^{n+1} values to rho_fp_temp and the J_i^{n+1/2} values to diff --git a/Source/Fields.H b/Source/Fields.H index b07661254c4..77589c4675e 100644 --- a/Source/Fields.H +++ b/Source/Fields.H @@ -61,8 +61,8 @@ namespace warpx::fields E_external_particle_field, /**< Stores external particle fields provided by the user as through an openPMD file */ B_external_particle_field, /**< Stores external particle fields provided by the user as through an openPMD file */ distance_to_eb, /**< Only used with embedded boundaries (EB). Stores the distance to the nearest EB */ - edge_lengths, /**< Only used with embedded boundaries (EB). Indicates the length of the cell edge that is covered by the EB, in SI units */ - face_areas, /**< Only used with embedded boundaries (EB). Indicates the area of the cell face that is covered by the EB, in SI units */ + edge_lengths, /**< Only used with the ECT solver. Indicates the length of the cell edge that is covered by the EB, in SI units */ + face_areas, /**< Only used with the ECT solver. Indicates the area of the cell face that is covered by the EB, in SI units */ area_mod, pml_E_fp, pml_B_fp, diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index c8c7d1f0337..fbbe18817ac 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1253,24 +1253,27 @@ 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)); + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { - auto face_areas_lev = m_fields.get_alldirs(FieldType::face_areas, lev); - ComputeFaceAreas(face_areas_lev, eb_fact); - ScaleAreas(face_areas_lev, CellSize(lev)); + auto edge_lengths_lev = m_fields.get_alldirs(FieldType::edge_lengths, lev); + ComputeEdgeLengths(edge_lengths_lev, eb_fact); + ScaleEdges(edge_lengths_lev, CellSize(lev)); + + auto face_areas_lev = m_fields.get_alldirs(FieldType::face_areas, lev); + ComputeFaceAreas(face_areas_lev, eb_fact); + ScaleAreas(face_areas_lev, CellSize(lev)); - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { // 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( m_eb_update_E[lev], @@ -1281,6 +1284,7 @@ void WarpX::InitializeEBGridData (int lev) 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 ad28c672ad5..4d0da668891 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -927,8 +927,7 @@ public: * \param[in] fx_parser parser function to initialize x-field * \param[in] fy_parser parser function to initialize y-field * \param[in] fz_parser parser function to initialize z-field - * \param[in] edge_lengths edge lengths information - * \param[in] face_areas face areas information + * \param[in] eb_update_field flag indicating which gridpoints should be modified by this functions * \param[in] topology flag indicating if field is edge-based or face-based * \param[in] lev level of the Multifabs that is initialized * \param[in] patch_type PatchType on which the field is initialized (fine or coarse) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 265603eae24..2f56f524b63 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2317,8 +2317,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm 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]"); + } + if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { - // 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); @@ -2334,8 +2335,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); m_fields.alloc_init(FieldType::face_areas, Direction{2}, lev, amrex::convert(ba, Bz_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, 0.0_rt); - } - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) { + AllocInitMultiFab(m_flag_info_face[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, guard_cells.ng_FieldSolver, lev, "m_flag_info_face[x]"); AllocInitMultiFab(m_flag_info_face[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps,