diff --git a/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py b/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
index 33c9eef071b..18ac5ee9736 100755
--- a/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
+++ b/Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
@@ -76,15 +76,30 @@
 E_decay = 0.0918984*MeV_to_Joule # Energy released during Be -> 2*alpha
 E_fusion_total = E_fusion + E_decay # Energy released during p + B -> 3*alpha
 
+## Checks whether this is the 2D or the 3D test
+is_2D = "2D" in sys.argv[1]
+
 ## Some numerical parameters for this test
 size_x = 8
-size_y = 8
+if is_2D:
+    size_y = 1
+else:
+    size_y = 8
 size_z = 16
 dV_total = size_x*size_y*size_z # Total simulation volume
 # Volume of a slice corresponding to a single cell in the z direction. In tests 1 and 2, all the
 # particles of a given species in the same slice have the exact same momentum
 dV_slice = size_x*size_y
-dt = 1./(scc.c*np.sqrt(3.))
+if is_2D:
+    dt = 1./(scc.c*np.sqrt(2.))
+    yt_z_string = "particle_position_y"
+    nppcell_1 = 10000*8
+    nppcell_2 = 900*8
+else:
+    dt = 1./(scc.c*np.sqrt(3.))
+    yt_z_string = "particle_position_z"
+    nppcell_1 = 10000
+    nppcell_2 = 900
 # In test 1 and 2, the energy in cells number i (in z direction) is typically Energy_step * i**2
 Energy_step = 22.*keV_to_Joule
 
@@ -547,10 +562,10 @@ def check_initial_energy2(data):
         ## Tolerance is quite high below because we don't have a lot of alphas to produce good
         ## statistics and an event like alpha1 emitted exactly in direction of proton & alpha2
         ## emitted exactly in direction opposite to Beryllium is somewhat rare.
-        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=1.e-1))
-        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, rtol=1.e-1))
-        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=1.e-1))
-        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, rtol=1.e-1))
+        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=1.1e-1))
+        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, rtol=1.1e-1))
+        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=1.1e-1))
+        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, rtol=1.1e-1))
 
 def check_xy_isotropy(data):
     ## Checks that the alpha particles are emitted isotropically in x and y
@@ -624,7 +639,7 @@ def specific_check1(data):
     check_isotropy(data, relative_tolerance = 3.e-2)
     expected_fusion_number = check_macroparticle_number(data,
                                                         fusion_probability_target_value = 0.002,
-                                                        num_pair_per_cell = 10000)
+                                                        num_pair_per_cell = nppcell_1)
     E_com = compute_E_com1(data)
     check_alpha_yield(data, expected_fusion_number, E_com, proton_density = 1.,
                                                            boron_density = 1.)
@@ -635,7 +650,7 @@ def specific_check2(data):
     ## Only 900 particles pairs per cell here because we ignore the 10% of protons that are at rest
     expected_fusion_number = check_macroparticle_number(data,
                                                         fusion_probability_target_value = 0.02,
-                                                        num_pair_per_cell = 900)
+                                                        num_pair_per_cell = nppcell_2)
     E_com = compute_E_com2(data)
     check_alpha_yield(data, expected_fusion_number, E_com, proton_density = 1.e20,
                                                            boron_density = 1.e26)
