Hybrid-PIC: Only calculate grad Pe when longitudinal part of E matters (

ECP-WarpX#5224) * rename `include_resistivity_term` to `solve_for_Faraday` * only include grad Pe term when the longitudinal part of E matters Signed-off-by: roelof-groenewald <[email protected]> * reset 2d benchmarks --------- Signed-off-by: roelof-groenewald <[email protected]>
JustinRayAngus · Sep 9, 2024 · 50e1326 · 50e1326
1 parent d254be6
commit 50e1326
Show file tree

Hide file tree

Showing 7 changed files with 84 additions and 88 deletions.
diff --git a/Regression/Checksum/benchmarks_json/test_2d_ohm_solver_landau_damping_picmi.json b/Regression/Checksum/benchmarks_json/test_2d_ohm_solver_landau_damping_picmi.json
@@ -1,21 +1,21 @@
 {
   "lev=0": {
     "Bx": 0.0,
-    "By": 7.079679609748623e-06,
+    "By": 0.0,
     "Bz": 0.0,
-    "Ex": 2726044.0536666242,
+    "Ex": 2726044.0531522455,
     "Ey": 0.0,
-    "Ez": 4060168.6414095857,
-    "jx": 177543428.8941278,
-    "jy": 187432087.03814715,
-    "jz": 594259755.4658135
+    "Ez": 4060168.641739453,
+    "jx": 177543428.93227622,
+    "jy": 187432087.0391676,
+    "jz": 594259755.4796929
   },
   "ions": {
-    "particle_momentum_x": 9.141594694084731e-17,
+    "particle_momentum_x": 9.14159469411586e-17,
     "particle_momentum_y": 9.135546407258978e-17,
-    "particle_momentum_z": 9.137866220861254e-17,
-    "particle_position_x": 1197.3344862524336,
-    "particle_position_y": 153269.17690371818,
+    "particle_momentum_z": 9.137866220831626e-17,
+    "particle_position_x": 1197.3344862525519,
+    "particle_position_y": 153269.1769037183,
     "particle_weight": 8.032598963696067e+16
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/test_2d_ohm_solver_magnetic_reconnection_picmi.json b/Regression/Checksum/benchmarks_json/test_2d_ohm_solver_magnetic_reconnection_picmi.json
@@ -1,18 +1,18 @@
 {
   "lev=0": {
-    "Bx": 1524.8789585503644,
+    "Bx": 1524.8789586064377,
     "By": 639.8314135126764,
-    "Bz": 7.476021044429126,
-    "Ex": 56499974.120022975,
-    "Ey": 75064070.52054195,
-    "Ez": 309548487.2375785
+    "Bz": 7.475443144751569,
+    "Ex": 56498215.028332025,
+    "Ey": 75045238.66546318,
+    "Ez": 309589044.22395706
   },
   "ions": {
-    "particle_momentum_x": 7.142955224072218e-15,
-    "particle_momentum_y": 7.138078059799475e-15,
-    "particle_momentum_z": 7.141134511766018e-15,
-    "particle_position_x": 11170689.202742986,
-    "particle_position_y": 5585328.083267269,
+    "particle_momentum_x": 7.142955285069654e-15,
+    "particle_momentum_y": 7.138077993721551e-15,
+    "particle_momentum_z": 7.141133937530391e-15,
+    "particle_position_x": 11170689.202379484,
+    "particle_position_y": 5585328.091626547,
     "particle_weight": 9.036667901693183e+18
   }
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -149,7 +149,7 @@ class FiniteDifferenceSolver
           * \param[in] edge_lengths length of edges along embedded boundaries
           * \param[in] lev  level number for the calculation
           * \param[in] hybrid_model instance of the hybrid-PIC model
-          * \param[in] include_resistivity_term boolean flag for whether to include resistivity
+          * \param[in] solve_for_Faraday boolean flag for whether the E-field is solved to be used in Faraday's equation
           */
         void HybridPICSolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3>& Jfield,
@@ -160,7 +160,7 @@ class FiniteDifferenceSolver
                       std::unique_ptr<amrex::MultiFab> const& Pefield,
                       std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
                       int lev, HybridPICModel const* hybrid_model,
-                      bool include_resistivity_term );
+                      bool solve_for_Faraday );
 
         /**
           * \brief Calculation of total current using Ampere's law (without
@@ -245,7 +245,7 @@ class FiniteDifferenceSolver
             std::unique_ptr<amrex::MultiFab> const& Pefield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
             int lev, HybridPICModel const* hybrid_model,
-            bool include_resistivity_term );
+            bool solve_for_Faraday );
 
         template<typename T_Algo>
         void CalculateCurrentAmpereCylindrical (
@@ -350,7 +350,7 @@ class FiniteDifferenceSolver
             std::unique_ptr<amrex::MultiFab> const& Pefield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
             int lev, HybridPICModel const* hybrid_model,
-            bool include_resistivity_term );
+            bool solve_for_Faraday );
 
         template<typename T_Algo>
         void CalculateCurrentAmpereCartesian (

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
@@ -91,23 +91,23 @@ 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,
-        bool include_resistivity_term);
+        bool solve_for_Faraday);
 
     void HybridPICSolveE (
         std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
         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,
-        int lev, bool include_resistivity_term);
+        int lev, bool solve_for_Faraday);
 
     void HybridPICSolveE (
         std::array< std::unique_ptr<amrex::MultiFab>, 3>& Efield,
         std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield,
         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,
-        int lev, PatchType patch_type, bool include_resistivity_term);
+        int lev, PatchType patch_type, bool solve_for_Faraday);
 
     void BfieldEvolveRK (
         amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
@@ -138,12 +138,11 @@ public:
 
     /**
      * \brief
-     * Function to calculate the electron pressure at a given timestep type
-     * using the simulation charge density. Used in the Ohm's law solver
-     * (kinetic-fluid hybrid model).
+     * Function to calculate the electron pressure using the simulation charge
+     * density. Used in the Ohm's law solver (kinetic-fluid hybrid model).
      */
-    void CalculateElectronPressure (          DtType a_dt_type);
-    void CalculateElectronPressure (int lev,  DtType a_dt_type);
+    void CalculateElectronPressure ();
+    void CalculateElectronPressure (int lev);
 
     /**
      * \brief Fill the electron pressure multifab given the kinetic particle

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
@@ -426,14 +426,14 @@ 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,
-    const bool include_resistivity_term)
+    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, include_resistivity_term
+            edge_lengths[lev], lev, solve_for_Faraday
         );
     }
 }
@@ -444,13 +444,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,
-    const int lev, const bool include_resistivity_term)
+    const int lev, const bool solve_for_Faraday)
 {
     WARPX_PROFILE("WarpX::HybridPICSolveE()");
 
     HybridPICSolveE(
         Efield, Jfield, Bfield, rhofield, edge_lengths, lev,
-        PatchType::fine, include_resistivity_term
+        PatchType::fine, solve_for_Faraday
     );
     if (lev > 0)
     {
@@ -466,7 +466,7 @@ void HybridPICModel::HybridPICSolveE (
     std::unique_ptr<amrex::MultiFab> const& rhofield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths,
     const int lev, PatchType patch_type,
-    const bool include_resistivity_term)
+    const bool solve_for_Faraday)
 {
     auto& warpx = WarpX::GetInstance();
 
@@ -475,36 +475,29 @@ void HybridPICModel::HybridPICSolveE (
         Efield, current_fp_ampere[lev], Jfield, current_fp_external[lev],
         Bfield, rhofield,
         electron_pressure_fp[lev],
-        edge_lengths, lev, this, include_resistivity_term
+        edge_lengths, lev, this, solve_for_Faraday
     );
     warpx.ApplyEfieldBoundary(lev, patch_type);
 }
 
-void HybridPICModel::CalculateElectronPressure(DtType a_dt_type)
+void HybridPICModel::CalculateElectronPressure()
 {
     auto& warpx = WarpX::GetInstance();
     for (int lev = 0; lev <= warpx.finestLevel(); ++lev)
     {
-        CalculateElectronPressure(lev, a_dt_type);
+        CalculateElectronPressure(lev);
     }
 }
 
-void HybridPICModel::CalculateElectronPressure(const int lev, DtType a_dt_type)
+void HybridPICModel::CalculateElectronPressure(const int lev)
 {
     WARPX_PROFILE("WarpX::CalculateElectronPressure()");
 
     auto& warpx = WarpX::GetInstance();
-    // The full step uses rho^{n+1}, otherwise use the old or averaged
-    // charge density.
-    if (a_dt_type == DtType::Full) {
-        FillElectronPressureMF(
-            electron_pressure_fp[lev], warpx.getFieldPointer(FieldType::rho_fp, lev)
-        );
-    } else {
-        FillElectronPressureMF(
-            electron_pressure_fp[lev], rho_fp_temp[lev].get()
-        );
-    }
+    // Calculate the electron pressure using rho^{n+1}.
+    FillElectronPressureMF(
+        electron_pressure_fp[lev], warpx.getFieldPointer(FieldType::rho_fp, lev)
+    );
     warpx.ApplyElectronPressureBoundary(lev, PatchType::fine);
     electron_pressure_fp[lev]->FillBoundary(warpx.Geom(lev).periodicity());
 }