From b3ba07df16de4e9736b01a8cc7024740714dcad6 Mon Sep 17 00:00:00 2001 From: David Grote Date: Sat, 4 Nov 2023 00:03:20 -0700 Subject: [PATCH] Allow negative mean velocity with NFluxPerCell (#4397) * Allow negative drift velocity with Gaussian flux injection * Simplify by using sign of u_m * Add negative u_m test to flux_injection * Update FluxInjection3D benchmark * Save abs(u_m) as a temporary const --- .../analysis_flux_injection_3d.py | 48 ++++++++++++++----- Examples/Tests/flux_injection/inputs_3d | 33 ++++++++++++- .../benchmarks_json/FluxInjection3D.json | 46 ++++++++++++------ Source/Initialization/InjectorMomentum.H | 36 +++++++------- 4 files changed, 115 insertions(+), 48 deletions(-) diff --git a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py index 048ef70f9cc..0998de25a7b 100755 --- a/Examples/Tests/flux_injection/analysis_flux_injection_3d.py +++ b/Examples/Tests/flux_injection/analysis_flux_injection_3d.py @@ -57,11 +57,8 @@ def gaussian_dist(u, u_th): return 1./((2*np.pi)**.5*u_th) * np.exp(-u**2/(2*u_th**2) ) def gaussian_flux_dist(u, u_th, u_m): - au_m = np.abs(u_m) - normalization_factor = u_th**2 * np.exp(-au_m**2/(2*u_th**2)) + (np.pi/2)**.5*au_m*u_th * (1 + erf(au_m/(2**.5*u_th))) - result = 1./normalization_factor * np.where( u>0, u * np.exp(-(u-au_m)**2/(2*u_th**2)), 0 ) - if u_m < 0.: - result = result[::-1] + normalization_factor = u_th**2 * np.exp(-u_m**2/(2*u_th**2)) + (np.pi/2)**.5*u_m*u_th * (1 + erf(u_m/(2**.5*u_th))) + result = 1./normalization_factor * np.where( u>0, u * np.exp(-(u-u_m)**2/(2*u_th**2)), 0 ) return result def compare_gaussian(u, w, u_th, label=''): @@ -84,10 +81,10 @@ def compare_gaussian_flux(u, w, u_th, u_m, label=''): # Load data and perform check -plt.figure(figsize=(5,7)) +plt.figure(figsize=(8,7)) -plt.subplot(211) -plt.title('Electrons') +plt.subplot(221) +plt.title('Electrons u_m=0.07') ux = ad['electron','particle_momentum_x'].to_ndarray()/(m_e*c) uy = ad['electron','particle_momentum_y'].to_ndarray()/(m_e*c) @@ -97,21 +94,46 @@ def compare_gaussian_flux(u, w, u_th, u_m, label=''): compare_gaussian(ux, w, u_th=0.1, label='u_x') compare_gaussian_flux(uy, w, u_th=0.1, u_m=0.07, label='u_y') compare_gaussian(uz, w, u_th=0.1, label='u_z') -plt.legend(loc=0) -plt.subplot(212) -plt.title('Protons') +plt.subplot(223) +plt.title('Protons u_m=0.05') ux = ad['proton','particle_momentum_x'].to_ndarray()/(m_p*c) uy = ad['proton','particle_momentum_y'].to_ndarray()/(m_p*c) uz = ad['proton','particle_momentum_z'].to_ndarray()/(m_p*c) w = ad['proton', 'particle_weight'].to_ndarray() -compare_gaussian_flux(ux, w, u_th=0.1, u_m=-0.05, label='u_x') +compare_gaussian_flux(-ux, w, u_th=0.1, u_m=0.05, label='u_x') compare_gaussian(uy, w, u_th=0.1, label='u_y') compare_gaussian(uz, w, u_th=0.1, label='u_z') -plt.legend(loc=0) +plt.subplot(222) +plt.title('Electrons u_m=-0.07') + +ux = ad['electron_negative','particle_momentum_x'].to_ndarray()/(m_e*c) +uy = ad['electron_negative','particle_momentum_y'].to_ndarray()/(m_e*c) +uz = ad['electron_negative','particle_momentum_z'].to_ndarray()/(m_e*c) +w = ad['electron_negative', 'particle_weight'].to_ndarray() + +compare_gaussian(ux, w, u_th=0.1, label='u_x') +compare_gaussian(uy, w, u_th=0.1, label='u_y') +compare_gaussian_flux(uz, w, u_th=0.1, u_m=-0.07, label='u_z') +plt.legend(loc=(1.02, 0.5)) + +plt.subplot(224) +plt.title('Protons u_m=-0.05') + +ux = ad['proton_negative','particle_momentum_x'].to_ndarray()/(m_p*c) +uy = ad['proton_negative','particle_momentum_y'].to_ndarray()/(m_p*c) +uz = ad['proton_negative','particle_momentum_z'].to_ndarray()/(m_p*c) +w = ad['proton_negative', 'particle_weight'].to_ndarray() + +compare_gaussian(ux, w, u_th=0.1, label='u_x') +compare_gaussian(uy, w, u_th=0.1, label='u_y') +compare_gaussian_flux(-uz, w, u_th=0.1, u_m=-0.05, label='u_z') +#plt.legend(loc=0) + +plt.tight_layout() plt.savefig('Distribution.png') # Verify checksum diff --git a/Examples/Tests/flux_injection/inputs_3d b/Examples/Tests/flux_injection/inputs_3d index 7d277e7c292..80ee23feff9 100644 --- a/Examples/Tests/flux_injection/inputs_3d +++ b/Examples/Tests/flux_injection/inputs_3d @@ -28,7 +28,7 @@ boundary.field_lo = periodic periodic periodic boundary.field_hi = periodic periodic periodic # particles -particles.species_names = electron proton +particles.species_names = electron proton electron_negative proton_negative algo.particle_shape = 3 electron.charge = -q_e @@ -61,6 +61,37 @@ proton.ux_m = 0.05 proton.uy_th = 0.1 proton.uz_th = 0.1 +# "negative" is negative u_m +electron_negative.charge = -q_e +electron_negative.mass = m_e +electron_negative.injection_style = NFluxPerCell +electron_negative.num_particles_per_cell = 100 +electron_negative.surface_flux_pos = -4. +electron_negative.flux_normal_axis = z +electron_negative.flux_direction = +1 +electron_negative.flux_profile = parse_flux_function +electron_negative.flux_function(x,y,z,t) = "1." +electron_negative.momentum_distribution_type = gaussianflux +electron_negative.ux_th = 0.1 +electron_negative.uy_th = 0.1 +electron_negative.uz_th = 0.1 +electron_negative.uz_m = -0.07 + +proton_negative.charge = +q_e +proton_negative.mass = m_p +proton_negative.injection_style = NFluxPerCell +proton_negative.num_particles_per_cell = 100 +proton_negative.surface_flux_pos = 4. +proton_negative.flux_normal_axis = z +proton_negative.flux_direction = -1 +proton_negative.flux_profile = constant +proton_negative.flux = 1. +proton_negative.momentum_distribution_type = gaussianflux +proton_negative.ux_th = 0.1 +proton_negative.uy_th = 0.1 +proton_negative.uz_th = 0.1 +proton_negative.uz_m = -0.05 + # Diagnostics diagnostics.diags_names = diag1 diag1.intervals = 1000 diff --git a/Regression/Checksum/benchmarks_json/FluxInjection3D.json b/Regression/Checksum/benchmarks_json/FluxInjection3D.json index b2b3733737e..aa30f988f68 100644 --- a/Regression/Checksum/benchmarks_json/FluxInjection3D.json +++ b/Regression/Checksum/benchmarks_json/FluxInjection3D.json @@ -1,21 +1,39 @@ { - "electron": { - "particle_momentum_x": 1.1192116199394354e-18, - "particle_momentum_y": 2.238114590066897e-18, - "particle_momentum_z": 1.1156457989239732e-18, - "particle_position_x": 102495.14197173176, - "particle_position_y": 188132.22608016344, - "particle_position_z": 102423.13701045913, + "lev=0": {}, + "electron_negative": { + "particle_momentum_x": 1.1222699783863554e-18, + "particle_momentum_y": 1.1202176725070554e-18, + "particle_momentum_z": 1.3925955132362978e-18, + "particle_position_x": 102352.09026544492, + "particle_position_y": 102418.88243172191, + "particle_position_z": 194298.7949373403, "particle_weight": 8.959999999999998e-07 }, - "lev=0": {}, "proton": { - "particle_momentum_x": 3.835423016604918e-15, - "particle_momentum_y": 2.0468371931479925e-15, - "particle_momentum_z": 2.055186547721331e-15, - "particle_position_x": 189256.1546041931, - "particle_position_y": 102293.00576740496, - "particle_position_z": 102314.93877691089, + "particle_momentum_x": 3.8338884590187296e-15, + "particle_momentum_y": 2.0442156829943128e-15, + "particle_momentum_z": 2.045804260395492e-15, + "particle_position_x": 189238.69249885075, + "particle_position_y": 102242.91543133644, + "particle_position_z": 102297.92915049737, + "particle_weight": 8.959999999999998e-07 + }, + "electron": { + "particle_momentum_x": 1.1150196665556376e-18, + "particle_momentum_y": 2.2311586451156107e-18, + "particle_momentum_z": 1.1115069298757383e-18, + "particle_position_x": 102477.68142952351, + "particle_position_y": 188137.20906834095, + "particle_position_z": 102443.44417709563, + "particle_weight": 8.959999999999998e-07 + }, + "proton_negative": { + "particle_momentum_x": 2.051429985038282e-15, + "particle_momentum_y": 2.053711846305655e-15, + "particle_momentum_z": 2.7212135003240815e-15, + "particle_position_x": 102307.73649638034, + "particle_position_y": 102475.13878406698, + "particle_position_z": 193638.89296895845, "particle_weight": 8.959999999999998e-07 } } diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index b1a3ac7b0c8..d8409258fa0 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -104,37 +104,40 @@ namespace { // Momentum to be returned at the end of this function amrex::Real u = 0._rt; + const amrex::Real abs_u_m = std::abs(u_m); + if (u_th == 0._rt) { u = u_m; // Trivial case ; avoids division by 0 in the rest of the code below - } else if (u_m < 0.6*u_th) { - // Mean velocity is lower than thermal velocity - // Use the distribution u*exp(-u**2*(1-u_m/u_th)/(2*u_th**2)) as an approximation + } else if (abs_u_m < 0.6*u_th) { + // Mean velocity magnitude is less than thermal velocity + // Use the distribution u*exp(-u**2*(1-abs(u_m)/u_th)/(2*u_th**2)) as an approximation // and then use the rejection method to correct it - // ( stop rejecting with probability exp(-u_m/(2*u_th**3)*(u-u_th)**2) ) + // ( stop rejecting with probability exp(-abs(u_m)/(2*u_th**3)*(u-sign(u_m)*u_th)**2) ) // Note that this is the method that is used in the common case u_m=0 - const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - u_m/u_th ); - const amrex::Real reject_prefactor = (u_m/u_th)/(2._rt*u_th*u_th); // To save computation + const amrex::Real umsign = std::copysign(1._rt, u_m); + const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - abs_u_m/u_th ); + const amrex::Real reject_prefactor = (abs_u_m/u_th)/(2._rt*u_th*u_th); // To save computation bool reject = true; while (reject) { - // Generates u according to u*exp(-u**2/(2*approx_u_th**2), + // Generates u according to u*exp(-u**2/(2*approx_u_th**2)), // using the method of the inverse cumulative function amrex::Real xrand = 1._rt - amrex::Random(engine); // ensures urand > 0 u = approx_u_th * std::sqrt(2._rt*std::log(1._rt/xrand)); // Rejection method xrand = amrex::Random(engine); - if (xrand < std::exp(-reject_prefactor*(u-u_th)*(u-u_th))) reject = false; + if (xrand < std::exp(-reject_prefactor*(u - umsign*u_th)*(u - umsign*u_th))) reject = false; } } else { - // Mean velocity is greater than thermal velocity - // Use the distribution exp(-(u-u_m-u_th**2/u_m)**2/(2*u_th**2)) as an approximation + // Mean velocity magnitude is greater than thermal velocity + // Use the distribution exp(-(u-u_m-u_th**2/abs(u_m))**2/(2*u_th**2)) as an approximation // and then use the rejection method to correct it - // ( stop rejecting with probability (u/u_m)*exp(1-(u/u_m)) ; note + // ( stop rejecting with probability (u/abs(u_m))*exp(1-(u/abs(u_m))) ; note // that this number is always between 0 and 1 ) // Note that in the common case `u_m = 0`, this rejection method // is not used, and the above rejection method is used instead. bool reject = true; - const amrex::Real approx_u_m = u_m + u_th*u_th/u_m; - const amrex::Real inv_um = 1._rt/u_m; // To save computation + const amrex::Real approx_u_m = u_m + u_th*u_th/abs_u_m; + const amrex::Real inv_um = 1._rt/abs_u_m; // To save computation while (reject) { // Approximate distribution: normal distribution, where we only retain positive u u = -1._rt; @@ -166,13 +169,6 @@ struct InjectorMomentumGaussianFlux m_flux_normal_axis(a_flux_normal_axis), m_flux_direction(a_flux_direction) { - // For now, do not allow negative `u_m` along the flux axis - bool raise_error = false; - if ((m_flux_normal_axis == 0) && (m_ux_m < 0)) raise_error = true; - if ((m_flux_normal_axis == 1) && (m_uy_m < 0)) raise_error = true; - if ((m_flux_normal_axis == 2) && (m_uz_m < 0)) raise_error = true; - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( raise_error==false, - "When using the `gaussianflux` distribution, the central momentum along the flux axis must be positive or zero." ); } AMREX_GPU_HOST_DEVICE