Energy-preserving explicit PIC-MCC algorithm: electromagnetic

Examples/Tests/collision/CMakeLists.txt

@@ -96,7 +96,17 @@ add_warpx_test(
     2  # dims
     1  # nprocs
     inputs_test_2d_collisions_split_position_push_electrostatic  # inputs
-    "analysis_test_2d_collisions_split_position_push_electrostatic.py --path diags/reducedfiles/"  # analysis
+    "analysis_test_2d_collisions_split_position_push.py --path diags/reducedfiles/"  # analysis
+    "analysis_default_regression.py --path diags/diag1/"  # checksum
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_2d_collisions_split_position_push_electromagnetic  # name
+    2  # dims
+    1  # nprocs
+    inputs_test_2d_collisions_split_position_push_electromagnetic  # inputs
+    "analysis_test_2d_collisions_split_position_push.py --path diags/reducedfiles/"  # analysis
     "analysis_default_regression.py --path diags/diag1/"  # checksum
     OFF  # dependency
 )

Examples/Tests/collision/analysis_test_2d_collisions_split_position_push_electrostatic.py → Examples/Tests/collision/analysis_test_2d_collisions_split_position_push.py

@@ -43,7 +43,18 @@ def analyze(args: argparse.Namespace) -> None:
         [int(nppc) for nppc in num_particles_per_cell_list]
     )
     num_particles_per_cell = np.prod(num_particles_per_cell_array)
-    equipartition_value = 1 / (6 * num_particles_per_cell + 1)
+    electrostatic = (
+        True
+        if (
+            "warpx.do_electrostatic" in input_dict
+            and input_dict["warpx.do_electrostatic"] != "none"
+        )
+        else False
+    )
+    if electrostatic:
+        equipartition_value = 1 / (6 * num_particles_per_cell + 1)
+    else:
+        equipartition_value = 5 / (6 * num_particles_per_cell + 5)
 
     # normalize the field energy variation
     electron_temperature = float(input_dict["my_constants.Te"][0])
@@ -99,14 +110,14 @@ def analyze(args: argparse.Namespace) -> None:
 
     # verify the total energy conservation
     total_energy_error_norm = np.max(np.abs(total_energy_error))
-    relative_tolerance = 1e-5
+    relative_tolerance = 1e-5 if electrostatic else 5e-5
     if total_energy_error_norm >= relative_tolerance:
         raise ValueError(
             f"Total energy conservation failed with a maximum relative error of {total_energy_error_norm}"
         )
 
     # compare the final value of the field energy variation to the reference equipartition value
