Various code cleanups

roelof-groenewald · Dec 2, 2023 · 4e94c15 · 4e94c15
1 parent 5716f3f
commit 4e94c15
Show file tree

Hide file tree

Showing 6 changed files with 91 additions and 104 deletions.
diff --git a/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py b/Examples/Tests/ohm_solver_ion_Landau_damping/PICMI_inputs.py
@@ -83,7 +83,7 @@ def __init__(self, test, dim, m, T_ratio, verbose):
         self.total_steps = int(np.ceil(self.LT / self.DT))
         # if this is a test case run for only 100 steps
         if self.test:
-            self.total_steps = 150
+            self.total_steps = 100
 
         self.dt = self.DT / self.w_ci # self.DT * self.t_ci
 

diff --git a/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py b/Examples/Tests/ohm_solver_ion_beam_instability/PICMI_inputs.py
@@ -98,7 +98,7 @@ def __init__(self, test, dim, resonant, verbose):
         # if this is a test case run for only 25 cyclotron periods and use
         # fewer substeps to speed up the simulation
         if self.test:
-            self.LT = 35.0
+            self.LT = 30.0
 
         self.total_steps = int(np.ceil(self.LT / self.DT))
 

diff --git a/Regression/Checksum/benchmarks_json/Python_ohms_law_solver_magnetic_reconnection_2d.json b/Regression/Checksum/benchmarks_json/Python_ohms_law_solver_magnetic_reconnection_2d.json
@@ -1,18 +1,18 @@
 {
   "lev=0": {
-    "Bx": 1524.8787030963545,
+    "Bx": 1524.8789585554075,
     "By": 639.8314135126764,
-    "Bz": 7.426827947249349,
-    "Ex": 55861538.51676546,
-    "Ey": 72378555.53423825,
-    "Ez": 304141253.3165592
+    "Bz": 7.476000395458839,
+    "Ex": 56499935.20497879,
+    "Ey": 75068107.08471295,
+    "Ez": 309535670.18599623
   },
   "ions": {
-    "particle_momentum_x": 7.139858023357858e-15,
-    "particle_momentum_y": 7.137973106773408e-15,
-    "particle_momentum_z": 7.136961183942627e-15,
-    "particle_position_x": 11170688.447856156,
-    "particle_position_y": 5585378.664524556,
+    "particle_momentum_x": 7.14295522755638e-15,
+    "particle_momentum_y": 7.1380780610425e-15,
+    "particle_momentum_z": 7.141134469227045e-15,
+    "particle_position_x": 11170689.203149632,
+    "particle_position_y": 5585328.083196239,
     "particle_weight": 9.036667901693183e+18
   }
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
@@ -117,13 +117,22 @@ public:
         amrex::IntVect ng, std::optional<bool> nodal_sync);
 
     void BfieldEvolveRK (
-            amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
-            amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
-            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::IntVect ng, std::optional<bool> nodal_sync);
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
+        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::IntVect ng, std::optional<bool> nodal_sync);
+
+    void BfieldPush (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
+        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::IntVect ng, std::optional<bool> nodal_sync);
 
     /**
      * \brief

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
@@ -563,56 +563,32 @@ void HybridPICModel::BfieldEvolveRK (
         amrex::Real const dt, int lev, DtType dt_type,
         IntVect ng, std::optional<bool> nodal_sync )
 {
-    // make copies of the B-field multifabs at t = n
+    // Make copies of the B-field multifabs at t = n and create multifabs for
+    // each direction to store the Runge-Kutta intermediate terms. Each
+    // multifab has 2 components for the different terms that need to be stored.
     std::array< MultiFab, 3 > B_old;
-    B_old[0] = MultiFab(
-        Bfield[lev][0]->boxArray(), Bfield[lev][0]->DistributionMap(), 1,
-        Bfield[lev][0]->nGrowVect()
-    );
-    MultiFab::Copy(B_old[0], *Bfield[lev][0], 0, 0, 1, ng);
-    B_old[1] = MultiFab(
-        Bfield[lev][1]->boxArray(), Bfield[lev][1]->DistributionMap(), 1,
-        Bfield[lev][1]->nGrowVect()
-    );
-    MultiFab::Copy(B_old[1], *Bfield[lev][1], 0, 0, 1, ng);
-    B_old[2] = MultiFab(
-        Bfield[lev][2]->boxArray(), Bfield[lev][2]->DistributionMap(), 1,
-        Bfield[lev][2]->nGrowVect()
-    );
-    MultiFab::Copy(B_old[2], *Bfield[lev][2], 0, 0, 1, ng);
-
-    // Create multifabs for each direction to store the Runge-Kutta
-    // intermediate terms. Each multifab has 2 components for the different
-    // terms that need to be stored.
-    std::array< amrex::MultiFab, 3 > K;
-    K[0] = amrex::MultiFab(
-        Bfield[lev][0]->boxArray(), Bfield[lev][0]->DistributionMap(), 2,
-        Bfield[lev][0]->nGrowVect()
-    );
-    K[0].setVal(0.0);
-    K[1] = amrex::MultiFab(
-        Bfield[lev][1]->boxArray(), Bfield[lev][1]->DistributionMap(), 2,
-        Bfield[lev][1]->nGrowVect()
-    );
-    K[1].setVal(0.0);
-    K[2] = amrex::MultiFab(
-        Bfield[lev][2]->boxArray(), Bfield[lev][2]->DistributionMap(), 2,
-        Bfield[lev][2]->nGrowVect()
-    );
-    K[2].setVal(0.0);
+    std::array< MultiFab, 3 > K;
+    for (int ii = 0; ii < 3; ii++)
+    {
+        B_old[ii] = MultiFab(
+            Bfield[lev][ii]->boxArray(), Bfield[lev][ii]->DistributionMap(), 1,
+            Bfield[lev][ii]->nGrowVect()
+        );
+        MultiFab::Copy(B_old[ii], *Bfield[lev][ii], 0, 0, 1, ng);
 
-    auto& warpx = WarpX::GetInstance();
+        K[ii] = MultiFab(
+            Bfield[lev][ii]->boxArray(), Bfield[lev][ii]->DistributionMap(), 2,
+            Bfield[lev][ii]->nGrowVect()
+        );
+        K[ii].setVal(0.0);
+    }
 
     // The Runge-Kutta scheme begins here.
     // Step 1:
-    // 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);
-    warpx.FillBoundaryE(ng, nodal_sync);
-    // Push forward the B-field using Faraday's law
-    warpx.EvolveB(0.5_rt*dt, dt_type);
-    warpx.FillBoundaryB(ng, nodal_sync);
+    BfieldPush(
+        Bfield, Efield, Jfield, rhofield, edge_lengths,
+        0.5_rt*dt, dt_type, ng, nodal_sync
+    );
 
     // The Bfield is now given by:
     // B_new = B_old + 0.5 * dt * [-curl x E(B_old)] = B_old + 0.5 * dt * K0.
@@ -625,14 +601,10 @@ void HybridPICModel::BfieldEvolveRK (
     }
 
     // Step 2:
-    // 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);
-    warpx.FillBoundaryE(ng, nodal_sync);
-    // Push forward the B-field using Faraday's law
-    warpx.EvolveB(0.5_rt*dt, dt_type);
-    warpx.FillBoundaryB(ng, nodal_sync);
+    BfieldPush(
+        Bfield, Efield, Jfield, rhofield, edge_lengths,
+        0.5_rt*dt, dt_type, ng, nodal_sync
+    );
 
     // The Bfield is now given by:
     // B_new = B_old + 0.5 * dt * K0 + 0.5 * dt * [-curl x E(B_old + 0.5 * dt * K1)]
@@ -649,18 +621,13 @@ void HybridPICModel::BfieldEvolveRK (
     }
 
     // Step 3:
-    // Calculate J = curl x B / mu0
-    CalculateCurrentAmpere(Bfield, edge_lengths);
-    // Calculate the E-field form Ohm's law
-    HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
-    warpx.FillBoundaryE(ng, nodal_sync);
-    // Push forward the B-field using Faraday's law.
-    warpx.EvolveB(dt, dt_type);
-    warpx.FillBoundaryB(ng, nodal_sync);
+    BfieldPush(
+        Bfield, Efield, Jfield, rhofield, edge_lengths,
+        dt, dt_type, ng, nodal_sync
+    );
 
     // The Bfield is now given by:
-    // B_new = B_old + 0.5 * dt * K1
-    //         + dt * [-curl  x E(B_old + 0.5 * dt * K1)]
+    // B_new = B_old + 0.5 * dt * K1 + dt * [-curl  x E(B_old + 0.5 * dt * K1)]
     //       = B_old + 0.5 * dt * K1 + dt * K2
     for (int ii = 0; ii < 3; ii++)
     {
@@ -670,14 +637,10 @@ void HybridPICModel::BfieldEvolveRK (
     }
 
     // Step 4:
-    // Calculate J = curl x B / mu0
-    CalculateCurrentAmpere(Bfield, edge_lengths);
-    // Calculate the E-field form Ohm's law
-    HybridPICSolveE(Efield, Jfield, Bfield, rhofield, edge_lengths, true);
-    warpx.FillBoundaryE(ng, nodal_sync);
-    // Push forward the B-field using Faraday's law.
-    warpx.EvolveB(0.5_rt*dt, dt_type);
-    warpx.FillBoundaryB(ng, nodal_sync);
+    BfieldPush(
+        Bfield, Efield, Jfield, rhofield, edge_lengths,
+        0.5_rt*dt, dt_type, ng, nodal_sync
+    );
 
     // The Bfield is now given by:
     // B_new = B_old + dt * K2 + 0.5 * dt * [-curl x E(B_old + dt * K2)]
@@ -703,3 +666,24 @@ void HybridPICModel::BfieldEvolveRK (
         );
     }
 }
+
+void HybridPICModel::BfieldPush (
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Bfield,
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>>& Efield,
+        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 const 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);
+    warpx.FillBoundaryE(ng, nodal_sync);
+    // Push forward the B-field using Faraday's law
+    warpx.EvolveB(dt, dt_type);
+    warpx.FillBoundaryB(ng, nodal_sync);
+}
diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
@@ -94,14 +94,11 @@ void WarpX::HybridPICEvolveFields ()
     // momentum equation
     for (int sub_step = 0; sub_step < sub_steps; sub_step++)
     {
-        // m_hybrid_pic_model->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
-        // m_hybrid_pic_model->HybridPICSolveE(
-        //     Efield_fp, current_fp_temp, Bfield_fp, rho_fp_temp, m_edge_lengths,
-        //     true
+        // m_hybrid_pic_model->BfieldPush(
+        //     Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp, m_edge_lengths,
+        //     0.5_rt/sub_steps*dt[0], DtType::FirstHalf,
+        //     guard_cells.ng_FieldSolver, WarpX::sync_nodal_points
         // );
-        // FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
-        // EvolveB(0.5_rt / sub_steps * dt[0], DtType::FirstHalf);
-        // FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
 
         m_hybrid_pic_model->BfieldEvolveRK(
             Bfield_fp, Efield_fp, current_fp_temp, rho_fp_temp,
@@ -129,14 +126,11 @@ void WarpX::HybridPICEvolveFields ()
     // 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->CalculateCurrentAmpere(Bfield_fp, m_edge_lengths);
-        // m_hybrid_pic_model->HybridPICSolveE(
-        //     Efield_fp, current_fp, Bfield_fp, rho_fp_temp, m_edge_lengths,
-        //     true
+        // m_hybrid_pic_model->BfieldPush(
+        //     Bfield_fp, Efield_fp, current_fp, rho_fp_temp, m_edge_lengths,
+        //     0.5_rt/sub_steps*dt[0], DtType::SecondHalf,
+        //     guard_cells.ng_FieldSolver, WarpX::sync_nodal_points
         // );
-        // FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
-        // EvolveB(0.5_rt / sub_steps * dt[0], DtType::SecondHalf);
-        // FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
 
         m_hybrid_pic_model->BfieldEvolveRK(
             Bfield_fp, Efield_fp, current_fp, rho_fp_temp,