diff --git a/Examples/Modules/nuclear_fusion/inputs_2d b/Examples/Modules/nuclear_fusion/inputs_2d
new file mode 100644
index 00000000000..aaff599105a
--- /dev/null
+++ b/Examples/Modules/nuclear_fusion/inputs_2d
@@ -0,0 +1,196 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+## With these parameters, each cell has a size of exactly 1 by 1
+max_step = 1
+amr.n_cell = 8 16
+amr.max_grid_size = 8
+amr.blocking_factor = 8
+amr.max_level = 0
+geometry.coord_sys   = 0
+geometry.prob_lo     = 0.    0.
+geometry.prob_hi     = 8.   16.
+
+#################################
+###### Boundary Condition #######
+#################################
+boundary.field_lo = periodic periodic
+boundary.field_hi = periodic periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_ics = 1
+warpx.verbose = 1
+warpx.cfl = 1.0
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = proton1 boron1 alpha1 proton2 boron2 alpha2 proton3 boron3 alpha3
+                          proton4 boron4 alpha4 proton5 boron5 alpha5
+
+my_constants.m_b11 =  10.9298*m_p # Boron 11 mass
+my_constants.m_reduced = m_p*m_b11/(m_p+m_b11)
+my_constants.keV_to_J = 1.e3*q_e
+my_constants.Energy_step = 22. * keV_to_J
+
+proton1.species_type = proton
+proton1.injection_style = "NRandomPerCell"
+proton1.num_particles_per_cell = 80000
+proton1.profile = constant
+proton1.density = 1.
+proton1.momentum_distribution_type = "parse_momentum_function"
+proton1.momentum_function_ux(x,y,z) = 0.
+proton1.momentum_function_uy(x,y,z) = 0.
+## Thanks to the floor, all particles in the same cell have the exact same momentum
+proton1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_p*clight)
+proton1.do_not_push = 1
+
+boron1.species_type = boron11
+boron1.injection_style = "NRandomPerCell"
+boron1.num_particles_per_cell = 80000
+boron1.profile = constant
+boron1.density = 1.
+boron1.momentum_distribution_type = "parse_momentum_function"
+boron1.momentum_function_ux(x,y,z) = 0.
+boron1.momentum_function_uy(x,y,z) = 0.
+## Thanks to the floor, all particles in the same cell have the exact same momentum
+boron1.momentum_function_uz(x,y,z) = -sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_b11*clight)
+boron1.do_not_push = 1
+
+alpha1.species_type = alpha
+alpha1.do_not_push = 1
+
+my_constants.background_dens = 1.e26
+my_constants.beam_dens = 1.e20
+
+proton2.species_type = proton
+proton2.injection_style = "NRandomPerCell"
+proton2.num_particles_per_cell = 8000
+proton2.profile = "parse_density_function"
+## A tenth of the macroparticles in each cell is made of immobile high-density background protons.
+## The other nine tenths are made of fast low-density beam protons.
+proton2.density_function(x,y,z) = if(x - floor(x) < 0.1, 10.*background_dens, 10./9.*beam_dens)
+proton2.momentum_distribution_type = "parse_momentum_function"
+proton2.momentum_function_ux(x,y,z) = 0.
+proton2.momentum_function_uy(x,y,z) = 0.
+proton2.momentum_function_uz(x,y,z) = "if(x - floor(x) < 0.1,
+                                          0., sqrt(2*m_p*Energy_step*(floor(z)**2))/(m_p*clight))"
+proton2.do_not_push = 1
+
+boron2.species_type = boron11
+boron2.injection_style = "NRandomPerCell"
+boron2.num_particles_per_cell = 800
+boron2.profile = constant
+boron2.density = background_dens
+boron2.momentum_distribution_type = "constant"
+boron2.do_not_push = 1
+
+alpha2.species_type = alpha
+alpha2.do_not_push = 1
+
+my_constants.temperature = 44. * keV_to_J
+
+proton3.species_type = proton
+proton3.injection_style = "NRandomPerCell"
+proton3.num_particles_per_cell = 4800
+proton3.profile = constant
+proton3.density = 1.e28
+proton3.momentum_distribution_type = "maxwell_boltzmann"
+proton3.theta = temperature/(m_p*clight**2)
+proton3.do_not_push = 1
+
+boron3.species_type = boron11
+boron3.injection_style = "NRandomPerCell"
+boron3.num_particles_per_cell = 8000
+boron3.profile = constant
+boron3.density = 5.e28
+boron3.momentum_distribution_type = "maxwell_boltzmann"
+boron3.theta = temperature/(m_b11*clight**2)
+boron3.do_not_push = 1
+
+alpha3.species_type = alpha
+alpha3.do_not_push = 1
+
+my_constants.proton4_energy = 550*keV_to_J
+
+proton4.species_type = proton
+proton4.injection_style = "NRandomPerCell"
+proton4.num_particles_per_cell = 800
+proton4.profile = "constant"
+proton4.density = 1.e35
+proton4.momentum_distribution_type = "constant"
+proton4.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton4.do_not_push = 1
+
+boron4.species_type = boron11
+boron4.injection_style = "NRandomPerCell"
+boron4.num_particles_per_cell = 800
+boron4.profile = constant
+boron4.density = 1.
+boron4.momentum_distribution_type = "constant"
+boron4.do_not_push = 1
+
+alpha4.species_type = alpha
+alpha4.do_not_push = 1
+
+proton5.species_type = proton
+proton5.injection_style = "NRandomPerCell"
+proton5.num_particles_per_cell = 800
+proton5.profile = "constant"
+proton5.density = 1.e35
+proton5.momentum_distribution_type = "constant"
+proton5.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton5.do_not_push = 1
+
+boron5.species_type = boron11
+boron5.injection_style = "NRandomPerCell"
+boron5.num_particles_per_cell = 800
+boron5.profile = constant
+boron5.density = 1.
+boron5.momentum_distribution_type = "constant"
+boron5.do_not_push = 1
+
+alpha5.species_type = alpha
+alpha5.do_not_push = 1
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = PBF1 PBF2 PBF3 PBF4 PBF5
+PBF1.species = proton1 boron1
+PBF1.product_species = alpha1
+PBF1.type = nuclearfusion
+PBF1.fusion_multiplier = 1.e50
+
+PBF2.species = proton2 boron2
+PBF2.product_species = alpha2
+PBF2.type = nuclearfusion
+PBF2.fusion_multiplier = 1.e15
+PBF2.fusion_probability_target_value = 0.02
+
+PBF3.species = proton3 boron3
+PBF3.product_species = alpha3
+PBF3.type = nuclearfusion
+PBF3.fusion_multiplier = 1.e15
+
+PBF4.species = proton4 boron4
+PBF4.product_species = alpha4
+PBF4.type = nuclearfusion
+PBF4.fusion_multiplier = 1.e21
+
+PBF5.species = proton5 boron5
+PBF5.product_species = alpha5
+PBF5.type = nuclearfusion
+PBF5.fusion_multiplier = 1.e21
+PBF5.fusion_probability_threshold = 1e123
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = rho
diff --git a/Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_2D.json b/Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_2D.json
new file mode 100644
index 00000000000..da391839ec6
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_2D.json
@@ -0,0 +1,145 @@
+{
+  "alpha1": {
+    "particle_cpu": 61116.0,
+    "particle_id": 1374948741264.0,
+    "particle_momentum_x": 4.688186295441445e-15,
+    "particle_momentum_y": 4.7208784582798594e-15,
+    "particle_momentum_z": 4.7027632008939775e-15,
+    "particle_position_x": 462511.1275563426,
+    "particle_position_y": 978224.9637545305,
+    "particle_weight": 2.982275352174261e-28
+  },
+  "alpha2": {
+    "particle_cpu": 54822.0,
+    "particle_id": 1218578759043.0,
+    "particle_momentum_x": 4.157834044969911e-15,
+    "particle_momentum_y": 4.10466425159846e-15,
+    "particle_momentum_z": 4.1936423888353025e-15,
+    "particle_position_x": 411510.6935411271,
+    "particle_position_y": 872044.4445281675,
+    "particle_weight": 2.9506572774666424e+18
+  },
+  "alpha3": {
+    "particle_cpu": 6450.0,
+    "particle_id": 158818703625.0,
+    "particle_momentum_x": 5.052515917170244e-16,
+    "particle_momentum_y": 4.950394146977726e-16,
+    "particle_momentum_z": 4.896418624715863e-16,
+    "particle_position_x": 54341.27656876971,
+    "particle_position_y": 104940.19225774865,
+    "particle_weight": 1.5087701480905637e+27
+  },
+  "alpha4": {
+    "particle_cpu": 307200.0,
+    "particle_id": 7432143974400.0,
+    "particle_momentum_x": 2.335145239158059e-14,
+    "particle_momentum_y": 2.348323075031563e-14,
+    "particle_momentum_z": 2.3477654672120372e-14,
+    "particle_position_x": 2458256.1229268126,
+    "particle_position_y": 4915207.19087103,
+    "particle_weight": 384.0000000000002
+  },
+  "alpha5": {
+    "particle_cpu": 307200.0,
+    "particle_id": 7620887654400.0,
+    "particle_momentum_x": 2.3402087717265325e-14,
+    "particle_momentum_y": 2.3415304986163325e-14,
+    "particle_momentum_z": 2.3501245625975474e-14,
+    "particle_position_x": 2457680.608473036,
+    "particle_position_y": 4915446.400371841,
+    "particle_weight": 3.839999999999998e-19
+  },
+  "boron1": {
+    "particle_cpu": 5120000.0,
+    "particle_id": 78643205120000.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.524243076150406e-13,
+    "particle_position_x": 40958301.591654316,
+    "particle_position_y": 81921136.14476715,
+    "particle_weight": 128.00000000000261
+  },
+  "boron2": {
+    "particle_cpu": 51200.0,
+    "particle_id": 1103626291200.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 409798.015821768,
+    "particle_position_y": 819270.9858143466,
+    "particle_weight": 1.2799999999016446e+28
+  },
+  "boron3": {
+    "particle_cpu": 512000.0,
+    "particle_id": 11639194112000.0,
+    "particle_momentum_x": 9.274252125064723e-15,
+    "particle_momentum_y": 9.27253250840722e-15,
+    "particle_momentum_z": 9.276115198445575e-15,
+    "particle_position_x": 4095810.08236732,
+    "particle_position_y": 8192326.250658387,
+    "particle_weight": 6.399497076617303e+30
+  },
+  "boron5": {
+    "particle_cpu": 51200.0,
+    "particle_id": 1208483891200.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 409516.8988702324,
+    "particle_position_y": 819323.2418191985,
+    "particle_weight": 127.99999999999999
+  },
+  "lev=0": {
+    "rho": 4.1015775274022543e+18
+  },
+  "proton1": {
+    "particle_cpu": 5120000.0,
+    "particle_id": 26214405120000.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.524243076150406e-13,
+    "particle_position_x": 40960140.72983792,
+    "particle_position_y": 81919772.69310114,
+    "particle_weight": 128.00000000000261
+  },
+  "proton2": {
+    "particle_cpu": 512000.0,
+    "particle_id": 10747904512000.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.371679089706304e-14,
+    "particle_position_x": 4095630.6981353555,
+    "particle_position_y": 8192073.551798361,
+    "particle_weight": 1.2885512789516445e+28
+  },
+  "proton3": {
+    "particle_cpu": 307170.0,
+    "particle_id": 6731210814406.0,
+    "particle_momentum_x": 1.6817613332220226e-15,
+    "particle_momentum_y": 1.6822707864863416e-15,
+    "particle_momentum_z": 1.6798193398170112e-15,
+    "particle_position_x": 2457375.867254385,
+    "particle_position_y": 4914855.8979505645,
+    "particle_weight": 1.2794970766173041e+30
+  },
+  "proton4": {
+    "particle_cpu": 51200.0,
+    "particle_id": 1192755251200.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.7581275870306353e-15,
+    "particle_position_x": 409842.79569749814,
+    "particle_position_y": 819245.2000431926,
+    "particle_weight": 1.2800000000000004e+37
+  },
+  "proton5": {
+    "particle_cpu": 51200.0,
+    "particle_id": 1203241011200.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.7581275870306353e-15,
+    "particle_position_x": 409709.9706207799,
+    "particle_position_y": 819158.891638082,
+    "particle_weight": 1.2800000000000004e+37
+  }
+}
\ No newline at end of file
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index a7d70499bfa..01fbf20a509 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -1803,6 +1803,21 @@ compileTest = 0
 doVis = 0
 analysisRoutine = Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
 
+[Proton_Boron_Fusion_2D]
+buildDir = .
+inputFile = Examples/Modules/nuclear_fusion/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString =
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 2
+compileTest = 0
+doVis = 0
+analysisRoutine = Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
+
 [Maxwell_Hybrid_QED_solver]
 buildDir = .
 inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d