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])

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 = 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
+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

Source/Evolve/WarpXEvolve.cpp

@@ -441,34 +441,40 @@ void WarpX::OneStep (
         }
         // electromagnetic solver
         else {
-            // perform particle collisions
-            ExecutePythonCallback("beforecollisions");
-            mypc->doCollisions(a_step, a_cur_time, a_dt);
-            ExecutePythonCallback("aftercollisions");
-
             // without mesh refinement
             if (finest_level == 0) {
                 // standard PIC loop
                 if (!m_JRhom) {
-                    OneStep_nosub(a_cur_time);
+                    OneStep_nosub(a_cur_time, a_dt, a_step);
                 }
                 // JRhom PIC loop
                 else {
+                    // perform particle collisions
+                    ExecutePythonCallback("beforecollisions");
+                    mypc->doCollisions(a_step, a_cur_time, a_dt);
+                    ExecutePythonCallback("aftercollisions");
+
                     OneStep_JRhom(a_cur_time);
                 }
             }
             // with mesh refinement
             else {
                 // without subcycling
                 if (!m_do_subcycling) {
-                    OneStep_nosub(a_cur_time);
+                    OneStep_nosub(a_cur_time, a_dt, a_step);
                 }
                 // with subcycling
                 else {
                     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
                         finest_level == 1,
                         "Subcycling not implemented with more than 1 mesh refinement level"
                     );
+
+                    // perform particle collisions
+                    ExecutePythonCallback("beforecollisions");
+                    mypc->doCollisions(a_step, a_cur_time, a_dt);
+                    ExecutePythonCallback("aftercollisions");
+
                     OneStep_sub1(a_cur_time);
                 }
             }
@@ -482,7 +488,11 @@ void WarpX::OneStep (
  * for the field advance and particle pusher.
  */
 void
-WarpX::OneStep_nosub (Real cur_time)
+WarpX::OneStep_nosub (
+    amrex::Real a_cur_time,
+    amrex::Real a_dt,
+    int a_step
+)
 {
     WARPX_PROFILE("WarpX::OneStep_nosub()");
 
@@ -494,7 +504,29 @@ WarpX::OneStep_nosub (Real cur_time)
     ExecutePythonCallback("particlescraper");
     ExecutePythonCallback("beforedeposition");
 
-    PushParticlesandDeposit(cur_time);
+    // push particles (half position and full momentum)
+    PushParticlesandDeposit(
+        a_cur_time,
+        /*skip_current=*/true,
+        PositionPushType::FirstHalf,
+        MomentumPushType::Full
+    );
+
+    // communicate particle data
+    mypc->Redistribute();
+
+    // perform particle collisions
+    ExecutePythonCallback("beforecollisions");
+    mypc->doCollisions(a_step, a_cur_time, a_dt);
+    ExecutePythonCallback("aftercollisions");
+
+    // push particles (half position)
+    PushParticlesandDeposit(
+        a_cur_time,
+        /*skip_current=*/false,
+        PositionPushType::SecondHalf,
+        MomentumPushType::None
+    );
 
     ExecutePythonCallback("afterdeposition");
 
@@ -521,7 +553,7 @@ WarpX::OneStep_nosub (Real cur_time)
             WarpX::Hybrid_QED_Push(dt);
             FillBoundaryE(guard_cells.ng_alloc_EB);
         }
-        PushPSATD(cur_time);
+        PushPSATD(a_cur_time);
 
         if (do_pml) {
             DampPML();
@@ -549,23 +581,23 @@ WarpX::OneStep_nosub (Real cur_time)
         FillBoundaryF(guard_cells.ng_FieldSolverF);
         FillBoundaryG(guard_cells.ng_FieldSolverG);
 
-        EvolveB(0.5_rt * dt[0], SubcyclingHalf::FirstHalf, cur_time); // We now have B^{n+1/2}
+        EvolveB(0.5_rt * dt[0], SubcyclingHalf::FirstHalf, a_cur_time); // We now have B^{n+1/2}
         FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
 
         if (m_em_solver_medium == MediumForEM::Vacuum) {
             // vacuum medium
-            EvolveE(dt[0], cur_time); // We now have E^{n+1}
+            EvolveE(dt[0], a_cur_time); // We now have E^{n+1}
         } else if (m_em_solver_medium == MediumForEM::Macroscopic) {
             // macroscopic medium
-            MacroscopicEvolveE(dt[0], cur_time); // We now have E^{n+1}
+            MacroscopicEvolveE(dt[0], a_cur_time); // We now have E^{n+1}
         } else {
             WARPX_ABORT_WITH_MESSAGE("Medium for EM is unknown");
         }
         FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);
 
         EvolveF(0.5_rt * dt[0], /*rho_comp=*/1);
         EvolveG(0.5_rt * dt[0]);
-        EvolveB(0.5_rt * dt[0], SubcyclingHalf::SecondHalf, cur_time + 0.5_rt * dt[0]); // We now have B^{n+1}
+        EvolveB(0.5_rt * dt[0], SubcyclingHalf::SecondHalf, a_cur_time + 0.5_rt * dt[0]); // We now have B^{n+1}
 
         if (do_pml) {
             DampPML();

Source/Particles/MultiParticleContainer.H

@@ -69,7 +69,10 @@ class MultiParticleContainer
 
 public:
 
-    MultiParticleContainer (amrex::AmrCore* amr_core);
+    MultiParticleContainer (
+        amrex::AmrCore* amr_core,
+        bool const collisions_split_position_push
+    );
 
     ~MultiParticleContainer() = default;