Use new stair-case approximation in hybrid solver

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<amrex::iMultiFab>,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<amrex::iMultiFab>,3> const& eb_update_E,
                       int lev );
 
     private:
@@ -243,15 +243,15 @@ 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<amrex::iMultiFab>,3> const& eb_update_E,
             int lev, HybridPICModel const* hybrid_model,
             bool solve_for_Faraday );
 
         template<typename T_Algo>
         void CalculateCurrentAmpereCylindrical (
             ablastr::fields::VectorField& Jfield,
             ablastr::fields::VectorField const& Bfield,
-            ablastr::fields::VectorField const& edge_lengths,
+            std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
             int lev
         );
 
@@ -347,15 +347,15 @@ 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<amrex::iMultiFab>,3> const& eb_update_E,
             int lev, HybridPICModel const* hybrid_model,
             bool solve_for_Faraday );
 
         template<typename T_Algo>
         void CalculateCurrentAmpereCartesian (
             ablastr::fields::VectorField& Jfield,
             ablastr::fields::VectorField const& Bfield,
-            ablastr::fields::VectorField const& edge_lengths,
+            std::array< std::unique_ptr<amrex::iMultiFab>,3> const& eb_update_E,
             int lev
         );
 #endif

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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E
     );
     void CalculatePlasmaCurrent (
         ablastr::fields::VectorField const& Bfield,
-        ablastr::fields::VectorField const& edge_lengths,
+        std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
         int lev
     );
 
@@ -83,31 +83,31 @@ public:
         ablastr::fields::MultiLevelVectorField const& Jfield,
         ablastr::fields::MultiLevelVectorField const& Bfield,
         ablastr::fields::MultiLevelScalarField const& rhofield,
-        ablastr::fields::MultiLevelVectorField const& edge_lengths,
+        amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
         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,
+        std::array< std::unique_ptr<amrex::iMultiFab>,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,
+        std::array< std::unique_ptr<amrex::iMultiFab>,3 >& eb_update_E,
         int lev, PatchType patch_type, bool solve_for_Faraday) const;
 
     void BfieldEvolveRK (
         ablastr::fields::MultiLevelVectorField const& Bfield,
         ablastr::fields::MultiLevelVectorField const& Efield,
         ablastr::fields::MultiLevelVectorField const& Jfield,
         ablastr::fields::MultiLevelScalarField const& rhofield,
-        ablastr::fields::MultiLevelVectorField const& edge_lengths,
+        amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
         amrex::Real dt, DtType a_dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
@@ -116,7 +116,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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
         amrex::Real dt, int lev, DtType dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
@@ -125,7 +125,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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
         amrex::Real dt, DtType dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp

@@ -250,26 +250,26 @@ void HybridPICModel::GetCurrentExternal ()
 
 void HybridPICModel::CalculatePlasmaCurrent (
     ablastr::fields::MultiLevelVectorField const& Bfield,
-    ablastr::fields::MultiLevelVectorField const& edge_lengths)
+    amrex::Vector<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& 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<amrex::iMultiFab>,3 >& eb_update_E,
     const int lev)
 {
     WARPX_PROFILE("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,15 +293,15 @@ 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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
     const bool solve_for_Faraday) const
 {
     auto& warpx = WarpX::GetInstance();
     for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
     {
         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<amrex::iMultiFab>,3 >& 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<amrex::iMultiFab>,3 >& 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,15 +411,15 @@ 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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
     amrex::Real dt, DtType dt_type,
     IntVect ng, std::optional<bool> nodal_sync )
 {
     auto& warpx = WarpX::GetInstance();
     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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
     amrex::Real dt, int lev, DtType dt_type,
     IntVect ng, std::optional<bool> 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<std::array< std::unique_ptr<amrex::iMultiFab>,3 > >& eb_update_E,
     amrex::Real dt, DtType dt_type,
     IntVect ng, std::optional<bool> 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);