-    relative_tolerance = 1e-1
+    relative_tolerance = 1e-1 if electrostatic else 5e-1
     # average over the last 100 time steps to reduce noise
     field_energy_error_avg = np.mean(field_energy_error[-100:])
     if not np.isclose(

Examples/Tests/collision/inputs_base_2d_collisions_split_position_push

@@ -0,0 +1,95 @@
+# algo
+algo.field_gathering = energy-conserving
+algo.particle_shape = 2
+
+# amr
+amr.max_level = 0
+amr.n_cell = 16 16
+
+# boundary
+boundary.field_hi = periodic periodic
+boundary.field_lo = periodic periodic
+boundary.particle_hi = periodic periodic
+boundary.particle_lo = periodic periodic
+
+# collisions
+collision1.species = electrons protons
+collision2.species = electrons electrons
+collision3.species = protons protons
+collisions.collision_names = collision1 collision2 collision3
+
+# diagnostics
+diagnostics.diags_names = diag1
+diag1.diag_type = Full
+diag1.format = openpmd
+diag1.intervals = 2000
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+diag1.write_species = 1
+diag1.species = electrons protons
+diag1.electrons.variables = w x z ux uy uz
+diag1.protons.variables = w x z ux uy uz
+
+# electrons
+electrons.density = n0
+electrons.injection_style = "NUniformPerCell"
+electrons.momentum_distribution_type = gaussian
+electrons.num_particles_per_cell_each_dim = 8 8
+electrons.profile = constant
+electrons.species_type = electron
+electrons.ux_th = sqrt(Te*q_e/m_e)/clight
+electrons.uy_th = sqrt(Te*q_e/m_e)/clight
+electrons.uz_th = sqrt(Te*q_e/m_e)/clight
+electrons.xmax = 10*de0
+electrons.xmin = 0
+electrons.zmax = 10*de0
+electrons.zmin = 0
+
+# reduced diagnostics
+field_energy.intervals = 1
+field_energy.type = FieldEnergy
+particle_energy.intervals = 1
+particle_energy.type = ParticleEnergy
+
+# geometry
+geometry.dims = 2
+geometry.prob_hi = 10*de0 10*de0
+geometry.prob_lo = 0. 0.
+
+# interpolation
+interpolation.galerkin_scheme = 1
+
+# max_step
+max_step = 2000
+
+# my_constants
+my_constants.de0 = clight/wpe  # skin depth, m
+my_constants.dt = 0.1/wpe  # time step, s
+my_constants.n0 = 1e30  # plasma density, 1/m**3
+my_constants.Te = 100.  # electron temperature, eV
+my_constants.Ti = 100.  # ion temperature, eV
+my_constants.wpe = q_e*sqrt(n0/(m_e*epsilon0))  # electron plasma frequency, radians/s
+
+# particles
+particles.species_names = protons electrons
+
+# protons
+protons.density = n0
+protons.injection_style = "NUniformPerCell"
+protons.momentum_distribution_type = gaussian
+protons.num_particles_per_cell_each_dim = 8 8
+protons.profile = constant
+protons.species_type = proton
+protons.ux_th = sqrt(Ti*q_e/m_p)/clight
+protons.uy_th = sqrt(Ti*q_e/m_p)/clight
+protons.uz_th = sqrt(Ti*q_e/m_p)/clight
+protons.xmax = 10*de0
+protons.xmin = 0
+protons.zmax = 10*de0
+protons.zmin = 0
+
+# warpx
+warpx.const_dt = dt
+warpx.reduced_diags_names = field_energy particle_energy
+warpx.serialize_initial_conditions = 1
+warpx.use_filter = 0
+warpx.verbose = 1

Examples/Tests/collision/inputs_test_2d_collisions_split_position_push_electromagnetic

@@ -0,0 +1,5 @@
+# base input parameters
+FILE = inputs_base_2d_collisions_split_position_push
+
+# test input parameters
+algo.current_deposition = esirkepov

Examples/Tests/collision/inputs_test_2d_collisions_split_position_push_electrostatic

@@ -1,96 +1,5 @@
-# algo
-algo.field_gathering = energy-conserving
-algo.particle_shape = 2
+# base input parameters
+FILE = inputs_base_2d_collisions_split_position_push
 
-# amr
-amr.max_level = 0
-amr.n_cell = 16 16
-
-# boundary
-boundary.field_hi = periodic periodic
-boundary.field_lo = periodic periodic
-boundary.particle_hi = periodic periodic
-boundary.particle_lo = periodic periodic
-
-# collisions
-collision1.species = electrons protons
-collision2.species = electrons electrons
-collision3.species = protons protons
-collisions.collision_names = collision1 collision2 collision3
-
-# diagnostics
-diagnostics.diags_names = diag1
-diag1.diag_type = Full
-diag1.format = openpmd
-diag1.intervals = 1000
-diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
-diag1.write_species = 1
-diag1.species = electrons protons
-diag1.electrons.variables = w x z ux uy uz
-diag1.protons.variables = w x z ux uy uz
-
-# electrons
-electrons.density = n0
-electrons.injection_style = "NUniformPerCell"
-electrons.momentum_distribution_type = gaussian
-electrons.num_particles_per_cell_each_dim = 8 8
-electrons.profile = constant
-electrons.species_type = electron
-electrons.ux_th = sqrt(Te*q_e/m_e)/clight
-electrons.uy_th = sqrt(Te*q_e/m_e)/clight
-electrons.uz_th = sqrt(Te*q_e/m_e)/clight
-electrons.xmax = 10*de0
-electrons.xmin = 0
-electrons.zmax = 10*de0
-electrons.zmin = 0
-
-# reduced diagnostics
-field_energy.intervals = 1
-field_energy.type = FieldEnergy
-particle_energy.intervals = 1
-particle_energy.type = ParticleEnergy
-
-# geometry
-geometry.dims = 2
-geometry.prob_hi = 10*de0 10*de0
-geometry.prob_lo = 0. 0.
-
-# interpolation
-interpolation.galerkin_scheme = 1
-
-# max_step
-max_step = 2000
-
-# my_constants
-my_constants.de0 = clight/wpe  # skin depth, m
-my_constants.dt = 0.1/wpe  # time step, s
-my_constants.n0 = 1e30  # plasma density, 1/m**3
-my_constants.Te = 100.  # electron temperature, eV
-my_constants.Ti = 100.  # ion temperature, eV
-my_constants.wpe = q_e*sqrt(n0/(m_e*epsilon0))  # electron plasma frequency, radians/s
-
-# particles
-particles.species_names = protons electrons
-
-# protons
-protons.density = n0
-protons.injection_style = "NUniformPerCell"
-protons.momentum_distribution_type = gaussian
-protons.num_particles_per_cell_each_dim = 8 8
-protons.profile = constant
-protons.species_type = proton
-protons.ux_th = sqrt(Ti*q_e/m_p)/clight
-protons.uy_th = sqrt(Ti*q_e/m_p)/clight
-protons.uz_th = sqrt(Ti*q_e/m_p)/clight
-protons.xmax = 10*de0
-protons.xmin = 0
-protons.zmax = 10*de0
-protons.zmin = 0
-
-# warpx
-warpx.const_dt = dt
+# test input parameters
 warpx.do_electrostatic = labframe
-warpx.reduced_diags_names = field_energy particle_energy
-warpx.serialize_initial_conditions = 1
-warpx.use_filter = 0
-warpx.verbose = 1

Regression/Checksum/benchmarks_json/test_2d_collisions_split_position_push_electromagnetic.json

@@ -0,0 +1,32 @@
+{
+  "electrons": {
+    "particle_momentum_x": 4.947658759425637e-20,
+    "particle_momentum_y": 4.9853212329450156e-20,
+    "particle_momentum_z": 4.956761347800454e-20,
+    "particle_position_x": 0.00043383932455626456,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.00043499970649608455,
+    "particle_weight": 2823958725036869.5
+  },
+  "lev=0": {
+    "Bx": 52191.341644683744,
+    "By": 79896.69628240584,
+    "Bz": 52185.73368927454,
+    "Ex": 23175956467799.395,
+    "Ey": 16495293322825.637,
+    "Ez": 22906776598246.027,
+    "jx": 8.039958285202928e+18,
+    "jy": 6.834686318583636e+18,
+    "jz": 8.821193573182705e+18,
+    "rho": 67024055594.50891
+  },
+  "protons": {
+    "particle_momentum_x": 2.125626411655415e-18,
+    "particle_momentum_y": 2.1646633593783105e-18,
+    "particle_momentum_z": 2.135106915660081e-18,
+    "particle_position_x": 0.0004364883212149715,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.00043564917834130193,
+    "particle_weight": 2823958725036869.5
+  }
+}

Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp

@@ -400,6 +400,10 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir,
         for (std::size_t index = PIdx::nattribs; index < rnames.size(); ++index) {
             real_names[index] = rnames[index];
             real_flags[index] = tmp.h_redistribute_real_comp[index];
+            // Do not write the average momentum runtime components added in PhysicalParticleContainer
+            if (rnames[index] == "ux_avg" || rnames[index] == "uy_avg" || rnames[index] == "uz_avg") {
+                real_flags[index] = 0;
+            }
         }
 
         //   note: skip the mandatory AMREX_SPACEDIM positions for pure SoA

Source/Diagnostics/WarpXOpenPMD.cpp

@@ -648,9 +648,12 @@ for (const auto & particle_diag : particle_diags) {
     real_flags.resize(tmp.NumRealComps());
     for (size_t index = PIdx::nattribs; index < rn.size(); ++index) {
         real_flags[index] = tmp.h_redistribute_real_comp[index];
+        // Do not write the average momentum runtime components added in PhysicalParticleContainer
+        if (rn[index] == "ux_avg" || rn[index] == "uy_avg" || rn[index] == "uz_avg") {
+            real_flags[index] = 0;
+        }
     }
 
-
     // and the int components
     amrex::Vector<int> int_flags(tmp.NumIntComps());
     for (size_t i = 0; i < in.size(); ++i)