Fixing bug that only applied last field in the list of external field… (#5690)

…s. Changed to add all fields.

This bug cropped up during code review and refactoring. I have
re-implemented the field accumulations for multiply defined external
fields. This was not caught by CI since the cylinder compression test
only has a single uniform compression field.

---------

Signed-off-by: S. Eric Clark <25495882+clarkse@users.noreply.github.com>

Author: clarkse

Examples/Tests/ohm_solver_cylinder_compression/inputs_test_3d_ohm_solver_cylinder_compression_picmi.py

@@ -95,8 +95,8 @@ def __init__(self, test, verbose):
 
             RM = np.sqrt(XM**2 + YM**2)
 
-            Ax_data = -0.5 * YM * self.dB
-            Ay_data = 0.5 * XM * self.dB
+            Ax_data = -0.25 * YM * self.dB
+            Ay_data = 0.25 * XM * self.dB
             Az_data = np.zeros_like(RM)
 
             # Write vector potential to file to exercise field loading via OpenPMD
@@ -261,11 +261,17 @@ def setup_run(self):
         #######################################################################
         # External Field definition. Sigmoid starting around 2.5 us
         A_ext = {
-            "uniform": {
+            "uniform_file": {
                 "read_from_file": True,
                 "path": "Afield.h5",
                 "A_time_external_function": "1/(1+exp(5*(1-(t-t0_ramp)*sqrt(2)/tau_ramp)))",
-            }
+            },
+            "uniform_analytical": {
+                "Ax_external_function": f"-0.25*y*{self.dB}",
+                "Ay_external_function": f"0.25*x*{self.dB}",
+                "Az_external_function": "0",
+                "A_time_external_function": "1/(1+exp(5*(1-(t-t0_ramp)*sqrt(2)/tau_ramp)))",
+            },
         }
 
         self.solver = picmi.HybridPICSolver(

Examples/Tests/ohm_solver_cylinder_compression/inputs_test_rz_ohm_solver_cylinder_compression_picmi.py

@@ -93,8 +93,8 @@ def __init__(self, test, verbose):
             Ar_data = np.zeros_like(RM)
             Az_data = np.zeros_like(RM)
 
-            # Zero padded outside of domain
-            At_data = 0.5 * RM * self.dB
+            # Only include half of the compression field here
+            At_data = 0.25 * RM * self.dB
 
             # Write vector potential to file to exercise field loading via
             series = io.Series("Afield.h5", io.Access.create)
@@ -255,11 +255,17 @@ def setup_run(self):
         #######################################################################
         # External Field definition. Sigmoid starting around 2.5 us
         A_ext = {
-            "uniform": {
+            "uniform_file": {
                 "read_from_file": True,
                 "path": "Afield.h5",
                 "A_time_external_function": "1/(1+exp(5*(1-(t-t0_ramp)*sqrt(2)/tau_ramp)))",
-            }
+            },
+            "uniform_analytical": {
+                "Ax_external_function": f"-0.25*y*{self.dB}",
+                "Ay_external_function": f"0.25*x*{self.dB}",
+                "Az_external_function": "0",
+                "A_time_external_function": "1/(1+exp(5*(1-(t-t0_ramp)*sqrt(2)/tau_ramp)))",
+            },
         }
 
         self.solver = picmi.HybridPICSolver(

Regression/Checksum/benchmarks_json/test_3d_ohm_solver_cylinder_compression_picmi.json

@@ -1,20 +1,20 @@
 {
   "lev=0": {
-    "Bx": 0.5334253070691776,
-    "By": 0.5318560243634998,
-    "Bz": 2252.108905639938,
-    "Ex": 10509838.331420777,
-    "Ey": 10512676.798857061,
-    "Ez": 8848.113963901804,
-    "rho": 384112.2912140536
+    "Bx": 0.5334251406746063,
+    "By": 0.5318559382056761,
+    "Bz": 2252.108858363565,
+    "Ex": 10509826.248254433,
+    "Ey": 10512665.210439455,
+    "Ez": 8848.110632517377,
+    "rho": 384112.2912140535
   },
   "ions": {
-    "particle_momentum_x": 2.161294367543349e-16,
-    "particle_momentum_y": 2.161870747294985e-16,
-    "particle_momentum_z": 2.0513400435256855e-16,
-    "particle_position_x": 769864.202585846,
-    "particle_position_y": 769908.6569812088,
-    "particle_position_z": 620721.1900338201,
+    "particle_momentum_x": 2.1612944051250247e-16,
+    "particle_momentum_y": 2.1618707843260381e-16,
+    "particle_momentum_z": 2.0513400435259629e-16,
+    "particle_position_x": 769864.2025705199,
+    "particle_position_y": 769908.6569665589,
+    "particle_position_z": 620721.1900338195,
     "particle_weight": 1.008292384042714e+19
   }
-}
+}

Regression/Checksum/benchmarks_json/test_rz_ohm_solver_cylinder_compression_picmi.json

@@ -1,20 +1,20 @@
 {
   "lev=0": {
-    "Br": 0.01190012639573578,
-    "Bt": 0.011313481779415917,
-    "Bz": 11.684908684984164,
-    "Er": 154581.58512851578,
-    "Et": 4798.276941148807,
-    "Ez": 193.22344271401872,
-    "rho": 7968.182346905438
+    "Br": 0.011900125915334049,
+    "Bt": 0.011313482081775999,
+    "Bz": 11.684907956225278,
+    "Er": 154581.64325434464,
+    "Et": 4797.794963571249,
+    "Ez": 193.22336541793413,
+    "rho": 7968.182346878508
   },
   "ions": {
-    "particle_momentum_x": 3.1125151786241107e-18,
-    "particle_momentum_y": 3.119385993047207e-18,
-    "particle_momentum_z": 3.0289560038617916e-18,
-    "particle_position_x": 13628.662686419664,
-    "particle_position_y": 2285.6952310457755,
-    "particle_theta": 115055.48935725243,
+    "particle_momentum_x": 3.112515238421571e-18,
+    "particle_momentum_y": 3.1193860531312334e-18,
+    "particle_momentum_z": 3.0289560038609835e-18,
+    "particle_position_x": 13628.662686094893,
+    "particle_position_y": 2285.6952310456554,
+    "particle_theta": 115055.48935714104,
     "particle_weight": 2.525423582445981e+18
   }
 }

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/ExternalVectorPotential.H

@@ -86,7 +86,7 @@ public:
     void CalculateExternalCurlA (std::string& coil_name);
 
     AMREX_FORCE_INLINE
-    void PopulateExternalFieldFromVectorPotential (
+    void AddExternalFieldFromVectorPotential (
         ablastr::fields::VectorField const& dstField,
         amrex::Real scale_factor,
         ablastr::fields::VectorField const& srcField,

Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/ExternalVectorPotential.cpp

@@ -46,7 +46,7 @@ ExternalVectorPotential::ReadParameters ()
     m_A_external.resize(m_nFields);
 
     m_A_ext_time_function.resize(m_nFields);
-    for (std::string & field_time : m_A_ext_time_function) {field_time = "1.0"; }
+    for (std::string & field_time : m_A_ext_time_function) { field_time = "1.0"; }
 
     m_A_external_time_parser.resize(m_nFields);
     m_A_time_scale.resize(m_nFields);
@@ -272,7 +272,7 @@ ExternalVectorPotential::CalculateExternalCurlA (std::string& coil_name)
 
 AMREX_FORCE_INLINE
 void
-ExternalVectorPotential::PopulateExternalFieldFromVectorPotential (
+ExternalVectorPotential::AddExternalFieldFromVectorPotential (
     ablastr::fields::VectorField const& dstField,
     amrex::Real scale_factor,
     ablastr::fields::VectorField const& srcField,
@@ -313,21 +313,21 @@ ExternalVectorPotential::PopulateExternalFieldFromVectorPotential (
                 // Skip field update in the embedded boundaries
                 if (update_Fx_arr && update_Fx_arr(i, j, k) == 0) { return; }
 
-                Fx(i,j,k) = scale_factor * Sx(i,j,k);
+                Fx(i,j,k) += scale_factor * Sx(i,j,k);
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 // Skip field update in the embedded boundaries
                 if (update_Fy_arr && update_Fy_arr(i, j, k) == 0) { return; }
 
-                Fy(i,j,k) = scale_factor * Sy(i,j,k);
+                Fy(i,j,k) += scale_factor * Sy(i,j,k);
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
                 // Skip field update in the embedded boundaries
                 if (update_Fz_arr && update_Fz_arr(i, j, k) == 0) { return; }
 
-                Fz(i,j,k) = scale_factor * Sz(i,j,k);
+                Fz(i,j,k) += scale_factor * Sz(i,j,k);
             }
         );
     }
@@ -339,12 +339,20 @@ ExternalVectorPotential::UpdateHybridExternalFields (const amrex::Real t, const
     using ablastr::fields::Direction;
     auto& warpx = WarpX::GetInstance();
 
-
     ablastr::fields::MultiLevelVectorField B_ext =
         warpx.m_fields.get_mr_levels_alldirs(FieldType::hybrid_B_fp_external, warpx.finestLevel());
     ablastr::fields::MultiLevelVectorField E_ext =
         warpx.m_fields.get_mr_levels_alldirs(FieldType::hybrid_E_fp_external, warpx.finestLevel());
 
+    // Zero E and B external fields prior to accumulating external fields
+    for (int lev = 0; lev <= warpx.finestLevel(); ++lev) {
+        for (int idir = 0; idir < 3; ++idir) {
+            B_ext[lev][Direction{idir}]->setVal(0.0_rt);
+            E_ext[lev][Direction{idir}]->setVal(0.0_rt);
+        }
+    }
+
+    // Iterate over external fields and add together with individual time functions.
     for (int i = 0; i < m_nFields; ++i) {
         const std::string Aext_field = m_field_names[i] + std::string{"_Aext"};
         const std::string curlAext_field = m_field_names[i] + std::string{"_curlAext"};
@@ -363,8 +371,8 @@ ExternalVectorPotential::UpdateHybridExternalFields (const amrex::Real t, const
             warpx.m_fields.get_mr_levels_alldirs(curlAext_field, warpx.finestLevel());
 
         for (int lev = 0; lev <= warpx.finestLevel(); ++lev) {
-            PopulateExternalFieldFromVectorPotential(E_ext[lev], scale_factor_E, A_ext[lev], warpx.GetEBUpdateEFlag()[lev]);
-            PopulateExternalFieldFromVectorPotential(B_ext[lev], scale_factor_B, curlA_ext[lev], warpx.GetEBUpdateBFlag()[lev]);
+            AddExternalFieldFromVectorPotential(E_ext[lev], scale_factor_E, A_ext[lev], warpx.GetEBUpdateEFlag()[lev]);
+            AddExternalFieldFromVectorPotential(B_ext[lev], scale_factor_B, curlA_ext[lev], warpx.GetEBUpdateBFlag()[lev]);
 
             for (int idir = 0; idir < 3; ++idir) {
                 E_ext[lev][Direction{idir}]->FillBoundary(warpx.Geom(lev).periodicity());