Rebasing to latest development branch

Adding time value for particle field lookup in field substepping.
ECP-WarpX · Sep 18, 2024 · 5f75428 · 5f75428
1 parent 98fd8af
commit 5f75428
Show file tree

Hide file tree

Showing 5 changed files with 51 additions and 22 deletions.
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -159,6 +159,7 @@ class FiniteDifferenceSolver
                       std::unique_ptr<amrex::MultiFab> const& rhofield,
                       std::unique_ptr<amrex::MultiFab> const& Pefield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+                      amrex::Real t,
                       int lev, HybridPICModel const* hybrid_model,
                       bool solve_for_Faraday );
 
@@ -244,6 +245,7 @@ class FiniteDifferenceSolver
             std::unique_ptr<amrex::MultiFab> const& rhofield,
             std::unique_ptr<amrex::MultiFab> const& Pefield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+            amrex::Real t,
             int lev, HybridPICModel const* hybrid_model,
             bool solve_for_Faraday );
 
@@ -349,6 +351,7 @@ class FiniteDifferenceSolver
             std::unique_ptr<amrex::MultiFab> const& rhofield,
             std::unique_ptr<amrex::MultiFab> const& Pefield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+            amrex::Real t,
             int lev, HybridPICModel const* hybrid_model,
             bool solve_for_Faraday );
 

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
@@ -91,6 +91,7 @@ public:
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
         amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
