From 1482f8bcf163aa2263d8c17354c77eefd8089b82 Mon Sep 17 00:00:00 2001 From: "S. Eric Clark" <25495882+clarkse@users.noreply.github.com> Date: Tue, 24 Sep 2024 12:48:32 -0700 Subject: [PATCH] Adding grid based split field advance for Hybrid solver to improve performance. --- Python/pywarpx/fields.py | 34 ++++ Python/pywarpx/picmi.py | 64 +++++++ .../HybridPICModel/HybridPICModel.H | 78 ++++++-- .../HybridPICModel/HybridPICModel.cpp | 177 ++++++++++++------ .../HybridPICSolveE.cpp | 152 ++++----------- .../FieldSolver/WarpXPushFieldsHybridPIC.cpp | 67 ++++++- Source/WarpX.cpp | 5 +- 7 files changed, 388 insertions(+), 189 deletions(-) diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py index cbdd8d4517a..518c13e9bb9 100644 --- a/Python/pywarpx/fields.py +++ b/Python/pywarpx/fields.py @@ -696,6 +696,40 @@ def BzFPExternalWrapper(level=0, include_ghosts=False): mf_name="Bfield_fp_external[z]", level=level, include_ghosts=include_ghosts ) +def ExHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Efield_hyb_external[x]", level=level, include_ghosts=include_ghosts + ) + + +def EyHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Efield_hyb_external[y]", level=level, include_ghosts=include_ghosts + ) + + +def EzHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Efield_hyb_external[z]", level=level, include_ghosts=include_ghosts + ) + + +def BxHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Bfield_hyb_external[x]", level=level, include_ghosts=include_ghosts + ) + + +def ByHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Bfield_hyb_external[y]", level=level, include_ghosts=include_ghosts + ) + + +def BzHybridExternalWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper( + mf_name="Bfield_hyb_external[z]", level=level, include_ghosts=include_ghosts + ) def JxFPWrapper(level=0, include_ghosts=False): return _MultiFABWrapper( diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 0d51a8723b4..c8c4c0433d3 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1783,6 +1783,9 @@ class HybridPICSolver(picmistandard.base._ClassWithInit): Jx/y/z_external_function: str Function of space and time specifying external (non-plasma) currents. + + Ex/y/z_external_function: str + Function of space and time specifying external (non-plasma) E-fields. """ def __init__( @@ -1798,6 +1801,12 @@ def __init__( Jx_external_function=None, Jy_external_function=None, Jz_external_function=None, + Ex_expression=None, + Ey_expression=None, + Ez_expression=None, + Bx_expression=None, + By_expression=None, + Bz_expression=None, **kw, ): self.grid = grid @@ -1816,6 +1825,24 @@ def __init__( self.Jy_external_function = Jy_external_function self.Jz_external_function = Jz_external_function + self.add_external_fields = None + + self.Ex_external_function = Ex_expression + self.Ey_external_function = Ey_expression + self.Ez_external_function = Ez_expression + + self.Bx_external_function = Bx_expression + self.By_external_function = By_expression + self.Bz_external_function = Bz_expression + + if (Ex_expression is not None + or Ey_expression is not None + or Ez_expression is not None + or Bx_expression is not None + or By_expression is not None + or Bz_expression is not None): + self.add_external_fields = True + # Handle keyword arguments used in expressions self.user_defined_kw = {} for k in list(kw.keys()): @@ -1864,6 +1891,43 @@ def solver_initialize_inputs(self): self.Jz_external_function, self.mangle_dict ), ) + pywarpx.hybridpicmodel.add_external_fields = self.add_external_fields + pywarpx.hybridpicmodel.__setattr__( + "Bx_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.Bx_external_function, self.mangle_dict + ), + ) + pywarpx.hybridpicmodel.__setattr__( + "By_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.By_external_function, self.mangle_dict + ), + ) + pywarpx.hybridpicmodel.__setattr__( + "Bz_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.Bz_external_function, self.mangle_dict + ), + ) + pywarpx.hybridpicmodel.__setattr__( + "Ex_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.Ex_external_function, self.mangle_dict + ), + ) + pywarpx.hybridpicmodel.__setattr__( + "Ey_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.Ey_external_function, self.mangle_dict + ), + ) + pywarpx.hybridpicmodel.__setattr__( + "Ez_external_grid_function(x,y,z,t)", + pywarpx.my_constants.mangle_expression( + self.Ez_external_function, self.mangle_dict + ), + ) class ElectrostaticSolver(picmistandard.PICMI_ElectrostaticSolver): diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H index 584125ac683..eaadf3e7a95 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H @@ -21,6 +21,9 @@ #include #include +#include +#include +#include #include @@ -38,10 +41,26 @@ public: /** Allocate hybrid-PIC specific multifabs. Called in constructor. */ void AllocateMFs (int nlevs_max); - void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, - int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho, - const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag, - const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag); + void AllocateLevelMFs ( + int lev, + const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, + const int ncomps, + const amrex::IntVect& ngJ, + const amrex::IntVect& ngRho, + const amrex::IntVect& ngE, + const amrex::IntVect& ngB, + const amrex::IntVect& jx_nodal_flag, + const amrex::IntVect& jy_nodal_flag, + const amrex::IntVect& jz_nodal_flag, + const amrex::IntVect& rho_nodal_flag, + const amrex::IntVect& Ex_nodal_flag, + const amrex::IntVect& Ey_nodal_flag, + const amrex::IntVect& Ez_nodal_flag, + const amrex::IntVect& Bx_nodal_flag, + const amrex::IntVect& By_nodal_flag, + const amrex::IntVect& Bz_nodal_flag + ); /** Helper function to clear values from hybrid-PIC specific multifabs. */ void ClearLevel (int lev); @@ -57,11 +76,39 @@ public: void GetCurrentExternal ( amrex::Vector, 3>> const& edge_lengths ); + void GetCurrentExternal ( std::array< std::unique_ptr, 3> const& edge_lengths, int lev ); + void GetFieldsExternal ( + amrex::Vector, 3>> const& edge_lengths + ); + + void GetFieldsExternal ( + amrex::Vector, 3>> const& edge_lengths, + amrex::Real t); + + void GetFieldsExternal ( + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev + ); + + void GetExternalFieldFromExpression ( + std::array< std::unique_ptr, 3> const& field, + std::array< amrex::ParserExecutor<4>, 3> const& expression, + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev + ); + + void GetExternalFieldFromExpression ( + std::array< std::unique_ptr, 3> const& field, + std::array< amrex::ParserExecutor<4>, 3> const& expression, + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev, amrex::Real t + ); + /** * \brief * Function to calculate the total current based on Ampere's law while @@ -152,7 +199,7 @@ public: * charge density (and assumption of quasi-neutrality) using the user * specified electron equation of state. * - * \param[out] Pe_filed scalar electron pressure MultiFab at a given level + * \param[out] Pe_field scalar electron pressure MultiFab at a given level * \param[in] rho_field scalar ion charge density Multifab at a given level */ void FillElectronPressureMF ( @@ -191,17 +238,17 @@ public: bool m_external_field_has_time_dependence = false; /** External E/B fields */ - bool m_add_ext_particle_B_field = false; - std::string m_Bx_ext_part_function = "0.0"; - std::string m_By_ext_part_function = "0.0"; - std::string m_Bz_ext_part_function = "0.0"; + bool m_add_external_fields = false; + + std::string m_Bx_ext_grid_function = "0.0"; + std::string m_By_ext_grid_function = "0.0"; + std::string m_Bz_ext_grid_function = "0.0"; std::array< std::unique_ptr, 3> m_B_external_parser; std::array< amrex::ParserExecutor<4>, 3> m_B_external; - bool m_add_ext_particle_E_field = false; - std::string m_Ex_ext_part_function = "0.0"; - std::string m_Ey_ext_part_function = "0.0"; - std::string m_Ez_ext_part_function = "0.0"; + std::string m_Ex_ext_grid_function = "0.0"; + std::string m_Ey_ext_grid_function = "0.0"; + std::string m_Ez_ext_grid_function = "0.0"; std::array< std::unique_ptr, 3> m_E_external_parser; std::array< amrex::ParserExecutor<4>, 3> m_E_external; @@ -212,6 +259,11 @@ public: amrex::Vector, 3 > > current_fp_external; amrex::Vector< std::unique_ptr > electron_pressure_fp; + amrex::Vector, 3 > > Bfield_hyb_external; + amrex::Vector, 3 > > Efield_hyb_external; + // amrex::Vector, 3 > > Bfield_hyb_self; + // amrex::Vector, 3 > > Efield_hyb_self; + // Helper functions to retrieve hybrid-PIC multifabs [[nodiscard]] amrex::MultiFab* get_pointer_current_fp_ampere (int lev, int direction) const diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp index 236dd708725..64e2ab0d445 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -57,13 +57,16 @@ void HybridPICModel::ReadParameters () pp_hybrid.query("Jz_external_grid_function(x,y,z,t)", m_Jz_ext_grid_function); // external fields - const ParmParse pp_part("particles"); - pp_hybrid.query("Bx_external_particle_function(x,y,z,t)", m_Bx_ext_part_function); - pp_hybrid.query("By_external_particle_function(x,y,z,t)", m_By_ext_part_function); - pp_hybrid.query("Bz_external_particle_function(x,y,z,t)", m_Bz_ext_part_function); - pp_hybrid.query("Ex_external_particle_function(x,y,z,t)", m_Ex_ext_part_function); - pp_hybrid.query("Ey_external_particle_function(x,y,z,t)", m_Ey_ext_part_function); - pp_hybrid.query("Ez_external_particle_function(x,y,z,t)", m_Ez_ext_part_function); + pp_hybrid.query("add_external_fields", m_add_external_fields); + + if (m_add_external_fields) { + pp_hybrid.query("Bx_external_grid_function(x,y,z,t)", m_Bx_ext_grid_function); + pp_hybrid.query("By_external_grid_function(x,y,z,t)", m_By_ext_grid_function); + pp_hybrid.query("Bz_external_grid_function(x,y,z,t)", m_Bz_ext_grid_function); + pp_hybrid.query("Ex_external_grid_function(x,y,z,t)", m_Ex_ext_grid_function); + pp_hybrid.query("Ey_external_grid_function(x,y,z,t)", m_Ey_ext_grid_function); + pp_hybrid.query("Ez_external_grid_function(x,y,z,t)", m_Ez_ext_grid_function); + } } void HybridPICModel::AllocateMFs (int nlevs_max) @@ -73,14 +76,30 @@ void HybridPICModel::AllocateMFs (int nlevs_max) current_fp_temp.resize(nlevs_max); current_fp_ampere.resize(nlevs_max); current_fp_external.resize(nlevs_max); + + if (m_add_external_fields) { + Bfield_hyb_external.resize(nlevs_max); + Efield_hyb_external.resize(nlevs_max); + // Bfield_hyb_self.resize(nlevs_max); + // Efield_hyb_self.resize(nlevs_max); + } } -void HybridPICModel::AllocateLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, - const int ncomps, const IntVect& ngJ, const IntVect& ngRho, - const IntVect& jx_nodal_flag, - const IntVect& jy_nodal_flag, - const IntVect& jz_nodal_flag, - const IntVect& rho_nodal_flag) +void HybridPICModel::AllocateLevelMFs ( + int lev, const BoxArray& ba, const DistributionMapping& dm, + const int ncomps, + const IntVect& ngJ, const IntVect& ngRho, + const IntVect& ngE, const IntVect& ngB, + const IntVect& jx_nodal_flag, + const IntVect& jy_nodal_flag, + const IntVect& jz_nodal_flag, + const IntVect& rho_nodal_flag, + const IntVect& Ex_nodal_flag, + const IntVect& Ey_nodal_flag, + const IntVect& Ez_nodal_flag, + const IntVect& Bx_nodal_flag, + const IntVect& By_nodal_flag, + const IntVect& Bz_nodal_flag) { // The "electron_pressure_fp" multifab stores the electron pressure calculated // from the specified equation of state. @@ -120,6 +139,22 @@ void HybridPICModel::AllocateLevelMFs (int lev, const BoxArray& ba, const Distri WarpX::AllocInitMultiFab(current_fp_external[lev][2], amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))), dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[z]", 0.0_rt); + if (m_add_external_fields) { + // These are nodal to match when B-field is added in evaluation of Ohm's law + WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][0], amrex::convert(ba, Bx_nodal_flag), + dm, ncomps, ngB, lev, "Bfield_hyb_external[x]", 0.0_rt); + WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][1], amrex::convert(ba, By_nodal_flag), + dm, ncomps, ngB, lev, "Bfield_hyb_external[y]", 0.0_rt); + WarpX::AllocInitMultiFab(Bfield_hyb_external[lev][2], amrex::convert(ba, Bz_nodal_flag), + dm, ncomps, ngB, lev, "Bfield_hyb_external[z]", 0.0_rt); + WarpX::AllocInitMultiFab(Efield_hyb_external[lev][0], amrex::convert(ba, Ex_nodal_flag), + dm, ncomps, ngE, lev, "Efield_hyb_external[x]", 0.0_rt); + WarpX::AllocInitMultiFab(Efield_hyb_external[lev][1], amrex::convert(ba, Ey_nodal_flag), + dm, ncomps, ngE, lev, "Efield_hyb_external[y]", 0.0_rt); + WarpX::AllocInitMultiFab(Efield_hyb_external[lev][2], amrex::convert(ba, Ez_nodal_flag), + dm, ncomps, ngE, lev, "Efield_hyb_external[z]", 0.0_rt); + } + #ifdef WARPX_DIM_RZ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (ncomps == 1), @@ -135,6 +170,10 @@ void HybridPICModel::ClearLevel (int lev) current_fp_temp[lev][i].reset(); current_fp_ampere[lev][i].reset(); current_fp_external[lev][i].reset(); + if (m_add_external_fields) { + Bfield_hyb_external[lev][i].reset(); + Efield_hyb_external[lev][i].reset(); + } } } @@ -165,33 +204,25 @@ void HybridPICModel::InitData () auto & warpx = WarpX::GetInstance(); const auto& mypc = warpx.GetPartContainer(); - if ( mypc.m_B_ext_particle_s == "parse_b_ext_particle_function") { - m_B_external_parser[0] = std::make_unique( - utils::parser::makeParser(m_Bx_ext_part_function,{"x","y","z","t"})); - m_B_external_parser[1] = std::make_unique( - utils::parser::makeParser(m_By_ext_part_function,{"x","y","z","t"})); - m_B_external_parser[2] = std::make_unique( - utils::parser::makeParser(m_Bz_ext_part_function,{"x","y","z","t"})); - m_B_external[0] = m_B_external_parser[0]->compile<4>(); - m_B_external[1] = m_B_external_parser[1]->compile<4>(); - m_B_external[2] = m_B_external_parser[2]->compile<4>(); - - m_add_ext_particle_B_field = true; - } - - if ( mypc.m_E_ext_particle_s == "parse_e_ext_particle_function") { - m_E_external_parser[0] = std::make_unique( - utils::parser::makeParser(m_Ex_ext_part_function,{"x","y","z","t"})); - m_E_external_parser[1] = std::make_unique( - utils::parser::makeParser(m_Ey_ext_part_function,{"x","y","z","t"})); - m_E_external_parser[2] = std::make_unique( - utils::parser::makeParser(m_Ez_ext_part_function,{"x","y","z","t"})); - m_E_external[0] = m_E_external_parser[0]->compile<4>(); - m_E_external[1] = m_E_external_parser[0]->compile<4>(); - m_E_external[2] = m_E_external_parser[0]->compile<4>(); - - m_add_ext_particle_E_field = true; - } + m_B_external_parser[0] = std::make_unique( + utils::parser::makeParser(m_Bx_ext_grid_function,{"x","y","z","t"})); + m_B_external_parser[1] = std::make_unique( + utils::parser::makeParser(m_By_ext_grid_function,{"x","y","z","t"})); + m_B_external_parser[2] = std::make_unique( + utils::parser::makeParser(m_Bz_ext_grid_function,{"x","y","z","t"})); + m_B_external[0] = m_B_external_parser[0]->compile<4>(); + m_B_external[1] = m_B_external_parser[1]->compile<4>(); + m_B_external[2] = m_B_external_parser[2]->compile<4>(); + + m_E_external_parser[0] = std::make_unique( + utils::parser::makeParser(m_Ex_ext_grid_function,{"x","y","z","t"})); + m_E_external_parser[1] = std::make_unique( + utils::parser::makeParser(m_Ey_ext_grid_function,{"x","y","z","t"})); + m_E_external_parser[2] = std::make_unique( + utils::parser::makeParser(m_Ez_ext_grid_function,{"x","y","z","t"})); + m_E_external[0] = m_E_external_parser[0]->compile<4>(); + m_E_external[1] = m_E_external_parser[1]->compile<4>(); + m_E_external[2] = m_E_external_parser[2]->compile<4>(); // Get the grid staggering of the fields involved in calculating E amrex::IntVect Jx_stag = warpx.getField(FieldType::current_fp, 0,0).ixType().toIntVect(); @@ -285,6 +316,7 @@ void HybridPICModel::InitData () } #endif GetCurrentExternal(edge_lengths, lev); + GetFieldsExternal(edge_lengths, lev); } } @@ -296,36 +328,73 @@ void HybridPICModel::GetCurrentExternal ( auto& warpx = WarpX::GetInstance(); for (int lev = 0; lev <= warpx.finestLevel(); ++lev) { - GetCurrentExternal(edge_lengths[lev], lev); + GetExternalFieldFromExpression(current_fp_external[lev], m_J_external, edge_lengths[lev], lev); } } - void HybridPICModel::GetCurrentExternal ( std::array< std::unique_ptr, 3> const& edge_lengths, int lev) { - // This logic matches closely to WarpX::InitializeExternalFieldsOnGridUsingParser - // except that the parsers include time dependence. - auto & warpx = WarpX::GetInstance(); + if (!m_external_field_has_time_dependence) { return; } + GetExternalFieldFromExpression(current_fp_external[lev], m_J_external, edge_lengths, lev); +} +void HybridPICModel::GetFieldsExternal ( + amrex::Vector, 3>> const& edge_lengths, + amrex::Real t) +{ + auto& warpx = WarpX::GetInstance(); + for (int lev = 0; lev <= warpx.finestLevel(); ++lev) + { + GetExternalFieldFromExpression(Bfield_hyb_external[lev], m_B_external, edge_lengths[lev], lev, t); + GetExternalFieldFromExpression(Efield_hyb_external[lev], m_E_external, edge_lengths[lev], lev, t); + } +} + +void HybridPICModel::GetFieldsExternal ( + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev) +{ + GetExternalFieldFromExpression(Bfield_hyb_external[lev], m_B_external, edge_lengths, lev); + GetExternalFieldFromExpression(Efield_hyb_external[lev], m_E_external, edge_lengths, lev); +} + +void HybridPICModel::GetExternalFieldFromExpression ( + std::array< std::unique_ptr, 3> const& field, + std::array< amrex::ParserExecutor<4>, 3> const& expression, + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev) +{ + auto & warpx = WarpX::GetInstance(); auto t = warpx.gett_new(lev); + GetExternalFieldFromExpression(field, expression, edge_lengths, lev, t); +} +void HybridPICModel::GetExternalFieldFromExpression ( + std::array< std::unique_ptr, 3> const& field, + std::array< amrex::ParserExecutor<4>, 3> const& expression, + std::array< std::unique_ptr, 3> const& edge_lengths, + int lev, amrex::Real t) +{ + // This logic matches closely to WarpX::InitializeExternalFieldsOnGridUsingParser + // except that the parsers include time dependence. + auto & warpx = WarpX::GetInstance(); auto dx_lev = warpx.Geom(lev).CellSizeArray(); const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - auto& mfx = current_fp_external[lev][0]; - auto& mfy = current_fp_external[lev][1]; - auto& mfz = current_fp_external[lev][2]; + auto& mfx = field[0]; + auto& mfy = field[1]; + auto& mfz = field[2]; const amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect(); const amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect(); const amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect(); // avoid implicit lambda capture - auto Jx_external = m_J_external[0]; - auto Jy_external = m_J_external[1]; - auto Jz_external = m_J_external[2]; + auto x_external = expression[0]; + auto y_external = expression[1]; + auto z_external = expression[2]; for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -371,7 +440,7 @@ void HybridPICModel::GetCurrentExternal ( const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; #endif // Initialize the x-component of the field. - mfxfab(i,j,k) = Jx_external(x,y,z,t); + mfxfab(i,j,k) = x_external(x,y,z,t); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // skip if node is covered by an embedded boundary @@ -397,7 +466,7 @@ void HybridPICModel::GetCurrentExternal ( const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; #endif // Initialize the y-component of the field. - mfyfab(i,j,k) = Jy_external(x,y,z,t); + mfyfab(i,j,k) = y_external(x,y,z,t); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // skip if node is covered by an embedded boundary @@ -423,7 +492,7 @@ void HybridPICModel::GetCurrentExternal ( const amrex::Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; #endif // Initialize the z-component of the field. - mfzfab(i,j,k) = Jz_external(x,y,z,t); + mfzfab(i,j,k) = z_external(x,y,z,t); } ); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp index ed47fee8433..924c3a17a88 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp @@ -422,21 +422,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( const bool include_hyper_resistivity_term = (eta_h > 0.0) && solve_for_Faraday; - const bool include_B_ext_part = hybrid_model->m_add_ext_particle_B_field; - const auto Br_part = hybrid_model->m_B_external[0]; - const auto Bt_part = hybrid_model->m_B_external[1]; - const auto Bz_part = hybrid_model->m_B_external[2]; - - const bool include_E_ext_part = hybrid_model->m_add_ext_particle_E_field; - const auto Er_part = hybrid_model->m_E_external[0]; - const auto Et_part = hybrid_model->m_E_external[1]; - const auto Ez_part = hybrid_model->m_E_external[2]; - - auto & warpx = WarpX::GetInstance(); - - auto dx_lev = warpx.Geom(lev).CellSizeArray(); - const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - const auto nodal_flag = IntVect::TheNodeVector(); + const bool include_external_fields = hybrid_model->m_add_external_fields; + auto const& Bfield_external = hybrid_model->Bfield_hyb_external[0]; // lev=0 + auto const& Efield_external = hybrid_model->Efield_hyb_external[0]; // lev=0 // Index type required for interpolating fields from their respective // staggering to the Ex, Ey, Ez locations @@ -497,6 +485,10 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( Array4 const& Bt = Bfield[1]->const_array(mfi); Array4 const& Bz = Bfield[2]->const_array(mfi); + Array4 const& Br_ext = Bfield_external[0]->const_array(mfi); + Array4 const& Bt_ext = Bfield_external[1]->const_array(mfi); + Array4 const& Bz_ext = Bfield_external[2]->const_array(mfi); + // Loop over the cells and update the nodal E field amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/){ @@ -515,17 +507,10 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( auto Bt_interp = Interp(Bt, Bt_stag, nodal, coarsen, i, j, 0, 0); auto Bz_interp = Interp(Bz, Bz_stag, nodal, coarsen, i, j, 0, 0); - if (include_B_ext_part) { - // Determine r and z on nodal mesh at i and j - const amrex::Real fac_x = (1._rt - nodal_flag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real yy = 0._rt; - const amrex::Real fac_z = (1._rt - nodal_flag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real zz = j*dx_lev[1] + real_box.lo(1) + fac_z; - - Br_interp += Br_part(xx,yy,zz,t); - Bt_interp += Bt_part(xx,yy,zz,t); - Bz_interp += Bz_part(xx,yy,zz,t); + if (include_external_fields) { + Br_interp += Interp(Br_ext, Br_stag, nodal, coarsen, i, j, 0, 0); + Bt_interp += Interp(Bt_ext, Bt_stag, nodal, coarsen, i, j, 0, 0); + Bz_interp += Interp(Bz_ext, Bz_stag, nodal, coarsen, i, j, 0, 0); } // calculate enE = (J - Ji) x B @@ -567,6 +552,9 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( Array4 const& Er = Efield[0]->array(mfi); Array4 const& Et = Efield[1]->array(mfi); Array4 const& Ez = Efield[2]->array(mfi); + Array4 const& Er_ext = Efield_external[0]->const_array(mfi); + Array4 const& Et_ext = Efield_external[1]->const_array(mfi); + Array4 const& Ez_ext = Efield_external[2]->const_array(mfi); Array4 const& Jr = Jfield[0]->const_array(mfi); Array4 const& Jt = Jfield[1]->const_array(mfi); Array4 const& Jz = Jfield[2]->const_array(mfi); @@ -639,15 +627,8 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( Er(i, j, 0) -= eta_h * nabla2Jr; } - if (include_E_ext_part) { - // Determine r and z on nodal mesh at i and j - const amrex::Real fac_x = (1._rt - Er_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real yy = 0._rt; - const amrex::Real fac_z = (1._rt - Er_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real zz = j*dx_lev[1] + real_box.lo(1) + fac_z; - - Er(i, j, 0) -= Er_part(xx,yy,zz,t); + if (include_external_fields) { + Er(i, j, 0) -= Er_ext(i, j, 0); } }, @@ -693,15 +674,8 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Note: Hyper-resisitivity should be revisited here when modal decomposition is implemented - if (include_E_ext_part) { - // Determine r and z on nodal mesh at i and j - const amrex::Real fac_x = (1._rt - Et_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real yy = 0._rt; - const amrex::Real fac_z = (1._rt - Et_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real zz = j*dx_lev[1] + real_box.lo(1) + fac_z; - - Et(i, j, 0) -= Et_part(xx,yy,zz,t); + if (include_external_fields) { + Et(i, j, 0) -= Et_ext(i, j, 0); } }, @@ -738,20 +712,13 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical ( // Add resistivity only if E field value is used to update B if (solve_for_Faraday) { Ez(i, j, 0) += eta(rho_val, jtot_val) * Jz(i, j, 0); } - if (include_hyper_resistivity_term) { + if (include_hyper_resistivity_term && solve_for_Faraday) { auto nabla2Jz = T_Algo::Dzz(Jz, coefs_z, n_coefs_z, i, j, 0, 0); Ez(i, j, 0) -= eta_h * nabla2Jz; } - if (include_E_ext_part) { - // Determine r and z on nodal mesh at i and j - const amrex::Real fac_x = (1._rt - Ez_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real yy = 0._rt; - const amrex::Real fac_z = (1._rt - Ez_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real zz = j*dx_lev[1] + real_box.lo(1) + fac_z; - - Ez(i, j, 0) -= Ez_part(xx,yy,zz,t); + if (include_external_fields) { + Ez(i, j, 0) -= Ez_ext(i, j, 0); } } ); @@ -794,21 +761,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( const bool include_hyper_resistivity_term = (eta_h > 0.) && solve_for_Faraday; - const bool include_B_ext_part = hybrid_model->m_add_ext_particle_B_field; - const auto Bx_part = hybrid_model->m_B_external[0]; - const auto By_part = hybrid_model->m_B_external[1]; - const auto Bz_part = hybrid_model->m_B_external[2]; - - const bool include_E_ext_part = hybrid_model->m_add_ext_particle_E_field; - const auto Ex_part = hybrid_model->m_E_external[0]; - const auto Ey_part = hybrid_model->m_E_external[1]; - const auto Ez_part = hybrid_model->m_E_external[2]; - - auto & warpx = WarpX::GetInstance(); - - auto dx_lev = warpx.Geom(lev).CellSizeArray(); - const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - const auto nodal_flag = IntVect::TheNodeVector(); + const bool include_external_fields = hybrid_model->m_add_external_fields; + auto const& Bfield_external = hybrid_model->Bfield_hyb_external[0]; // lev=0 + auto const& Efield_external = hybrid_model->Efield_hyb_external[0]; // lev=0 // Index type required for interpolating fields from their respective // staggering to the Ex, Ey, Ez locations @@ -869,6 +824,10 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Array4 const& By = Bfield[1]->const_array(mfi); Array4 const& Bz = Bfield[2]->const_array(mfi); + Array4 const& Bx_ext = Bfield_external[0]->const_array(mfi); + Array4 const& By_ext = Bfield_external[1]->const_array(mfi); + Array4 const& Bz_ext = Bfield_external[2]->const_array(mfi); + // Loop over the cells and update the nodal E field amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -887,18 +846,10 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( auto By_interp = Interp(By, By_stag, nodal, coarsen, i, j, k, 0); auto Bz_interp = Interp(Bz, Bz_stag, nodal, coarsen, i, j, k, 0); - if (include_B_ext_part) { - // Determine r and z on nodal mesh at i and j - const amrex::Real fac_x = (1._rt - nodal_flag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real fac_y = (1._rt - nodal_flag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real yy = j*dx_lev[1] + real_box.lo(1) + fac_y; - const amrex::Real fac_z = (1._rt - nodal_flag[2]) * dx_lev[2] * 0.5_rt; - const amrex::Real zz = k*dx_lev[2] + real_box.lo(2) + fac_z; - - Bx_interp += Bx_part(xx,yy,zz,t); - By_interp += By_part(xx,yy,zz,t); - Bz_interp += Bz_part(xx,yy,zz,t); + if (include_external_fields) { + Bx_interp += Interp(Bx_ext, Bx_stag, nodal, coarsen, i, j, k, 0); + By_interp += Interp(By_ext, By_stag, nodal, coarsen, i, j, k, 0); + Bz_interp += Interp(Bz_ext, Bz_stag, nodal, coarsen, i, j, k, 0); } // calculate enE = (J - Ji) x B @@ -940,6 +891,9 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Array4 const& Ex = Efield[0]->array(mfi); Array4 const& Ey = Efield[1]->array(mfi); Array4 const& Ez = Efield[2]->array(mfi); + Array4 const& Ex_ext = Efield_external[0]->const_array(mfi); + Array4 const& Ey_ext = Efield_external[1]->const_array(mfi); + Array4 const& Ez_ext = Efield_external[2]->const_array(mfi); Array4 const& Jx = Jfield[0]->const_array(mfi); Array4 const& Jy = Jfield[1]->const_array(mfi); Array4 const& Jz = Jfield[2]->const_array(mfi); @@ -1007,16 +961,8 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Ex(i, j, k) -= eta_h * nabla2Jx; } - if (include_E_ext_part) { - // Determine x, y, and z on nodal mesh at i, j, & k - const amrex::Real fac_x = (1._rt - Ex_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real fac_y = (1._rt - Ex_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real yy = j*dx_lev[1] + real_box.lo(1) + fac_y; - const amrex::Real fac_z = (1._rt - Ex_stag[2]) * dx_lev[2] * 0.5_rt; - const amrex::Real zz = k*dx_lev[2] + real_box.lo(2) + fac_z; - - Ex(i, j, k) -= Ex_part(xx,yy,zz,t); + if (include_external_fields) { + Ex(i, j, k) -= Ex_ext(i, j, k); } }, @@ -1063,16 +1009,8 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Ey(i, j, k) -= eta_h * nabla2Jy; } - if (include_E_ext_part) { - // Determine x, y, and z on nodal mesh at i, j, & k - const amrex::Real fac_x = (1._rt - Ey_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real fac_y = (1._rt - Ey_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real yy = j*dx_lev[1] + real_box.lo(1) + fac_y; - const amrex::Real fac_z = (1._rt - Ey_stag[2]) * dx_lev[2] * 0.5_rt; - const amrex::Real zz = k*dx_lev[2] + real_box.lo(2) + fac_z; - - Ey(i, j, k) -= Ey_part(xx,yy,zz,t); + if (include_external_fields) { + Ey(i, j, k) -= Ey_ext(i, j, k); } }, @@ -1115,16 +1053,8 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian ( Ez(i, j, k) -= eta_h * nabla2Jz; } - if (include_E_ext_part) { - // Determine x, y, and z on nodal mesh at i, j, & k - const amrex::Real fac_x = (1._rt - Ez_stag[0]) * dx_lev[0] * 0.5_rt; - const amrex::Real xx = i*dx_lev[0] + real_box.lo(0) + fac_x; - const amrex::Real fac_y = (1._rt - Ez_stag[1]) * dx_lev[1] * 0.5_rt; - const amrex::Real yy = j*dx_lev[1] + real_box.lo(1) + fac_y; - const amrex::Real fac_z = (1._rt - Ez_stag[2]) * dx_lev[2] * 0.5_rt; - const amrex::Real zz = k*dx_lev[2] + real_box.lo(2) + fac_z; - - Ez(i, j, k) -= Ez_part(xx,yy,zz,t); + if (include_external_fields) { + Ez(i, j, k) -= Ez_ext(i, j, k); } } ); diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index 0195c7a3356..97736b45de9 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -55,6 +55,31 @@ void WarpX::HybridPICEvolveFields () // Get requested number of substeps to use const int sub_steps = m_hybrid_pic_model->m_substeps; + + amrex::Real t_eval = gett_old(0); + amrex::Real sub_dt = 0.5_rt*dt[0]/sub_steps; + + const bool add_external_fields = m_hybrid_pic_model->m_add_external_fields; + auto const& Bfield_hyb_external = m_hybrid_pic_model->Bfield_hyb_external; + auto const& Efield_hyb_external = m_hybrid_pic_model->Efield_hyb_external; + + // Handle field splitting for Hybrid field push + if (add_external_fields) { + // Get the external fields + m_hybrid_pic_model->GetFieldsExternal(m_edge_lengths, t_eval); + + // If using split fields, subtract the external field at the old time + for (int lev = 0; lev <= finest_level; ++lev) { + for (int idim = 0; idim < 3; ++idim) { + MultiFab::Subtract( + *Bfield_fp[lev][idim], + *Bfield_hyb_external[lev][idim], + 0, 0, 1, + Bfield_hyb_external[lev][idim]->nGrowVect()); + } + } + FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); + } // Get the external current m_hybrid_pic_model->GetCurrentExternal(m_edge_lengths); @@ -88,9 +113,6 @@ void WarpX::HybridPICEvolveFields () } } - amrex::Real t_start = gett_old(0); - amrex::Real sub_dt = 0.5_rt/sub_steps*dt[0]; - // Push the B field from t=n to t=n+1/2 using the current and density // at t=n, while updating the E field along with B using the electron // momentum equation @@ -99,8 +121,7 @@ void WarpX::HybridPICEvolveFields () m_hybrid_pic_model->BfieldEvolveRK( Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp, m_edge_lengths, - t_start + static_cast(sub_step)*sub_dt, - sub_dt, + t_eval, sub_dt, DtType::FirstHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points ); @@ -118,7 +139,12 @@ void WarpX::HybridPICEvolveFields () ); } - t_start += 0.5_rt*dt[0]; + t_eval += 0.5_rt*dt[0]; + + if (add_external_fields) { + // Get the external fields + m_hybrid_pic_model->GetFieldsExternal(m_edge_lengths, t_eval); + } // Now push the B field from t=n+1/2 to t=n+1 using the n+1/2 quantities for (int sub_step = 0; sub_step < sub_steps; sub_step++) @@ -126,8 +152,7 @@ void WarpX::HybridPICEvolveFields () m_hybrid_pic_model->BfieldEvolveRK( Bfield_fp, Efield_fp, current_fp, rho_fp_temp, m_edge_lengths, - t_start + static_cast(sub_step)*sub_dt, - sub_dt, + t_eval, sub_dt, DtType::SecondHalf, guard_cells.ng_FieldSolver, WarpX::sync_nodal_points ); @@ -154,11 +179,35 @@ void WarpX::HybridPICEvolveFields () // Calculate the electron pressure at t=n+1 m_hybrid_pic_model->CalculateElectronPressure(); + t_eval = gett_new(0); + // Update the E field to t=n+1 using the extrapolated J_i^n+1 value m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths); m_hybrid_pic_model->HybridPICSolveE( - Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, gett_new(0), false + Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, t_eval, false ); + + // Handle field splitting for Hybrid field push + if (add_external_fields) { + m_hybrid_pic_model->GetFieldsExternal(m_edge_lengths, t_eval); + + // If using split fields, add the external field at the new time + for (int lev = 0; lev <= finest_level; ++lev) { + for (int idim = 0; idim < 3; ++idim) { + MultiFab::Add( + *Bfield_fp[lev][idim], + *Bfield_hyb_external[lev][idim], + 0, 0, 1, + Bfield_hyb_external[lev][idim]->nGrowVect()); + MultiFab::Add( + *Efield_fp[lev][idim], + *Efield_hyb_external[lev][idim], + 0, 0, 1, + Efield_hyb_external[lev][idim]->nGrowVect()); + } + } + FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); + } 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/WarpX.cpp b/Source/WarpX.cpp index 647b3e0a6b2..c78a5037212 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2357,8 +2357,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) { m_hybrid_pic_model->AllocateLevelMFs( - lev, ba, dm, ncomps, ngJ, ngRho, jx_nodal_flag, jy_nodal_flag, - jz_nodal_flag, rho_nodal_flag + lev, ba, dm, ncomps, ngJ, ngRho, ngEB, ngEB, jx_nodal_flag, jy_nodal_flag, + jz_nodal_flag, rho_nodal_flag, Ex_nodal_flag, Ey_nodal_flag, Ez_nodal_flag, + Bx_nodal_flag, By_nodal_flag, Bz_nodal_flag ); }