Add 2D tests for proton boron fusion

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
 
@@ -99,7 +114,7 @@ def add_existing_species_to_dict(yt_ad, data_dict, species_name, prefix, suffix)
     data_dict[prefix+"_w_"+suffix]   = yt_ad[species_name, "particle_weight"].v
     data_dict[prefix+"_id_"+suffix]  = yt_ad[species_name, "particle_id"].v
     data_dict[prefix+"_cpu_"+suffix] = yt_ad[species_name, "particle_cpu"].v
-    data_dict[prefix+"_z_"+suffix]   = yt_ad[species_name, "particle_position_z"].v
+    data_dict[prefix+"_z_"+suffix]   = yt_ad[species_name, yt_z_string].v
 
 def add_empty_species_to_dict(data_dict, species_name, prefix, suffix):
     data_dict[prefix+"_px_"+suffix]  = np.empty(0)
@@ -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)

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