Replace theta with temperature_in_eV for the Maxwell-Juttner momentum distribution

Docs/source/usage/parameters.rst

@@ -1807,10 +1807,13 @@ Particle initialization
 
       * ``<species_name>.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
+        ``<species_name>.maxwellian_temperature_in_eV_distribution_type``.
 
         * If ``constant``, the following are required: ``<species_name>.ux_std``,
           ``<species_name>.uy_std``, ``<species_name>.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.
 
+      * ``<species_name>.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 ``<species_name>.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: ``<species_name>.temperature_in_eV`` (`float`).
+        * If ``parser``, the following is required:
+          ``<species_name>.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:
+          ``<species_name>.read_temperature_in_eV_from_path`` (openPMD file path). The file must
+          contain a scalar openPMD mesh with the name given by
+          ``<species_name>.temperature_in_eV_mesh_name`` (default ``temperature_in_eV``). See
+          `this file <https://github.com/BLAST-WarpX/warpx/blob/development/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py>`__
+          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 <https://github.com/BLAST-WarpX/warpx/blob/development/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py>`__
           for an example of how to prepare the openPMD data file.
 
-      * ``<species_name>.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: ``<species_name>.theta`` (`float`).
-        * If ``parser``, the following is required: ``<species_name>.theta_function(x,y,z)``.
+      * ``<species_name>.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: ``<species_name>.temperature_in_eV`` (`float`).
+        * If ``parser``, the following is required:
+          ``<species_name>.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:
+          ``<species_name>.read_temperature_in_eV_from_path`` (openPMD file path). The file must
+          contain a scalar openPMD mesh with the name given by
+          ``<species_name>.temperature_in_eV_mesh_name`` (default ``temperature_in_eV``). See
+          `this file <https://github.com/BLAST-WarpX/warpx/blob/development/Examples/Tests/initial_distribution/inputs_test_3d_initial_distribution_prepare.py>`__
+          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

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