diff options
6 files changed, 290 insertions, 16 deletions
diff --git a/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py index f218df54d..b5f5da683 100755 --- a/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py +++ b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py @@ -6,6 +6,7 @@ # License: BSD-3-Clause-LBNL import os +import re import sys import yt @@ -66,15 +67,30 @@ barn_to_square_meter = 1.e-28 E_fusion = 17.5893*MeV_to_Joule # Energy released during the fusion reaction +## Checks whether this is the 2D or the 3D test +warpx_used_inputs = open('./warpx_used_inputs', 'r').read() +if re.search('geometry.dims = RZ', warpx_used_inputs): + is_RZ = True +else: + is_RZ = False + ## Some numerical parameters for this test size_x = 8 size_y = 8 size_z = 16 -dV_total = size_x*size_y*size_z # Total simulation volume +if is_RZ: + dV_slice = np.pi * size_x**2 + yt_z_string = "particle_position_y" + nppcell_1 = 10000*8 + nppcell_2 = 900*8 +else: + dV_slice = size_x*size_y + yt_z_string = "particle_position_z" + nppcell_1 = 10000 + nppcell_2 = 900 # 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.)) + # 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 @@ -89,7 +105,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) @@ -263,8 +279,11 @@ def expected_weight_com(E_com, reactant0_density, reactant1_density, dV, dt): def check_macroparticle_number(data, fusion_probability_target_value, num_pair_per_cell): ## Checks that the number of macroparticles is as expected for the first and second tests - ## The first slice 0 < z < 1 does not contribute to product species creation - numcells = dV_total - dV_slice + ## The first slice 0 < z < 1 does not contribute to alpha creation + if is_RZ: + numcells = size_x*(size_z-1) + else: + numcells = size_x*size_y*(size_z-1) ## In these tests, the fusion_multiplier is so high that the fusion probability per pair is ## equal to the parameter fusion_probability_target_value fusion_probability_per_pair = fusion_probability_target_value @@ -315,7 +334,7 @@ def compute_E_com2(data): p_reactant0_sq = 2.*mass[reactant_species[0]]*(Energy_step*np.arange(size_z)**2) return p_sq_reactant1_frame_to_E_COM_frame(p_reactant0_sq) -def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, reactant1_density): +def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, reactant1_density, dt): ## Checks that the fusion yield is as expected for the first and second tests. product_weight_theory = expected_weight_com(E_com/keV_to_Joule, reactant0_density, reactant1_density, dV_slice, dt) @@ -330,24 +349,25 @@ def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, r assert(np.all(is_close(product_weight_theory, product_weight_simulation, rtol = 5.*relative_std_weight))) -def specific_check1(data): - check_isotropy(data, relative_tolerance = 3.e-2) +def specific_check1(data, dt): + if not is_RZ: + 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_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1., - reactant1_density = 1.) + reactant1_density = 1., dt=dt) -def specific_check2(data): +def specific_check2(data, dt): check_xy_isotropy(data) ## Only 900 particles pairs per cell here because we ignore the 10% of reactants 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_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1.e20, - reactant1_density = 1.e26) + reactant1_density = 1.e26, dt=dt) def check_charge_conservation(rho_start, rho_end): assert(np.all(is_close(rho_start, rho_end, rtol=2.e-11))) @@ -359,6 +379,7 @@ def main(): ds_start = yt.load(filename_start) ad_end = ds_end.all_data() ad_start = ds_start.all_data() + dt = float(ds_end.current_time - ds_start.current_time) field_data_end = ds_end.covering_grid(level=0, left_edge=ds_end.domain_left_edge, dims=ds_end.domain_dimensions) field_data_start = ds_start.covering_grid(level=0, left_edge=ds_start.domain_left_edge, @@ -379,7 +400,7 @@ def main(): generic_check(data) # Checks that are specific to test number i - eval("specific_check"+str(i)+"(data)") + eval("specific_check"+str(i)+"(data, dt)") rho_start = field_data_start["rho"].to_ndarray() rho_end = field_data_end["rho"].to_ndarray() diff --git a/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz new file mode 100644 index 000000000..fb581c825 --- /dev/null +++ b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz @@ -0,0 +1,129 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +## With these parameters, each cell has a size of exactly 1 by 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.dims = RZ +geometry.prob_lo = 0. 0. +geometry.prob_hi = 8. 16. + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = none periodic +boundary.field_hi = pec periodic + +################################# +############ NUMERICS ########### +################################# +warpx.verbose = 1 +warpx.cfl = 1.0 + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = deuterium1 tritium1 helium1 neutron1 deuterium2 tritium2 helium2 neutron2 + +my_constants.m_deuterium = 2.01410177812*m_u +my_constants.m_tritium = 3.0160492779*m_u +my_constants.m_reduced = m_deuterium*m_tritium/(m_deuterium+m_tritium) +my_constants.keV_to_J = 1.e3*q_e +my_constants.Energy_step = 22. * keV_to_J + +deuterium1.species_type = deuterium +deuterium1.injection_style = "NRandomPerCell" +deuterium1.num_particles_per_cell = 80000 +deuterium1.profile = constant +deuterium1.density = 1. +deuterium1.momentum_distribution_type = parse_momentum_function +## Thanks to the floor, all particles in the same cell have the exact same momentum +deuterium1.momentum_function_ux(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight); if(x*x+y*y>0.0, -u*y/sqrt(x*x+y*y), 0.0)" # azimuthal velocity +deuterium1.momentum_function_uy(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight); if(x*x+y*y>0.0, u*x/sqrt(x*x+y*y), 0.0)" # azimuthal velocity +deuterium1.momentum_function_uz(x,y,z) = "0" +deuterium1.do_not_push = 1 +deuterium1.do_not_deposit = 1 + +tritium1.species_type = tritium +tritium1.injection_style = "NRandomPerCell" +tritium1.num_particles_per_cell = 80000 +tritium1.profile = constant +tritium1.density = 1. +tritium1.momentum_distribution_type = "parse_momentum_function" +## Thanks to the floor, all particles in the same cell have the exact same momentum +tritium1.momentum_function_ux(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight); if(x*x+y*y>0.0, u*y/sqrt(x*x+y*y), 0.0)" # counter-streaming azimuthal velocity +tritium1.momentum_function_uy(x,y,z) = "u = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight); if(x*x+y*y>0.0, -u*x/sqrt(x*x+y*y), 0.0)" # counter-streaming azimuthal velocity +tritium1.momentum_function_uz(x,y,z) = 0 +tritium1.do_not_push = 1 +tritium1.do_not_deposit = 1 + +helium1.species_type = helium4 +helium1.do_not_push = 1 +helium1.do_not_deposit = 1 + +neutron1.species_type = neutron +neutron1.do_not_push = 1 +neutron1.do_not_deposit = 1 + +my_constants.background_dens = 1.e26 +my_constants.beam_dens = 1.e20 + +deuterium2.species_type = deuterium +deuterium2.injection_style = "NRandomPerCell" +deuterium2.num_particles_per_cell = 8000 +deuterium2.profile = "parse_density_function" +## A tenth of the macroparticles in each cell is made of immobile high-density background deuteriums. +## The other nine tenths are made of fast low-density beam deuteriums. +deuterium2.density_function(x,y,z) = if(y - floor(y) < 0.1, 10.*background_dens, 10./9.*beam_dens) +deuterium2.momentum_distribution_type = "parse_momentum_function" +deuterium2.momentum_function_ux(x,y,z) = 0. +deuterium2.momentum_function_uy(x,y,z) = 0. +deuterium2.momentum_function_uz(x,y,z) = "if(y - floor(y) < 0.1, + 0., sqrt(2*m_deuterium*Energy_step*(floor(z)**2))/(m_deuterium*clight))" +deuterium2.do_not_push = 1 +deuterium2.do_not_deposit = 1 + +tritium2.species_type = tritium +tritium2.injection_style = "NRandomPerCell" +tritium2.num_particles_per_cell = 800 +tritium2.profile = constant +tritium2.density = background_dens +tritium2.momentum_distribution_type = "constant" +tritium2.do_not_push = 1 +tritium2.do_not_deposit = 1 + +helium2.species_type = helium4 +helium2.do_not_push = 1 +helium2.do_not_deposit = 1 + +neutron2.species_type = neutron +neutron2.do_not_push = 1 +neutron2.do_not_deposit = 1 + +################################# +############ COLLISION ########## +################################# +collisions.collision_names = DTF1 DTF2 + +DTF1.species = deuterium1 tritium1 +DTF1.product_species = helium1 neutron1 +DTF1.type = nuclearfusion +DTF1.fusion_multiplier = 1.e50 + +DTF2.species = deuterium2 tritium2 +DTF2.product_species = helium2 neutron2 +DTF2.type = nuclearfusion +DTF2.fusion_multiplier = 1.e15 +DTF2.fusion_probability_target_value = 0.02 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 1 +diag1.diag_type = Full +diag1.fields_to_plot = rho diff --git a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json index eec472698..2c6408468 100644 --- a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json +++ b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json @@ -74,4 +74,4 @@ "particle_position_z": 819255.8152412223, "particle_weight": 1.0239999999424347e+29 } -}
\ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json new file mode 100644 index 000000000..a5136dda6 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_RZ.json @@ -0,0 +1,77 @@ +{ + "deuterium1": { + "particle_momentum_x": 1.8388106511899905e-13, + "particle_momentum_y": 1.837868790009435e-13, + "particle_momentum_z": 0.0, + "particle_position_x": 40959919.499819286, + "particle_position_y": 81919224.48541151, + "particle_theta": 32166860.23003994, + "particle_weight": 3216.984554806547 + }, + "deuterium2": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 3.336364094249911e-14, + "particle_position_x": 4095908.9083809257, + "particle_position_y": 8192069.080030457, + "particle_theta": 3216444.348910214, + "particle_weight": 3.1898417901971444e+29 + }, + "helium1": { + "particle_momentum_x": 1.858124399143442e-15, + "particle_momentum_y": 1.876715110797694e-15, + "particle_momentum_z": 1.7098432207359157e-15, + "particle_position_x": 152920.23233108618, + "particle_position_y": 323733.9138644398, + "particle_theta": 120064.13771707338, + "particle_weight": 1.603083276067953e-27 + }, + "helium2": { + "particle_momentum_x": 1.5195006688950936e-15, + "particle_momentum_y": 1.52430083815551e-15, + "particle_momentum_z": 1.7654865863613367e-15, + "particle_position_x": 136867.63803188328, + "particle_position_y": 286903.30393175944, + "particle_theta": 107912.20520382549, + "particle_weight": 2.0862696876352987e+19 + }, + "lev=0": { + "rho": 0.0 + }, + "neutron1": { + "particle_momentum_x": 1.7160671487712845e-15, + "particle_momentum_y": 1.7154753069055672e-15, + "particle_momentum_z": 1.7098432207359157e-15, + "particle_position_x": 152920.23233108618, + "particle_position_y": 323733.9138644398, + "particle_theta": 120064.13771707338, + "particle_weight": 1.603083276067953e-27 + }, + "neutron2": { + "particle_momentum_x": 1.5195006688950936e-15, + "particle_momentum_y": 1.52430083815551e-15, + "particle_momentum_z": 1.5463311225724366e-15, + "particle_position_x": 136867.63803188328, + "particle_position_y": 286903.30393175944, + "particle_theta": 107912.20520382549, + "particle_weight": 2.0862696876352987e+19 + }, + "tritium1": { + "particle_momentum_x": 1.8384658063720362e-13, + "particle_momentum_y": 1.8381593257898129e-13, + "particle_momentum_z": 0.0, + "particle_position_x": 40961278.052658774, + "particle_position_y": 81919046.8061561, + "particle_theta": 32163925.891884565, + "particle_weight": 3217.0912552970394 + }, + "tritium2": { + "particle_momentum_x": 0.0, + "particle_momentum_y": 0.0, + "particle_momentum_z": 0.0, + "particle_position_x": 409793.9651940968, + "particle_position_y": 819237.3558155322, + "particle_theta": 321974.4557387621, + "particle_weight": 3.218514276139388e+29 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index e123efef2..e274b9058 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2291,6 +2291,22 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py +[Deuterium_Tritium_Fusion_RZ] +buildDir = . +inputFile = Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_rz +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 +dim = 2 +addToCompileString = USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS=RZ +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py + [Maxwell_Hybrid_QED_solver] buildDir = . inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index a807a4262..759462427 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -196,12 +196,36 @@ public: multiplier_ratio = max_N; } +#if (defined WARPX_DIM_RZ) + /* This momentum rotation is analogous to the one in ElasticCollisionPerez.H. */ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::n_rz_azimuthal_modes==1, + "RZ mode `warpx.n_rz_azimuthal_modes` must be 1 when using the binary nuclear fusion module."); + amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; + amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta]; +#endif + for (int k = 0; k < max_N; ++k) { // c1k : how many times the current particle of species 1 is paired with a particle // of species 2. Same for c2k. const int c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1; const int c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2; + +#if (defined WARPX_DIM_RZ) + /* In RZ geometry, macroparticles can collide with other macroparticles + * in the same *cylindrical* cell. For this reason, collisions between macroparticles + * are actually not local in space. In this case, the underlying assumption is that + * particles within the same cylindrical cell represent a cylindrically-symmetry + * momentum distribution function. Therefore, here, we temporarily rotate the + * momentum of one of the macroparticles in agreement with this cylindrical symmetry. + * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation; + * there is a corresponding assert statement at initialization.) */ + amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]]; + amrex::ParticleReal const u1xbuf = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta); + u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta); +#endif + SingleNuclearFusionEvent( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], @@ -211,6 +235,13 @@ public: m_probability_threshold, m_probability_target_value, m_fusion_type, engine); + +#if (defined WARPX_DIM_RZ) + amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta); + u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta); +#endif + p_pair_indices_1[pair_index] = I1[i1]; p_pair_indices_2[pair_index] = I2[i2]; ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; } |