+        amrex::Real t, 
         bool solve_for_Faraday);
 
     void HybridPICSolveE (
@@ -99,6 +100,7 @@ public:
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
         std::unique_ptr<amrex::MultiFab> const& rhofield,
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        amrex::Real t, 
         int lev, bool solve_for_Faraday);
 
     void HybridPICSolveE (
@@ -107,6 +109,7 @@ public:
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
         std::unique_ptr<amrex::MultiFab> const& rhofield,
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+        amrex::Real t, 
         int lev, PatchType patch_type, bool solve_for_Faraday);
 
     void BfieldEvolveRK (
@@ -115,7 +118,7 @@ public:
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
         amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-        amrex::Real dt, DtType a_dt_type,
+        amrex::Real t, amrex::Real dt, DtType a_dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
     void BfieldEvolveRK (
@@ -124,7 +127,7 @@ public:
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
         amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-        amrex::Real dt, int lev, DtType dt_type,
+        amrex::Real t, amrex::Real dt, int lev, DtType dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
     void FieldPush (
@@ -133,7 +136,7 @@ public:
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
         amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-        amrex::Real dt, DtType dt_type,
+        amrex::Real t, amrex::Real dt, DtType dt_type,
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
     /**

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
@@ -465,14 +465,15 @@ void HybridPICModel::HybridPICSolveE (
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield,
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
+    amrex::Real t,
     const bool solve_for_Faraday)
 {
     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
+            edge_lengths[lev], t, lev, solve_for_Faraday
         );
     }
 }
@@ -483,12 +484,13 @@ void HybridPICModel::HybridPICSolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+    amrex::Real t,
     const int lev, const bool solve_for_Faraday)
 {
     WARPX_PROFILE("WarpX::HybridPICSolveE()");
 
     HybridPICSolveE(
-        Efield, Jfield, Bfield, rhofield, edge_lengths, lev,
+        Efield, Jfield, Bfield, rhofield, edge_lengths, t, lev,
         PatchType::fine, solve_for_Faraday
     );
     if (lev > 0)
@@ -504,6 +506,7 @@ void HybridPICModel::HybridPICSolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield,
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
+    amrex::Real t,
     const int lev, PatchType patch_type,
     const bool solve_for_Faraday)
 {
@@ -514,7 +517,7 @@ void HybridPICModel::HybridPICSolveE (
         Efield, current_fp_ampere[lev], Jfield, current_fp_external[lev],
         Bfield, rhofield,
         electron_pressure_fp[lev],
-        edge_lengths, lev, this, solve_for_Faraday
+        edge_lengths, t, lev, this, solve_for_Faraday
     );
     warpx.ApplyEfieldBoundary(lev, patch_type);
 }
@@ -576,14 +579,14 @@ void HybridPICModel::BfieldEvolveRK (
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-    amrex::Real dt, DtType dt_type,
+    amrex::Real t, 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, edge_lengths, t, dt, lev, dt_type,
             ng, nodal_sync
         );
     }
@@ -595,7 +598,7 @@ void HybridPICModel::BfieldEvolveRK (
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-    amrex::Real dt, int lev, DtType dt_type,
+    amrex::Real t, amrex::Real dt, int lev, DtType dt_type,
     IntVect ng, std::optional<bool> nodal_sync )
 {
     // Make copies of the B-field multifabs at t = n and create multifabs for
@@ -618,11 +621,13 @@ void HybridPICModel::BfieldEvolveRK (
         K[ii].setVal(0.0);
     }
 
+    amrex::Real t_eval = t;
+
     // The Runge-Kutta scheme begins here.
     // Step 1:
     FieldPush(
         Bfield, Efield, Jfield, rhofield, edge_lengths,
-        0.5_rt*dt, dt_type, ng, nodal_sync
+        t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync
     );
 
     // The Bfield is now given by:
@@ -636,9 +641,10 @@ void HybridPICModel::BfieldEvolveRK (
     }
 
     // Step 2:
+    t_eval = t+0.5_rt*dt;
     FieldPush(
         Bfield, Efield, Jfield, rhofield, edge_lengths,
-        0.5_rt*dt, dt_type, ng, nodal_sync
+        t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync
     );
 
     // The Bfield is now given by:
@@ -658,7 +664,7 @@ void HybridPICModel::BfieldEvolveRK (
     // Step 3:
     FieldPush(
         Bfield, Efield, Jfield, rhofield, edge_lengths,
-        dt, dt_type, ng, nodal_sync
+        t_eval, dt, dt_type, ng, nodal_sync
     );
 
     // The Bfield is now given by:
@@ -672,9 +678,10 @@ void HybridPICModel::BfieldEvolveRK (
     }
 
     // Step 4:
+    t_eval = t + dt;
     FieldPush(
         Bfield, Efield, Jfield, rhofield, edge_lengths,
-        0.5_rt*dt, dt_type, ng, nodal_sync
+        t_eval, 0.5_rt*dt, dt_type, ng, nodal_sync
     );
 
     // The Bfield is now given by:
@@ -708,15 +715,15 @@ void HybridPICModel::FieldPush (
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield,
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield,
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths,
-    amrex::Real dt, DtType dt_type,
+    amrex::Real t, amrex::Real dt, DtType dt_type,
     IntVect ng, std::optional<bool> nodal_sync )
 {
     auto& warpx = WarpX::GetInstance();
 
     // Calculate J = curl x B / mu0
     CalculateCurrentAmpere(Bfield, edge_lengths);
     // Calculate the E-field from Ohm's law
-    HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
+    HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, t, true);
     warpx.FillBoundaryE(ng, nodal_sync);
     // Push forward the B-field using Faraday's law
     warpx.EvolveB(dt, dt_type);

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICSolveE.cpp
@@ -360,6 +360,7 @@ void FiniteDifferenceSolver::HybridPICSolveE (
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::unique_ptr<amrex::MultiFab> const& Pefield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    amrex::Real t,
     int lev, HybridPICModel const* hybrid_model,
     const bool solve_for_Faraday)
 {
@@ -370,14 +371,14 @@ void FiniteDifferenceSolver::HybridPICSolveE (
 
         HybridPICSolveECylindrical <CylindricalYeeAlgorithm> (
             Efield, Jfield, Jifield, Jextfield, Bfield, rhofield, Pefield,
-            edge_lengths, lev, hybrid_model, solve_for_Faraday
+            edge_lengths, t, lev, hybrid_model, solve_for_Faraday
         );
 
 #else
 
         HybridPICSolveECartesian <CartesianYeeAlgorithm> (
             Efield, Jfield, Jifield, Jextfield, Bfield, rhofield, Pefield,
-            edge_lengths, lev, hybrid_model, solve_for_Faraday
+            edge_lengths, t, lev, hybrid_model, solve_for_Faraday
         );
 
 #endif
@@ -398,6 +399,7 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::unique_ptr<amrex::MultiFab> const& Pefield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    amrex::Real t,
     int lev, HybridPICModel const* hybrid_model,
     const bool solve_for_Faraday )
 {
@@ -431,7 +433,6 @@ void FiniteDifferenceSolver::HybridPICSolveECylindrical (
     const auto Ez_part = hybrid_model->m_E_external[2];
 
     auto & warpx = WarpX::GetInstance();
-    auto t = warpx.gett_new(lev);
 
     auto dx_lev = warpx.Geom(lev).CellSizeArray();
     const RealBox& real_box = warpx.Geom(lev).ProbDomain();
@@ -776,6 +777,7 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::unique_ptr<amrex::MultiFab> const& Pefield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
+    amrex::Real t,
     int lev, HybridPICModel const* hybrid_model,
     const bool solve_for_Faraday )
 {
@@ -803,7 +805,6 @@ void FiniteDifferenceSolver::HybridPICSolveECartesian (
     const auto Ez_part = hybrid_model->m_E_external[2];
 
     auto & warpx = WarpX::GetInstance();
-    auto t = warpx.gett_new(lev);
 
     auto dx_lev = warpx.Geom(lev).CellSizeArray();
     const RealBox& real_box = warpx.Geom(lev).ProbDomain();

diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
@@ -88,14 +88,22 @@ void WarpX::HybridPICEvolveFields ()
         }
     }
 
+    // Calculate the electron pressure at t=n using rho^n
+    m_hybrid_pic_model->CalculateElectronPressure(DtType::FirstHalf);
+
+    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
     for (int sub_step = 0; sub_step < sub_steps; sub_step++)
     {
         m_hybrid_pic_model->BfieldEvolveRK(
             Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp,
-            m_edge_lengths, 0.5_rt/sub_steps*dt[0],
+            m_edge_lengths,
+            t_start + static_cast<amrex::Real>(sub_step)*sub_dt,
+            sub_dt,
             DtType::FirstHalf, guard_cells.ng_FieldSolver,
             WarpX::sync_nodal_points
         );
@@ -113,12 +121,19 @@ void WarpX::HybridPICEvolveFields ()
         );
     }
 
+    // Calculate the electron pressure at t=n+1/2
+    m_hybrid_pic_model->CalculateElectronPressure(DtType::SecondHalf);
+
+    t_start += 0.5_rt*dt[0];
+
     // 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++)
     {
         m_hybrid_pic_model->BfieldEvolveRK(
             Bfield_fp, Efield_fp, current_fp, rho_fp_temp,
-            m_edge_lengths, 0.5_rt/sub_steps*dt[0],
+            m_edge_lengths,
+            t_start + static_cast<amrex::Real>(sub_step)*sub_dt,
+            sub_dt,
             DtType::SecondHalf, guard_cells.ng_FieldSolver,
             WarpX::sync_nodal_points
         );
@@ -148,7 +163,7 @@ void WarpX::HybridPICEvolveFields ()
     // 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, false
+        Efield_fp, current_fp_temp, Bfield_fp, rho_fp, m_edge_lengths, gett_new(0), false
     );
     FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);