diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 57242047a58..2c0428b0bdf 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1807,10 +1807,13 @@ Particle initialization * ``.maxwellian_u_std_distribution_type`` (`string`, default ``constant``): Specifies the distribution type for the thermal spread (standard deviation) of the - particle momentum. Here, ``u_std`` is a 3D vector (with components ``ux_std``, + particle momentum. + Here, ``u_std`` is a 3D vector (with components ``ux_std``, ``uy_std``, ``uz_std``) representing the standard deviation of the normalized momentum :math:`u_\mathrm{std} = \sqrt{\theta}`, where :math:`\theta = \frac{k_\mathrm{B} \cdot T}{m \cdot c^2}`. + Mutually exclusive with + ``.maxwellian_temperature_in_eV_distribution_type``. * If ``constant``, the following are required: ``.ux_std``, ``.uy_std``, ``.uz_std`` (`float`, default ``0``). @@ -1831,6 +1834,26 @@ Particle initialization Particles may be relativistic in the lab frame, but the sampling model treats them as non-relativistic in the drift frame. For a relativistic thermal spread, use ``maxwell_juttner`` instead. + * ``.maxwellian_temperature_in_eV_distribution_type`` (`string`): + Specifies the temperature in eV for the thermal spread of the particle momentum. This is an alternative to and mutually + exclusive with ``.maxwellian_u_std_distribution_type`` for specifying the thermal spread. + Under the hood, the standard deviation of each normalized momentum component in the drift frame is computed as: + :math:`u_\mathrm{std} = \sqrt{\mathrm{temperature\_in\_eV}\, q_e / (m c^2)}`, where + :math:`m` is the species mass (from ``species_type`` or ``mass``). + + * If ``constant``, the following is required: ``.temperature_in_eV`` (`float`). + * If ``parser``, the following is required: + ``.temperature_in_eV_function(x,y,z)``. + * If ``read_from_file``, ``temperature_in_eV`` is read as a scalar function of position + from an openPMD file and converted to an isotropic vector + :math:`u_\mathrm{std}` at the particle positions (requires a WarpX build with openPMD; + not supported yet in ``RZ`` / ``RCYLINDER`` / ``RSPHERE``). The following is required: + ``.read_temperature_in_eV_from_path`` (openPMD file path). The file must + contain a scalar openPMD mesh with the name given by + ``.temperature_in_eV_mesh_name`` (default ``temperature_in_eV``). See + `this file `__ + for an example of how to prepare the openPMD data file. + * ``maxwell_juttner``: Maxwell-Juttner distribution for relativistic plasma. More specifically, the plasma is initialized with a Maxwell-Juttner distribution @@ -1870,12 +1893,23 @@ Particle initialization `this file `__ for an example of how to prepare the openPMD data file. - * ``.theta_distribution_type`` (`string`, default ``constant``): - Specifies the distribution type for the temperature :math:`\theta`. - Values less than zero are not allowed. - - * If ``constant``, the following is required: ``.theta`` (`float`). - * If ``parser``, the following is required: ``.theta_function(x,y,z)``. + * ``.maxwell_juttner_temperature_in_eV_distribution_type`` (`string`, default ``constant``): + Specifies the distribution type for the temperature in eV. + Internally converted to dimensionless :math:`\theta = \mathrm{temperature\_in\_eV}\, q_e / (m c^2)`, + where :math:`m` is the species mass (from ``species_type`` or ``mass``). + Values of ``temperature_in_eV`` less than zero are not allowed. + + * If ``constant``, the following is required: ``.temperature_in_eV`` (`float`). + * If ``parser``, the following is required: + ``.temperature_in_eV_function(x,y,z)``. + * If ``read_from_file``, ``temperature_in_eV`` is read as a scalar function of position + from an openPMD file (requires a WarpX build with openPMD; not supported yet in ``RZ`` / ``RCYLINDER`` / + ``RSPHERE``). The following is required: + ``.read_temperature_in_eV_from_path`` (openPMD file path). The file must + contain a scalar openPMD mesh with the name given by + ``.temperature_in_eV_mesh_name`` (default ``temperature_in_eV``). See + `this file `__ + for an example of how to prepare the openPMD data file. Sampling uses the Sobol and flipping methods described in :cite:t:`param-ZenitaniPOP2015`. For :math:`\theta \lesssim 0.1`, the Sobol method becomes inefficient (its acceptance diff --git a/Examples/Tests/initial_distribution/analysis.py b/Examples/Tests/initial_distribution/analysis.py index 6feacf58aac..cdc15e819ec 100755 --- a/Examples/Tests/initial_distribution/analysis.py +++ b/Examples/Tests/initial_distribution/analysis.py @@ -190,56 +190,64 @@ assert charge_error < tolerance # ============================================= -# maxwell-juttner with temperature from parser +# maxwell-juttner with temperature from parser (h5_neg/h5_pos folders) or read from file (h14_neg/h14_pos folders) # ============================================= -# load data -bin_value, bin_data_neg = read_reduced_diags_histogram("h5_neg.txt")[2:] -bin_data_pos = read_reduced_diags_histogram("h5_pos.txt")[3] +for i in [5, 14]: + # load data + bin_value, bin_data_neg = read_reduced_diags_histogram(f"h{i}_neg.txt")[2:] + + bin_data_pos = read_reduced_diags_histogram(f"h{i}_pos.txt")[3] + + # parameters of theory + # _neg denotes where x < 0, _pos where x > 0 + theta_neg = 1.0 + theta_pos = 2.0 + + K2_neg = scs.kn(2, 1.0 / theta_neg) + K2_pos = scs.kn(2, 1.0 / theta_pos) + + n = 1.0e21 + V = 8.0 / 2 # each histogram is for half the domain + db = 0.22 + + # compute the analytical solution for each half of the domain + f_neg = ( + n + * V + * db + * bin_value**2 + * np.sqrt(1.0 - 1.0 / bin_value**2) + / (theta_neg * K2_neg) + * np.exp(-bin_value / theta_neg) + ) -# parameters of theory -# _neg denotes where x<0, _pos where x>0 -theta_neg = 1.0 -theta_pos = 2.0 -K2_neg = scs.kn(2, 1.0 / theta_neg) -K2_pos = scs.kn(2, 1.0 / theta_pos) -n = 1.0e21 -V = 8.0 / 2 # because each of these are for half the domain -db = 0.22 + f_pos = ( + n + * V + * db + * bin_value**2 + * np.sqrt(1.0 - 1.0 / bin_value**2) + / (theta_pos * K2_pos) + * np.exp(-bin_value / theta_pos) + ) -# compute the analytical solution for each half of the domain -f_neg = ( - n - * V - * db - * bin_value**2 - * np.sqrt(1.0 - 1.0 / bin_value**2) - / (theta_neg * K2_neg) - * np.exp(-bin_value / theta_neg) -) -f_neg_peak = np.amax(f_neg) -f_pos = ( - n - * V - * db - * bin_value**2 - * np.sqrt(1.0 - 1.0 / bin_value**2) - / (theta_pos * K2_pos) - * np.exp(-bin_value / theta_pos) -) -f_pos_peak = np.amax(f_pos) -f_peak = max(f_neg_peak, f_pos_peak) + # common normalization + f_peak = max(np.amax(f_neg), np.amax(f_pos)) -# compute error -f5_error = ( - np.sum(np.abs(f_neg - bin_data_neg) + np.abs(f_pos - bin_data_pos)) - / bin_value.size - / f_peak -) + # compute error + error = ( + np.sum(np.abs(f_neg - bin_data_neg) + np.abs(f_pos - bin_data_pos)) + / bin_value.size + / f_peak + ) -print("Maxwell-Juttner parser temperature difference:", f5_error) + if i == 5: + print("Maxwell-Juttner parser temperature difference:", error) + elif i == 14: + print("Maxwell-Juttner read from file temperature difference:", error) + assert error < tolerance -assert f5_error < tolerance # ================================================= # maxwell-juttner with a constant asymmetric bulk drift @@ -313,43 +321,44 @@ assert f10_error < tolerance # ============================================== -# maxwellian with constant bulk velocity +# maxwellian with constant bulk (h6/h6uy folders) velocity and constant temperature in eV (h12/h12uy folders) # ============================================== -# load data -bin_value_g, bin_data_g = read_reduced_diags_histogram("h6.txt")[2:] -bin_value_uy, bin_data_uy = read_reduced_diags_histogram("h6uy.txt")[2:] - -# Expected values for beta and u = beta*gamma -beta_const = 0.2 -g_const = 1.0 / np.sqrt(1.0 - beta_const * beta_const) -uy_const = beta_const * g_const -g_bin_size = 0.004 -g_bin_min = 1.0 -uy_bin_size = 0.04 -uy_bin_min = -1.0 -V = 8.0 # volume in m^3 -n = 1.0e21 # number density in 1/m^3 - -f_g = np.zeros_like(bin_value_g) -i_g = int(np.floor((g_const - g_bin_min) / g_bin_size)) -f_g[i_g] = n * V -f_peak = np.amax(f_g) - -f_uy = np.zeros_like(bin_value_uy) -i_uy = int(np.floor((-uy_const - uy_bin_min) / uy_bin_size)) -f_uy[i_uy] = n * V - -f6_error = ( - np.sum(np.abs(f_g - bin_data_g) + np.abs(f_uy - bin_data_uy)) - / bin_value_g.size - / f_peak -) +for i in [6, 12]: + # load data + bin_value_g, bin_data_g = read_reduced_diags_histogram(f"h{i}.txt")[2:] + bin_value_uy, bin_data_uy = read_reduced_diags_histogram(f"h{i}uy.txt")[2:] + + # Expected values for beta and u = beta*gamma + beta_const = 0.2 + g_const = 1.0 / np.sqrt(1.0 - beta_const * beta_const) + uy_const = beta_const * g_const + g_bin_size = 0.004 + g_bin_min = 1.0 + uy_bin_size = 0.04 + uy_bin_min = -1.0 + V = 8.0 # volume in m^3 + n = 1.0e21 # number density in 1/m^3 -print("Maxwell-Boltzmann constant velocity difference:", f6_error) + f_g = np.zeros_like(bin_value_g) + i_g = int(np.floor((g_const - g_bin_min) / g_bin_size)) + f_g[i_g] = n * V + f_peak = np.amax(f_g) -assert f6_error < tolerance + f_uy = np.zeros_like(bin_value_uy) + i_uy = int(np.floor((-uy_const - uy_bin_min) / uy_bin_size)) + f_uy[i_uy] = n * V + error = ( + np.sum(np.abs(f_g - bin_data_g) + np.abs(f_uy - bin_data_uy)) + / bin_value_g.size + / f_peak + ) + if i == 6: + print("Maxwell-Boltzmann constant velocity difference:", error) + elif i == 12: + print("Maxwell-Boltzmann constant temperature_in_eV difference:", error) + assert error < tolerance # ============================================ # maxwellian with parser bulk velocity # ============================================ @@ -405,7 +414,9 @@ def check_standard_normal(u, mean_ref, std_ref, tolerance): r = (u - mean_ref) / std_ref r_mean = np.mean(r) r_std = np.std(r) + print(" abs(r_mean) = ", abs(r_mean)) assert abs(r_mean) < tolerance + print(" abs(r_std - 1.0) = ", abs(r_std - 1.0)) assert abs(r_std - 1.0) < tolerance @@ -434,6 +445,27 @@ def check_standard_normal(u, mean_ref, std_ref, tolerance): check_standard_normal(uz, uz_mean_interp, uz_std_interp, standard_normal_tolerance) +# ============================================== +# maxwellian with bulk velocity and temperature_in_eV from openPMD file +# (isotropic u_std = 0.2 * |z| from temperature_in_eV) +# ============================================== +standard_normal_tolerance = 8e-2 + +ux, uy, uz, z = ts.get_particle( + ["ux", "uy", "uz", "z"], + species="gaussian_temperature_in_eV_from_file", + iteration=0, +) + +ux_std_interp = np.interp(z, z_array, 0.2 * np.abs(z_array)) +uy_std_interp = np.interp(z, z_array, 0.2 * np.abs(z_array)) +uz_std_interp = np.interp(z, z_array, 0.2 * np.abs(z_array)) + +check_standard_normal(ux, ux_mean_interp, ux_std_interp, standard_normal_tolerance) +check_standard_normal(uy, uy_mean_interp, uy_std_interp, standard_normal_tolerance) +check_standard_normal(uz, uz_mean_interp, uz_std_interp, standard_normal_tolerance) + + # ============================================ # Cuboid distribution in momentum space # ============================================ @@ -528,14 +560,9 @@ def check_validity_uniform(bins, histogram, u_min, u_max, Ntrials=1000): check_validity_uniform(bin_value_z, h8z[timestep] / N0, uz_min, uz_max) # ================================================= -# Gaussian with parser mean and standard deviation +# Gaussian with parser mean and standard deviation (h9x,h9y,h9z folders) and parser temperature in eV (h13x/h13y/h13z folders) # ================================================= -# load data -bin_value_ux, bin_data_ux = read_reduced_diags_histogram("h9x.txt")[2:] -bin_value_uy, bin_data_uy = read_reduced_diags_histogram("h9y.txt")[2:] -bin_value_uz, bin_data_uz = read_reduced_diags_histogram("h9z.txt")[2:] - def Gaussian(mean, sigma, u): V = 8.0 # volume in m^3 @@ -545,20 +572,33 @@ def Gaussian(mean, sigma, u): ) -du = 2.0 / 50 -f_ux = Gaussian(0.1, 0.2, bin_value_ux) * du -f_uy = Gaussian(0.12, 0.21, bin_value_uy) * du -f_uz = Gaussian(0.14, 0.22, bin_value_uz) * du - -f9_error = ( - np.sum( - np.abs(f_ux - bin_data_ux) / f_ux.max() - + np.abs(f_uy - bin_data_uy) / f_ux.max() - + np.abs(f_uz - bin_data_uz) / f_uz.max() +for i, i_tolerance in [(9, tolerance), (13, tolerance)]: + # load data + bin_value_ux, bin_data_ux = read_reduced_diags_histogram(f"h{i}x.txt")[2:] + bin_value_uy, bin_data_uy = read_reduced_diags_histogram(f"h{i}y.txt")[2:] + bin_value_uz, bin_data_uz = read_reduced_diags_histogram(f"h{i}z.txt")[2:] + + du = 2.0 / 50 + + if i == 9: + f_ux = Gaussian(0.1, 0.2, bin_value_ux) * du + f_uy = Gaussian(0.12, 0.21, bin_value_uy) * du + f_uz = Gaussian(0.14, 0.22, bin_value_uz) * du + elif i == 13: + f_ux = Gaussian(0.1, 0.2, bin_value_ux) * du + f_uy = Gaussian(0.12, 0.2, bin_value_uy) * du + f_uz = Gaussian(0.14, 0.2, bin_value_uz) * du + + error = ( + np.sum( + np.abs(f_ux - bin_data_ux) / f_ux.max() + + np.abs(f_uy - bin_data_uy) / f_ux.max() + + np.abs(f_uz - bin_data_uz) / f_uz.max() + ) + / bin_value_ux.size ) - / bin_value_ux.size -) - -print("maxwellian parser mean/std velocity difference:", f9_error) - -assert f9_error < tolerance + if i == 9: + print("Maxwellian parser mean/std velocity difference:", error) + elif i == 13: + print("Maxwellian parser mean/temperature_in_eV difference:", error) + assert error < i_tolerance diff --git a/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution b/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution index c866110ab2a..4c2dd43548b 100644 --- a/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution +++ b/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution @@ -29,7 +29,7 @@ algo.particle_shape = 1 ################################# ############ PLASMA ############# ################################# -particles.species_names = gaussian maxwell_boltzmann maxwell_juttner maxwell_juttner_drift maxwell_juttner_low_theta beam maxwell_juttner_parser velocity_constant velocity_parser uniform gaussian_parser gaussian_momentum_from_file +particles.species_names = gaussian maxwell_boltzmann maxwell_juttner maxwell_juttner_drift maxwell_juttner_low_theta beam maxwell_juttner_parser velocity_constant velocity_parser uniform gaussian_parser gaussian_momentum_from_file temperature_in_eV_constant gaussian_temperature_in_eV_from_file gaussian_parser_temperature_in_eV_from_file maxwell_juttner_read_from_file particles.rigid_injected_species = beam gaussian.charge = -q_e @@ -66,7 +66,8 @@ maxwell_juttner.num_particles_per_cell = 1000 maxwell_juttner.profile = constant maxwell_juttner.density = 1.0e21 maxwell_juttner.momentum_distribution_type = "maxwell_juttner" -maxwell_juttner.theta = 1.0 +maxwell_juttner.maxwell_juttner_temperature_in_eV_distribution_type = "constant" +maxwell_juttner.temperature_in_eV = m_e * clight**2 / q_e # Maxwell-Juttner with an asymmetric bulk drift (normalized momentum gamma*v/c). # ux_mean=0.5, uy_mean=0.3, uz_mean=0.1 gives |u_mean|^2=0.35, gamma_bulk=sqrt(1.35). @@ -79,7 +80,8 @@ maxwell_juttner_drift.num_particles_per_cell = 1000 maxwell_juttner_drift.profile = constant maxwell_juttner_drift.density = 1.0e21 maxwell_juttner_drift.momentum_distribution_type = "maxwell_juttner" -maxwell_juttner_drift.theta = 1.0 +maxwell_juttner_drift.maxwell_juttner_temperature_in_eV_distribution_type = "constant" +maxwell_juttner_drift.temperature_in_eV = m_e * clight**2 / q_e maxwell_juttner_drift.maxwell_juttner_u_mean_distribution_type = "constant" maxwell_juttner_drift.ux_mean = 0.5 maxwell_juttner_drift.uy_mean = 0.3 @@ -92,7 +94,8 @@ maxwell_juttner_low_theta.num_particles_per_cell = 1000 maxwell_juttner_low_theta.profile = constant maxwell_juttner_low_theta.density = 1.0e21 maxwell_juttner_low_theta.momentum_distribution_type = "maxwell_juttner" -maxwell_juttner_low_theta.theta = 0.05 +maxwell_juttner_low_theta.maxwell_juttner_temperature_in_eV_distribution_type = "constant" +maxwell_juttner_low_theta.temperature_in_eV = 0.05 * m_e * clight**2 / q_e beam.charge = -q_e beam.mass = m_e @@ -123,8 +126,18 @@ maxwell_juttner_parser.num_particles_per_cell = 1000 maxwell_juttner_parser.profile = constant maxwell_juttner_parser.density = 1.0e21 maxwell_juttner_parser.momentum_distribution_type = "maxwell_juttner" -maxwell_juttner_parser.theta_distribution_type = "parser" -maxwell_juttner_parser.theta_function(x,y,z) = "1.0 + heaviside(x,0)" +maxwell_juttner_parser.maxwell_juttner_temperature_in_eV_distribution_type = "parser" +maxwell_juttner_parser.temperature_in_eV_function(x,y,z) = "(1.0 + heaviside(x,0)) * m_e * clight**2 / q_e" + +maxwell_juttner_read_from_file.charge = -q_e +maxwell_juttner_read_from_file.mass = m_e +maxwell_juttner_read_from_file.injection_style = "NRandomPerCell" +maxwell_juttner_read_from_file.num_particles_per_cell = 1000 +maxwell_juttner_read_from_file.profile = constant +maxwell_juttner_read_from_file.density = 1.0e21 +maxwell_juttner_read_from_file.momentum_distribution_type = "maxwell_juttner" +maxwell_juttner_read_from_file.maxwell_juttner_temperature_in_eV_distribution_type = "read_from_file" +maxwell_juttner_read_from_file.read_temperature_in_eV_from_path = "../test_3d_initial_distribution_prepare/example-temperature-in-eV-MJ.h5" velocity_constant.charge = -q_e velocity_constant.mass = m_e @@ -142,6 +155,20 @@ velocity_constant.ux_mean =0 velocity_constant.uy_mean = -0.20412414523193154 velocity_constant.uz_mean =0 +temperature_in_eV_constant.charge = -q_e +temperature_in_eV_constant.mass = m_e +temperature_in_eV_constant.injection_style = "NRandomPerCell" +temperature_in_eV_constant.num_particles_per_cell = 1000 +temperature_in_eV_constant.profile = constant +temperature_in_eV_constant.density = 1.0e21 +temperature_in_eV_constant.momentum_distribution_type = "maxwellian" +temperature_in_eV_constant.maxwellian_temperature_in_eV_distribution_type= "constant" +temperature_in_eV_constant.temperature_in_eV = 0.0005109989506917532 +temperature_in_eV_constant.maxwellian_u_mean_distribution_type = "constant" +temperature_in_eV_constant.ux_mean =0 +temperature_in_eV_constant.uy_mean = -0.20412414523193154 +temperature_in_eV_constant.uz_mean =0 + velocity_parser.charge = -q_e velocity_parser.mass = m_e velocity_parser.injection_style = "NRandomPerCell" @@ -200,6 +227,32 @@ gaussian_momentum_from_file.read_u_mean_from_path = "../test_3d_initial_distribu gaussian_momentum_from_file.maxwellian_u_std_distribution_type = "read_from_file" gaussian_momentum_from_file.read_u_std_from_path = "../test_3d_initial_distribution_prepare/example-u-std.h5" +gaussian_temperature_in_eV_from_file.charge = -q_e +gaussian_temperature_in_eV_from_file.mass = m_e +gaussian_temperature_in_eV_from_file.injection_style = "NRandomPerCell" +gaussian_temperature_in_eV_from_file.num_particles_per_cell = 1000 +gaussian_temperature_in_eV_from_file.profile = constant +gaussian_temperature_in_eV_from_file.density = 1.0e21 +gaussian_temperature_in_eV_from_file.momentum_distribution_type = "maxwellian" +gaussian_temperature_in_eV_from_file.maxwellian_u_mean_distribution_type = "read_from_file" +gaussian_temperature_in_eV_from_file.read_u_mean_from_path = "../test_3d_initial_distribution_prepare/example-u-mean.h5" +gaussian_temperature_in_eV_from_file.maxwellian_temperature_in_eV_distribution_type = "read_from_file" +gaussian_temperature_in_eV_from_file.read_temperature_in_eV_from_path = "../test_3d_initial_distribution_prepare/example-temperature-in-eV.h5" + +gaussian_parser_temperature_in_eV_from_file.charge = -q_e +gaussian_parser_temperature_in_eV_from_file.mass = m_e +gaussian_parser_temperature_in_eV_from_file.injection_style = "NRandomPerCell" +gaussian_parser_temperature_in_eV_from_file.num_particles_per_cell = 1000 +gaussian_parser_temperature_in_eV_from_file.profile = constant +gaussian_parser_temperature_in_eV_from_file.density = 1.0e21 +gaussian_parser_temperature_in_eV_from_file.momentum_distribution_type = "maxwellian" +gaussian_parser_temperature_in_eV_from_file.maxwellian_u_mean_distribution_type= "parser" +gaussian_parser_temperature_in_eV_from_file.ux_mean_function(x,y,z) = "0.1*z" +gaussian_parser_temperature_in_eV_from_file.uy_mean_function(x,y,z) = "0.12*z" +gaussian_parser_temperature_in_eV_from_file.uz_mean_function(x,y,z) = "0.14*z" +gaussian_parser_temperature_in_eV_from_file.maxwellian_temperature_in_eV_distribution_type= "parser" +gaussian_parser_temperature_in_eV_from_file.temperature_in_eV_function(x,y,z) = "20439.958027670127 * z**2" + ################################# ########## DIAGNOSTIC ########### ################################# @@ -214,7 +267,10 @@ gaussian_momentum_from_file.read_u_std_from_path = "../test_3d_initial_distribut # 10 for maxwell-juttner with low theta (Gaussian fallback) # 11 for maxwell-juttner with constant diagonal bulk drift # 12 for maxwellian with mean and standard deviation from openPMD files -warpx.reduced_diags_names = h1x h1y h1z h2x h2y h2z h3 h3_filtered h4x h4y h4z bmmntr h5_neg h5_pos h6 h6uy h7 h7uy_pos h7uy_neg h8x h8y h8z h9x h9y h9z h10 h11 +# 13 for maxwellian with parser mean and parser temperature in eV +# 14 for maxwell-juttner with temperature in eV read from file + +warpx.reduced_diags_names = h1x h1y h1z h2x h2y h2z h3 h3_filtered h4x h4y h4z bmmntr h5_neg h5_pos h6 h6uy h7 h7uy_pos h7uy_neg h8x h8y h8z h9x h9y h9z h10 h11 h12 h12uy h13x h13y h13z h14_neg h14_pos h1x.type = ParticleHistogram h1x.intervals = 1 @@ -336,16 +392,6 @@ h5_pos.bin_max = 12.0 h5_pos.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)" h5_pos.filter_function(t,x,y,z,ux,uy,uz) = "x > 0.0" -h5_pos.type = ParticleHistogram -h5_pos.intervals = 1 -h5_pos.path = "./" -h5_pos.species = maxwell_juttner_parser -h5_pos.bin_number = 50 -h5_pos.bin_min = 1.0 -h5_pos.bin_max = 12.0 -h5_pos.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)" -h5_pos.filter_function(t,x,y,z,ux,uy,uz) = "x > 0.0" - h10.type = ParticleHistogram h10.intervals = 1 h10.path = "./" @@ -468,6 +514,71 @@ h11.bin_min = 1.0 h11.bin_max = 12.0 h11.histogram_function(t,x,y,z,ux,uy,uz) = "1.161895003862225*sqrt(1.0 + ux*ux + uy*uy + uz*uz) - 0.5*ux - 0.3*uy - 0.1*uz" +h12.type = ParticleHistogram +h12.intervals = 1 +h12.path = "./" +h12.species = temperature_in_eV_constant +h12.bin_number = 50 +h12.bin_min = 1.0 +h12.bin_max = 1.2 +h12.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)" + +h12uy.type = ParticleHistogram +h12uy.intervals = 1 +h12uy.path = "./" +h12uy.species = temperature_in_eV_constant +h12uy.bin_number = 50 +h12uy.bin_min = -1.0 +h12uy.bin_max = 1.0 +h12uy.histogram_function(t,x,y,z,ux,uy,uz) = "uy" + +h13x.type = ParticleHistogram +h13x.intervals = 1 +h13x.path = "./" +h13x.species = gaussian_parser_temperature_in_eV_from_file +h13x.bin_number = 50 +h13x.bin_min = -1 +h13x.bin_max = 1 +h13x.histogram_function(t,x,y,z,ux,uy,uz) = "ux/z" + +h13y.type = ParticleHistogram +h13y.intervals = 1 +h13y.path = "./" +h13y.species = gaussian_parser_temperature_in_eV_from_file +h13y.bin_number = 50 +h13y.bin_min = -1 +h13y.bin_max = 1 +h13y.histogram_function(t,x,y,z,ux,uy,uz) = "uy/z" + +h13z.type = ParticleHistogram +h13z.intervals = 1 +h13z.path = "./" +h13z.species = gaussian_parser_temperature_in_eV_from_file +h13z.bin_number = 50 +h13z.bin_min = -1 +h13z.bin_max = 1 +h13z.histogram_function(t,x,y,z,ux,uy,uz) = "uz/z" + +h14_neg.type = ParticleHistogram +h14_neg.intervals = 1 +h14_neg.path = "./" +h14_neg.species = maxwell_juttner_read_from_file +h14_neg.bin_number = 50 +h14_neg.bin_min = 1.0 +h14_neg.bin_max = 12.0 +h14_neg.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)" +h14_neg.filter_function(t,x,y,z,ux,uy,uz) = "x < 0.0" + +h14_pos.type = ParticleHistogram +h14_pos.intervals = 1 +h14_pos.path = "./" +h14_pos.species = maxwell_juttner_read_from_file +h14_pos.bin_number = 50 +h14_pos.bin_min = 1.0 +h14_pos.bin_max = 12.0 +h14_pos.histogram_function(t,x,y,z,ux,uy,uz) = "sqrt(1.0 + ux*ux + uy*uy + uz*uz)" +h14_pos.filter_function(t,x,y,z,ux,uy,uz) = "x > 0.0" + # our little beam monitor bmmntr.type = BeamRelevant bmmntr.intervals = 1 diff --git a/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py b/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py index 9299ea9a0d1..af2a667c9b8 100644 --- a/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py +++ b/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py @@ -8,6 +8,9 @@ import numpy as np import openpmd_api as io +from scipy.constants import c as clight +from scipy.constants import electron_mass as m_e +from scipy.constants import elementary_charge as q_e # Define u_mean and u_std as functions of x, y, z, using numpy syntax # - Define the grid @@ -16,7 +19,7 @@ z_1d = np.linspace(-1.0, 1.0, 8) x, y, z = np.meshgrid(x_1d, y_1d, z_1d, indexing="ij") -# - Define the normalized momentum data, u = gamma * v / c +# Define the normalized momentum data, u = gamma * v / c ux_std_data = 0.2 * abs(z) uy_std_data = 0.21 * abs(z) uz_std_data = 0.22 * abs(z) @@ -25,6 +28,12 @@ uy_mean_data = 0.12 * z uz_mean_data = 0.14 * z +# Define temperature in eV to correspond above isotropic thermal spread u_std = 0.2 * |z| +# temperature [eV] = u_std^2 * m_e * c^2 / q_e +temperature_in_eV_data = 20439.95 * z**2 + +temperature_in_eV_data_MJ = (1.0 + np.heaviside(x, 0)) * m_e * clight**2 / q_e + grid_spacing = np.array( [ x_1d[1] - x_1d[0], @@ -58,6 +67,25 @@ def write_vector_mesh_file(filename, mesh_name, components): series.flush() +def write_scalar_mesh_file(filename, mesh_name, data): + series = io.Series(filename, io.Access.create) + it = series.iterations[1] + + mesh = it.meshes[mesh_name] + mesh.grid_spacing = grid_spacing + mesh.grid_global_offset = grid_offset + mesh.axis_labels = ["x", "y", "z"] + mesh.geometry = io.Geometry.cartesian + mesh.unit_dimension = {} + + component = mesh[io.Mesh_Record_Component.SCALAR] + component.position = [0.0, 0.0, 0.0] + component.reset_dataset(io.Dataset(data.dtype, data.shape)) + component.store_chunk(data) + + series.flush() + + write_vector_mesh_file( "example-u-std.h5", "u_std", @@ -77,3 +105,15 @@ def write_vector_mesh_file(filename, mesh_name, components): "z": uz_mean_data, }, ) + +write_scalar_mesh_file( + "example-temperature-in-eV.h5", + "temperature_in_eV", + temperature_in_eV_data, +) + +write_scalar_mesh_file( + "example-temperature-in-eV-MJ.h5", + "temperature_in_eV", + temperature_in_eV_data_MJ, +) diff --git a/Regression/Checksum/benchmarks_json/test_3d_initial_distribution.json b/Regression/Checksum/benchmarks_json/test_3d_initial_distribution.json index 9412de9b622..4281f297934 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_initial_distribution.json +++ b/Regression/Checksum/benchmarks_json/test_3d_initial_distribution.json @@ -117,5 +117,41 @@ "particle_position_y": 256011.49684468308, "particle_position_z": 256109.9226490142, "particle_weight": 8e+21 + }, + "temperature_in_eV_constant": { + "particle_momentum_x": 3.5285787492061906e-21, + "particle_momentum_y": 2.854131991544734e-17, + "particle_momentum_z": 3.52328932816798e-21, + "particle_position_x": 256013.9399284002, + "particle_position_y": 256011.5933546236, + "particle_position_z": 255905.84694351768, + "particle_weight": 8e+21 + }, + "gaussian_temperature_in_eV_from_file": { + "particle_momentum_x": 1.2956059014974661e-17, + "particle_momentum_y": 1.3566264317285339e-17, + "particle_momentum_z": 1.4185798889934373e-17, + "particle_position_x": 255862.65084565995, + "particle_position_y": 255935.2699936382, + "particle_position_z": 255968.0187925212, + "particle_weight": 8e+21 + }, + "gaussian_parser_temperature_in_eV_from_file": { + "particle_momentum_x": 1.2537256955983462e-17, + "particle_momentum_y": 1.3101027839551143e-17, + "particle_momentum_z": 1.3760320087390352e-17, + "particle_position_x": 256024.46311264794, + "particle_position_y": 256018.00713414286, + "particle_position_z": 256001.97548263433, + "particle_weight": 8e+21 + }, + "maxwell_juttner_read_from_file": { + "particle_momentum_x": 3.240820500502024e-16, + "particle_momentum_y": 3.244055276686748e-16, + "particle_momentum_z": 3.2461559322585413e-16, + "particle_position_x": 255954.94884787325, + "particle_position_y": 256074.99143837203, + "particle_position_z": 255916.95081755883, + "particle_weight": 8e+21 } } diff --git a/Source/Initialization/GetTemperature.H b/Source/Initialization/GetTemperature.H index 452ada94472..0baa1a02408 100644 --- a/Source/Initialization/GetTemperature.H +++ b/Source/Initialization/GetTemperature.H @@ -14,6 +14,21 @@ #include +#include + +/** Convert a scalar temperature in eV to an isotropic u_std vector. */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::XDim3 +ConvertTemperatureInEvToUStd (amrex::Real const T_eV, + amrex::Real const q_e_over_mc2) noexcept +{ + if (T_eV < 0.0) { + amrex::Abort("The Maxwellian temperature in eV must be non-negative."); + } + const amrex::Real u_std = std::sqrt(T_eV * q_e_over_mc2); + return amrex::XDim3{u_std, u_std, u_std}; +} + /** * \brief Get temperature at a point on the grid * @@ -21,16 +36,27 @@ * or a spatially varying value computed using the parser function (m_temperature_parser). * It provides the temperature information held by the TemperatureProperties instance * passed to the constructor. + * + * Returns temperature_in_eV, which is used to compute: + * - the thermal spread in InjectorMomentumJuttner for the Maxwell-Juttner + * momentum distribution (e.g., temperature_in_eV * q_e_over_mc2); + * - the isotropic u_std vector for the Maxwellian momentum distribution (e.g., sqrt(temperature_in_eV * q_e_over_mc2)). */ struct GetTemperature { /* Type of temperature initialization */ TemperatureInitType m_type; - /* Constant temperature value, if m_type == TempConstantValue */ + /* Constant temperature in eV, if m_type == TempConstantValue */ amrex::Real m_temperature; - /* Temperature parser function, if m_type == TempParserFunction */ + /* Parser temperature in eV, if m_type == TempParserFunction */ amrex::ParserExecutor<3> m_temperature_parser; + /* q_e / (m c^2) factor, used for temperature_in_eV conversion */ + amrex::Real m_q_e_over_mc2{0.0}; +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ + !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) + ExternalFieldView m_temperature_in_eV_from_file; +#endif /** * \brief Construct the functor with information provided by temp @@ -41,15 +67,8 @@ struct GetTemperature explicit GetTemperature (TemperatureProperties const& temp) noexcept; /** - * \brief Functor call. Returns the value of temperature at the location (x,y,z) + * \brief Returns temperature_in_eV at (x,y,z). * - * \param[in] x: x-coordinate of given location - * \param[in] y: y-coordinate of given location - * \param[in] z: z-cooridnate of given location - * - *\return value of temperature at (x,y,z). - * m_temperature if m_type is TempConstantValue - * m_temperature_parser(x,y,z) if m_type is TempParserFunction */ AMREX_GPU_HOST_DEVICE amrex::Real operator() (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept @@ -64,6 +83,30 @@ struct GetTemperature { return m_temperature_parser(x,y,z); } +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ + !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) + case (TempFromFileValue): + { +#if (AMREX_SPACEDIM < 3) + amrex::ignore_unused(x,y,z); +#endif +#if (AMREX_SPACEDIM == 1) + amrex::RealVect const pos{z}; +#elif defined(WARPX_DIM_XZ) + amrex::RealVect const pos{x,z}; +#else + amrex::RealVect const pos{x,y,z}; +#endif +#if defined(AMREX_USE_GPU) + AMREX_IF_ON_HOST((amrex::Abort( + "For GPU builds, TempFromFileValue only works on device.");)) + AMREX_IF_ON_DEVICE((return m_temperature_in_eV_from_file(pos);)) + return 0.0; +#else + return m_temperature_in_eV_from_file(pos); +#endif + } +#endif default: { amrex::Abort("Get initial temperature: unknown type"); @@ -79,9 +122,12 @@ struct GetTemperatureVector amrex::Real m_ux_std{0}, m_uy_std{0}, m_uz_std{0}; amrex::ParserExecutor<3> m_ux_std_parser, m_uy_std_parser, m_uz_std_parser; + GetTemperature m_temperature; + amrex::Real m_q_e_over_mc2{0}; // q_e / (m c^2) #if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) - ExternalFieldVectorView m_from_file; + ExternalFieldVectorView m_from_file; // Vector openPMD u_std field, if m_type == TempFromFileVector + ExternalFieldView m_temperature_in_eV_from_file; // scalar openPMD temperature_in_eV field, if m_type == TempFromFileValue #endif /** * \brief Construct the functor with information provided by temp @@ -92,15 +138,15 @@ struct GetTemperatureVector explicit GetTemperatureVector (TemperatureProperties const& temp) noexcept; /** - * \brief Functor call. Returns the value of temperature at the location (x,y,z) + * \brief Returns the thermal spread as normalized momentum standard deviations + * (ux_std, uy_std, uz_std) in the drift frame at (x,y,z). * * \param[in] x: x-coordinate of given location * \param[in] y: y-coordinate of given location * \param[in] z: z-cooridnate of given location * - *\return value of temperature at (x,y,z). - * m_ux_std, m_uy_std, m_uz_std if m_type is TempConstantVector - * m_ux_std_parser(x,y,z), m_uy_std_parser(x,y,z), m_uz_std_parser(x,y,z) if m_type is TempParserFunctionVector + * \return amrex::XDim3{ux_std, uy_std, uz_std} — isotropic or anisotropic + * depending on m_type (direct u_std, or converted from temperature_in_eV). */ AMREX_GPU_HOST_DEVICE amrex::XDim3 operator() (amrex::Real const x, amrex::Real const y, amrex::Real const z) const noexcept @@ -115,12 +161,41 @@ struct GetTemperatureVector { return amrex::XDim3{m_ux_std_parser(x,y,z), m_uy_std_parser(x,y,z), m_uz_std_parser(x,y,z)}; } + case (TempConstantValue): + case (TempParserFunction): + { + return ConvertTemperatureInEvToUStd( + m_temperature(x,y,z), m_q_e_over_mc2); + } #if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) case (TempFromFileVector): { return m_from_file.getValue(x, y, z); } + case (TempFromFileValue): + { +#if (AMREX_SPACEDIM < 3) + amrex::ignore_unused(x,y,z); +#endif +#if (AMREX_SPACEDIM == 1) + amrex::RealVect const pos{z}; +#elif defined(WARPX_DIM_XZ) + amrex::RealVect const pos{x,z}; +#else + amrex::RealVect const pos{x,y,z}; +#endif +#if defined(AMREX_USE_GPU) + AMREX_IF_ON_HOST((amrex::Abort( + "For GPU builds, TempFromFileValue only works on device.");)) + AMREX_IF_ON_DEVICE((return ConvertTemperatureInEvToUStd( + m_temperature_in_eV_from_file(pos), m_q_e_over_mc2);)) + return amrex::XDim3{0.0, 0.0, 0.0}; +#else + return ConvertTemperatureInEvToUStd( + m_temperature_in_eV_from_file(pos), m_q_e_over_mc2); +#endif + } #endif default: { diff --git a/Source/Initialization/GetTemperature.cpp b/Source/Initialization/GetTemperature.cpp index c7699f73027..0ad2841aa1a 100644 --- a/Source/Initialization/GetTemperature.cpp +++ b/Source/Initialization/GetTemperature.cpp @@ -11,6 +11,7 @@ // Constructor for single-component (scalar) temperature GetTemperature::GetTemperature (TemperatureProperties const& temp) noexcept : m_type{temp.m_type} + , m_q_e_over_mc2{temp.m_q_e_over_mc2} { if (m_type == TempConstantValue) { m_temperature = temp.m_temperature; @@ -18,11 +19,19 @@ GetTemperature::GetTemperature (TemperatureProperties const& temp) noexcept else if (m_type == TempParserFunction) { m_temperature_parser = temp.m_ptr_temperature_parser->compile<3>(); } +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ + !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) + else if (m_type == TempFromFileValue) { + m_temperature_in_eV_from_file = temp.m_temperature_in_eV_reader->getView(); + } +#endif } // Constructor for three-component (vector) temperature GetTemperatureVector::GetTemperatureVector (TemperatureProperties const& temp) noexcept : m_type{temp.m_type} + , m_temperature{temp} + , m_q_e_over_mc2{temp.m_q_e_over_mc2} #if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) , m_from_file{temp.m_u_std_x_reader.get(), temp.m_u_std_y_reader.get(), @@ -39,4 +48,10 @@ GetTemperatureVector::GetTemperatureVector (TemperatureProperties const& temp) n m_uy_std_parser = temp.m_ptr_uy_std_parser->compile<3>(); m_uz_std_parser = temp.m_ptr_uz_std_parser->compile<3>(); } +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ + !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) + else if (m_type == TempFromFileValue) { + m_temperature_in_eV_from_file = temp.m_temperature_in_eV_reader->getView(); + } +#endif } diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index 4db2e5ec8ab..7047404a894 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -264,7 +264,8 @@ struct InjectorMomentumJuttner using namespace amrex::literals; amrex::Real x1, x2, gamma; amrex::Real u [3]; - amrex::Real const theta = temperature(x,y,z); + amrex::Real const theta = + temperature(x,y,z) * temperature.m_q_e_over_mc2; if (theta < 0) { amrex::Abort("Temperature parameter theta is less than zero, " "which is not allowed for Maxwell-Juttner distribution."); diff --git a/Source/Initialization/TemperatureProperties.H b/Source/Initialization/TemperatureProperties.H index 48de42c6c84..52c9bc400df 100644 --- a/Source/Initialization/TemperatureProperties.H +++ b/Source/Initialization/TemperatureProperties.H @@ -25,7 +25,8 @@ enum TemperatureInitType { TempParserFunction, TempParserFunctionVector, TempConstantVector, - TempFromFileVector + TempFromFileVector, + TempFromFileValue }; /** @@ -51,7 +52,7 @@ struct TemperatureProperties /* Type of temperature initialization */ TemperatureInitType m_type; - /* Constant temperature value, if m_type == TempConstantValue */ + /* Constant temperature_in_eV, if m_type == TempConstantValue */ amrex::Real m_temperature; /* Storage of the parser function, if m_type == TempParserFunction */ std::unique_ptr m_ptr_temperature_parser; @@ -62,12 +63,16 @@ struct TemperatureProperties std::unique_ptr m_ptr_ux_std_parser, m_ptr_uy_std_parser, m_ptr_uz_std_parser; std::string m_read_u_std_path; + std::string m_read_temperature_in_eV_path; + /* q_e / (m c^2) factor, used for temperature_in_eV conversion */ + amrex::Real m_q_e_over_mc2{0.0}; #if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) std::unique_ptr m_u_std_x_reader; std::unique_ptr m_u_std_y_reader; std::unique_ptr m_u_std_z_reader; + std::unique_ptr m_temperature_in_eV_reader; #endif }; diff --git a/Source/Initialization/TemperatureProperties.cpp b/Source/Initialization/TemperatureProperties.cpp index e25dbd2f792..ee5cc173d49 100644 --- a/Source/Initialization/TemperatureProperties.cpp +++ b/Source/Initialization/TemperatureProperties.cpp @@ -8,124 +8,217 @@ #include "TemperatureProperties.H" #include "ExternalField.H" +#include "Particles/SpeciesPhysicalProperties.H" #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" +#include "Utils/WarpXConst.H" #include "WarpX.H" #include #include #include +#include #include +#include + +namespace { +/** Parse shared ``temperature_in_eV`` for ``maxwellian`` / ``maxwell_juttner`` (constant, parser, or file). */ +void +parse_temperature_in_eV ( + amrex::ParmParse const& pp, + std::string const& source_name, + [[maybe_unused]] amrex::Geometry const& geom, + std::string const& dist_type_param, + std::string const& mom_dist_s, + TemperatureProperties& temp) +{ + amrex::Real mass = 0.0; + std::string physical_species_s; + if (pp.query("species_type", physical_species_s)) { + const auto physical_species_from_string = species::from_string(physical_species_s); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(physical_species_from_string, + physical_species_s + " does not exist!"); + mass = species::get_mass(physical_species_from_string.value()); + } + utils::parser::queryWithParser(pp, "mass", mass); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(mass > 0.0, + "Need to specify species_type or mass > 0.0 for the " + mom_dist_s + + " temperature in eV initialization"); + + std::string temperature_in_eV_dist_s = "constant"; + utils::parser::query(pp, source_name, dist_type_param.c_str(), temperature_in_eV_dist_s); + + const amrex::Real q_e_over_mc2 = + PhysConst::q_e / (mass * PhysConst::c * PhysConst::c); + + if (temperature_in_eV_dist_s == "constant") { + amrex::Real temperature_in_eV = 0.0; + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + utils::parser::queryWithParser(pp, source_name, "temperature_in_eV", + temperature_in_eV), + "Temperature parameter temperature_in_eV not specified"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(temperature_in_eV >= 0.0, + "temperature_in_eV = " + std::to_string(temperature_in_eV) + + " is less than zero, which is not allowed"); + temp.m_temperature = temperature_in_eV; + temp.m_q_e_over_mc2 = q_e_over_mc2; + temp.m_type = TempConstantValue; + } + else if (temperature_in_eV_dist_s == "parser") { + std::string str_temperature_in_eV_function; + utils::parser::Store_parserString( + pp, source_name, "temperature_in_eV_function(x,y,z)", + str_temperature_in_eV_function); + temp.m_ptr_temperature_parser = std::make_unique( + utils::parser::makeParser(str_temperature_in_eV_function, {"x", "y", "z"})); + temp.m_q_e_over_mc2 = q_e_over_mc2; + temp.m_type = TempParserFunction; + } + else if (temperature_in_eV_dist_s == "read_from_file") { +#if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ + !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) + if (WarpX::gamma_boost > 1.0) { + WARPX_ABORT_WITH_MESSAGE( + dist_type_param + " = read_from_file is not " + "supported in boosted-frame simulations yet."); + } + utils::parser::get(pp, source_name, "read_temperature_in_eV_from_path", + temp.m_read_temperature_in_eV_path); + std::string field_name = "temperature_in_eV"; + utils::parser::query(pp, source_name, "temperature_in_eV_mesh_name", field_name); + amrex::GpuArray const problo = + geom.ProbLoArray(); + amrex::GpuArray const dx = + geom.CellSizeArray(); + amrex::Box const dombox = amrex::convert(geom.Domain(), amrex::IntVect(1)); + temp.m_temperature_in_eV_reader = std::make_unique( + temp.m_read_temperature_in_eV_path, field_name, "", problo, dx, dombox, false); + amrex::BoxArray const grids; + amrex::DistributionMapping const dmap; + temp.m_temperature_in_eV_reader->prepare(grids, dmap, amrex::IntVect(0)); + temp.m_q_e_over_mc2 = q_e_over_mc2; + temp.m_type = TempFromFileValue; +#else + WARPX_ABORT_WITH_MESSAGE( + dist_type_param + " = read_from_file requires " + "WarpX built with openPMD support and is not supported in " + "RZ/RCYLINDER/RSPHERE geometries."); +#endif + } + else { + std::stringstream ss; + ss << mom_dist_s << " temperature distribution type '" + << temperature_in_eV_dist_s << "' not recognized."; + WARPX_ABORT_WITH_MESSAGE(ss.str()); + } +} + +} // namespace /** Construct TemperatureProperties from the passed particle source parameters. * Parse the momentum distribution type and initialize the corresponding - * temperature parameters: thermal spread `ux_std`, `uy_std`, `uz_std` - * for `maxwellian` distribution, and `theta` for `maxwell_juttner`. + * temperature parameters: thermal spread `ux_std`, `uy_std`, `uz_std` or `temperature_in_eV` for the Maxwellian distributions; + * and `temperature_in_eV` for `maxwell_juttner` distributions. */ TemperatureProperties::TemperatureProperties (const amrex::ParmParse& pp, std::string const& source_name, amrex::Geometry const& geom) { - amrex::ignore_unused(geom); - std::string mom_dist_s; utils::parser::query(pp, source_name, "momentum_distribution_type", mom_dist_s); if (mom_dist_s == "maxwell_juttner") { - // Set defaults - amrex::Real theta = 0; // quiet GCC warning maybe-uninitialized - std::string temp_dist_s = "constant"; - utils::parser::query(pp, source_name, "theta_distribution_type", temp_dist_s); - - if (temp_dist_s == "constant") { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - utils::parser::queryWithParser(pp, source_name, "theta", theta), - "Temperature parameter theta not specified"); + parse_temperature_in_eV( + pp, source_name, geom, + "maxwell_juttner_temperature_in_eV_distribution_type", + mom_dist_s, + *this); + } + else if (mom_dist_s == "maxwellian") { + const bool use_temperature_in_eV = + pp.contains("maxwellian_temperature_in_eV_distribution_type") || + (!source_name.empty() && + pp.contains(source_name + ".maxwellian_temperature_in_eV_distribution_type")); + const bool use_u_std = + pp.contains("maxwellian_u_std_distribution_type") || + (!source_name.empty() && + pp.contains(source_name + ".maxwellian_u_std_distribution_type")); - // Do validation on theta value. - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta >= 0, - "Temperature parameter theta = " + std::to_string(theta) + - " is less than zero, which is not allowed"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !(use_temperature_in_eV && use_u_std), + "Cannot specify both maxwellian_u_std_distribution_type and " + "maxwellian_temperature_in_eV_distribution_type."); - m_type = TempConstantValue; - m_temperature = theta; - } - else if (temp_dist_s == "parser") { - std::string str_theta_function; - utils::parser::Store_parserString(pp, source_name, "theta_function(x,y,z)", str_theta_function); - m_ptr_temperature_parser = - std::make_unique( - utils::parser::makeParser(str_theta_function,{"x","y","z"})); - m_type = TempParserFunction; + if (use_temperature_in_eV) { + parse_temperature_in_eV( + pp, source_name, geom, + "maxwellian_temperature_in_eV_distribution_type", + mom_dist_s, + *this); } else { - std::stringstream ss; - ss << "Temperature distribution type '" << temp_dist_s << "' not recognized."; - WARPX_ABORT_WITH_MESSAGE(ss.str()); - } - } - else if (mom_dist_s == "maxwellian") { - // ``maxwellian`` distribution uses ``u_std_*`` - std::string u_std_dist_s = "constant"; - utils::parser::query(pp, source_name, "maxwellian_u_std_distribution_type", u_std_dist_s); + // ``maxwellian`` distribution uses ``u_std_*`` + std::string u_std_dist_s = "constant"; + utils::parser::query(pp, source_name, "maxwellian_u_std_distribution_type", u_std_dist_s); - if (u_std_dist_s == "constant") { - utils::parser::queryWithParser(pp, source_name, "ux_std", m_ux_std); - utils::parser::queryWithParser(pp, source_name, "uy_std", m_uy_std); - utils::parser::queryWithParser(pp, source_name, "uz_std", m_uz_std); - m_type = TempConstantVector; - } - else if (u_std_dist_s == "parser") { - std::string sx, sy, sz; - utils::parser::Store_parserString(pp, source_name, "ux_std_function(x,y,z)", sx); - utils::parser::Store_parserString(pp, source_name, "uy_std_function(x,y,z)", sy); - utils::parser::Store_parserString(pp, source_name, "uz_std_function(x,y,z)", sz); - m_ptr_ux_std_parser = - std::make_unique(utils::parser::makeParser(sx, {"x", "y", "z"})); - m_ptr_uy_std_parser = - std::make_unique(utils::parser::makeParser(sy, {"x", "y", "z"})); - m_ptr_uz_std_parser = - std::make_unique(utils::parser::makeParser(sz, {"x", "y", "z"})); - m_type = TempParserFunctionVector; - } - else if (u_std_dist_s == "read_from_file") { + if (u_std_dist_s == "constant") { + utils::parser::queryWithParser(pp, source_name, "ux_std", m_ux_std); + utils::parser::queryWithParser(pp, source_name, "uy_std", m_uy_std); + utils::parser::queryWithParser(pp, source_name, "uz_std", m_uz_std); + m_type = TempConstantVector; + } + else if (u_std_dist_s == "parser") { + std::string sx, sy, sz; + utils::parser::Store_parserString(pp, source_name, "ux_std_function(x,y,z)", sx); + utils::parser::Store_parserString(pp, source_name, "uy_std_function(x,y,z)", sy); + utils::parser::Store_parserString(pp, source_name, "uz_std_function(x,y,z)", sz); + m_ptr_ux_std_parser = + std::make_unique(utils::parser::makeParser(sx, {"x", "y", "z"})); + m_ptr_uy_std_parser = + std::make_unique(utils::parser::makeParser(sy, {"x", "y", "z"})); + m_ptr_uz_std_parser = + std::make_unique(utils::parser::makeParser(sz, {"x", "y", "z"})); + m_type = TempParserFunctionVector; + } + else if (u_std_dist_s == "read_from_file") { #if defined(WARPX_USE_OPENPMD) && !defined(WARPX_DIM_RZ) && \ !defined(WARPX_DIM_RCYLINDER) && !defined(WARPX_DIM_RSPHERE) - if (WarpX::gamma_boost > 1.0) { - WARPX_ABORT_WITH_MESSAGE( - "maxwellian_u_std_distribution_type = read_from_file is not " - "supported in boosted-frame simulations yet."); - } - utils::parser::get(pp, source_name, "read_u_std_from_path", m_read_u_std_path); - amrex::GpuArray const problo = - geom.ProbLoArray(); - amrex::GpuArray const dx = - geom.CellSizeArray(); - amrex::Box const dombox = amrex::convert(geom.Domain(), amrex::IntVect(1)); - m_u_std_x_reader = std::make_unique( - m_read_u_std_path, "u_std", "x", problo, dx, dombox, false); - m_u_std_y_reader = std::make_unique( - m_read_u_std_path, "u_std", "y", problo, dx, dombox, false); - m_u_std_z_reader = std::make_unique( - m_read_u_std_path, "u_std", "z", problo, dx, dombox, false); - amrex::BoxArray const grids; - amrex::DistributionMapping const dmap; - m_u_std_x_reader->prepare(grids, dmap, amrex::IntVect(0)); - m_u_std_y_reader->prepare(grids, dmap, amrex::IntVect(0)); - m_u_std_z_reader->prepare(grids, dmap, amrex::IntVect(0)); - m_type = TempFromFileVector; + if (WarpX::gamma_boost > 1.0) { + WARPX_ABORT_WITH_MESSAGE( + "maxwellian_u_std_distribution_type = read_from_file is not " + "supported in boosted-frame simulations yet."); + } + utils::parser::get(pp, source_name, "read_u_std_from_path", m_read_u_std_path); + amrex::GpuArray const problo = + geom.ProbLoArray(); + amrex::GpuArray const dx = + geom.CellSizeArray(); + amrex::Box const dombox = amrex::convert(geom.Domain(), amrex::IntVect(1)); + m_u_std_x_reader = std::make_unique( + m_read_u_std_path, "u_std", "x", problo, dx, dombox, false); + m_u_std_y_reader = std::make_unique( + m_read_u_std_path, "u_std", "y", problo, dx, dombox, false); + m_u_std_z_reader = std::make_unique( + m_read_u_std_path, "u_std", "z", problo, dx, dombox, false); + amrex::BoxArray const grids; + amrex::DistributionMapping const dmap; + m_u_std_x_reader->prepare(grids, dmap, amrex::IntVect(0)); + m_u_std_y_reader->prepare(grids, dmap, amrex::IntVect(0)); + m_u_std_z_reader->prepare(grids, dmap, amrex::IntVect(0)); + m_type = TempFromFileVector; #else - WARPX_ABORT_WITH_MESSAGE( - "maxwellian_u_std_distribution_type = read_from_file requires " - "WarpX built with openPMD support and is not supported in " - "RZ/RCYLINDER/RSPHERE geometries."); + WARPX_ABORT_WITH_MESSAGE( + "maxwellian_u_std_distribution_type = read_from_file requires " + "WarpX built with openPMD support and is not supported in " + "RZ/RCYLINDER/RSPHERE geometries."); #endif - } - else { - std::stringstream ss; - ss << "Maxwellian velocity standard deviation distribution type '" << u_std_dist_s - << "' not recognized."; - WARPX_ABORT_WITH_MESSAGE(ss.str()); + } + else { + std::stringstream ss; + ss << "Maxwellian velocity standard deviation distribution type '" << u_std_dist_s + << "' not recognized."; + WARPX_ABORT_WITH_MESSAGE(ss.str()); + } } } else { diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 848889281a6..330d77bf00f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -464,6 +464,28 @@ PhysicalParticleContainer::BackwardCompatibility () WARPX_ABORT_WITH_MESSAGE(msg); } } + + std::string mom_dist_s; + if (pp_species_name.query("momentum_distribution_type", mom_dist_s) && + mom_dist_s == "maxwell_juttner") { + const std::string juttner_temp_msg = + "Maxwell-Juttner thermal spread is now specified with temperature_in_eV. " + "Use .maxwell_juttner_temperature_in_eV_distribution_type = constant " + "(default), parser, or read_from_file, and provide .temperature_in_eV, " + ".temperature_in_eV_function(x,y,z), or " + ".read_temperature_in_eV_from_path. " + "Requires species_type or mass."; + for (const std::string old_param : + {"theta", "theta_distribution_type", "theta_function(x,y,z)"}) { + if (pp_species_name.query(old_param, backward_string)) { + std::string msg = "."; + msg += old_param; + msg += " is no longer supported. "; + msg += juttner_temp_msg; + WARPX_ABORT_WITH_MESSAGE(msg); + } + } + } } void PhysicalParticleContainer::InitData ()