Update constants to use values from CODATA 2022

.azure-pipelines.yml

@@ -1,7 +1,7 @@
 # -*- mode: yaml -*-
 
 pool:
-  vmImage: 'ubuntu-20.04'
+  vmImage: 'ubuntu-24.04'
 
 trigger:
   branches:
@@ -111,6 +111,10 @@ jobs:
           -DCMAKE_CXX_STANDARD=17                         \
           -Duse_cmake_find_lapack=ON -Dbuild_tests=OFF -DCMAKE_VERBOSE_MAKEFILE=ON
       fi
+      # Remove system copy of Matplotlib to avoid conflict
+      # with version set in the requirements file - see, e.g.,
+      # https://github.com/matplotlib/matplotlib/issues/28768.
+      sudo apt remove python3-matplotlib
       # Python modules required for test analysis
       python3 -m pip install --upgrade -r Regression/requirements.txt
       python3 -m pip cache purge

Docs/source/usage/workflows/ml_materials/run_warpx_training.py

@@ -107,7 +107,6 @@ def get_species_of_accelerator_stage(
         density_expression=f"n0*(1.+4.*(x**2+y**2)/(kp**2*Rc**4))*(0.5*(1.-cos(pi*(z-{stage_zmin})/Lplus)))*((z-{stage_zmin})<Lplus)"
         + f"+n0*(1.+4.*(x**2+y**2)/(kp**2*Rc**4))*((z-{stage_zmin})>=Lplus)*((z-{stage_zmin})<(Lplus+Lp))"
         + f"+n0*(1.+4.*(x**2+y**2)/(kp**2*Rc**4))*(0.5*(1.+cos(pi*((z-{stage_zmin})-Lplus-Lp)/Lminus)))*((z-{stage_zmin})>=(Lplus+Lp))*((z-{stage_zmin})<(Lplus+Lp+Lminus))",
-        pi=3.141592653589793,
         n0=n0,
         kp=q_e / c * math.sqrt(n0 / (m_e * ep0)),
         Rc=Rc,

Examples/Tests/collision/analysis_collision_1d.py

@@ -19,7 +19,7 @@
 
 import numpy as np
 import yt
-from scipy.constants import e
+from scipy.constants import e, m_e, m_u
 
 # this will be the name of the plot file
 last_fn = sys.argv[1]
@@ -29,7 +29,7 @@
 )
 
 # carbon 12 ion (mass = 12*amu - 6*me)
-mass = 1.992100316897910e-26
+mass = 12.0 * m_u - 6.0 * m_e
 
 # Separate macroparticles from group A (low weight) and group B (high weight)
 # by sorting based on weight

Examples/Tests/collision/analysis_collision_2d.py

@@ -31,6 +31,7 @@
 import numpy
 import post_processing_utils
 import yt
+from scipy.constants import c, m_e
 
 test_name = os.path.split(os.getcwd())[1]
 
@@ -41,9 +42,7 @@
 ni = ng * 200
 np = ne + ni
 
-c = 299792458.0
-me = 9.10938356e-31
-mi = me * 5.0
+mi = m_e * 5.0
 
 ## In the first part of the test we verify that the output data is consistent with the exponential
 ## fit.
@@ -76,7 +75,7 @@
     # get time index j
     j = int(fn[-5:])
     # compute error
-    vxe = numpy.mean(px[0:ne]) / me / c
+    vxe = numpy.mean(px[0:ne]) / m_e / c
     vxi = numpy.mean(px[ne:np]) / mi / c
     vxd = vxe - vxi
     fit = a * math.exp(b * j)

Examples/Tests/collision/analysis_collision_3d.py

@@ -30,6 +30,7 @@
 import numpy
 import post_processing_utils
 import yt
+from scipy.constants import c, m_e
 
 tolerance = 0.001
 
@@ -38,9 +39,7 @@
 ni = ng * 200
 np = ne + ni
 
-c = 299792458.0
-me = 9.10938356e-31
-mi = me * 5.0
+mi = m_e * 5.0
 
 ## In the first part of the test we verify that the output data is consistent with the exponential
 ## fit.
@@ -69,7 +68,7 @@
     # get time index j
     j = int(fn[-5:])
     # compute error
-    vxe = numpy.mean(pxe) / me / c
+    vxe = numpy.mean(pxe) / m_e / c
     vxi = numpy.mean(pxi) / mi / c
     vxd = vxe - vxi
     fit = a * math.exp(b * j)

Examples/Tests/collision/inputs_test_2d_collision_xz

@@ -47,7 +47,7 @@ electron.uz_th = 0.044237441120300
 electron.ux_m  = 0.044237441120300
 
 ion.charge = q_e
-ion.mass = 4.554691780000000e-30
+ion.mass = 5*m_e
 ion.injection_style = "NRandomPerCell"
 ion.num_particles_per_cell = 200
 ion.profile = constant

Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py

@@ -25,7 +25,7 @@
 import numpy as np
 import yt
 from openpmd_viewer import OpenPMDTimeSeries
-from scipy.constants import c
+from scipy.constants import c, e, epsilon_0, m_e
 from scipy.optimize import fsolve
 
 yt.funcs.mylog.setLevel(0)
@@ -41,10 +41,10 @@
 
 if emass_10:
     l2_tolerance = 0.096
-    m_e = 10
+    e_mass = 10
 else:
     l2_tolerance = 0.05
-    m_e = 9.10938356e-31  # Electron mass in kg
+    e_mass = m_e  # Electron mass in kg
 ndims = np.count_nonzero(ds.domain_dimensions > 1)
 
 if ndims == 2:
@@ -68,8 +68,7 @@
 iz0 = round((0.0 - zmin) / dz)
 
 # Constants
-eps_0 = 8.8541878128e-12  # Vacuum Permittivity in C/(V*m)
-q_e = -1.60217662e-19  # Electron charge in C
+q_e = -e  # Electron charge in C
 pi = np.pi  # Circular constant of the universe
 r_0 = 0.1  # Initial radius of sphere
 q_tot = -1e-15  # Total charge of sphere in C
@@ -81,15 +80,15 @@
 # v(r) and t(r) can be solved analytically.
 #
 # The solution r(t) solves the ODE: r''(t) = a/(r(t)**2) with initial conditions
-# r(0) = r_0, r'(0) = 0, and a = q_e*q_tot/(4*pi*eps_0*m_e)
+# r(0) = r_0, r'(0) = 0, and a = q_e*q_tot/(4*pi*epsilon_0*e_mass)
 #
 # The E was calculated at the end of the last time step
 def v_exact(r):
-    return np.sqrt(q_e * q_tot / (2 * pi * m_e * eps_0) * (1 / r_0 - 1 / r))
+    return np.sqrt(q_e * q_tot / (2 * pi * e_mass * epsilon_0) * (1 / r_0 - 1 / r))
 
 
 def t_exact(r):
-    return np.sqrt(r_0**3 * 2 * pi * m_e * eps_0 / (q_e * q_tot)) * (
+    return np.sqrt(r_0**3 * 2 * pi * e_mass * epsilon_0 / (q_e * q_tot)) * (
         np.sqrt(r / r_0 - 1) * np.sqrt(r / r_0)
         + np.log(np.sqrt(r / r_0 - 1) + np.sqrt(r / r_0))
     )
@@ -104,8 +103,8 @@ def func(rho):
 
 def E_exact(r):
     return np.sign(r) * (
-        q_tot / (4 * pi * eps_0 * r**2) * (abs(r) >= r_end)
-        + q_tot * abs(r) / (4 * pi * eps_0 * r_end**3) * (abs(r) < r_end)
+        q_tot / (4 * pi * epsilon_0 * r**2) * (abs(r) >= r_end)
+        + q_tot * abs(r) / (4 * pi * epsilon_0 * r_end**3) * (abs(r) < r_end)
     )
 
 

Examples/Tests/embedded_boundary_cube/inputs_test_3d_embedded_boundary_cube_macroscopic

@@ -3,6 +3,6 @@ FILE = inputs_base_3d
 
 # test input parameters
 algo.em_solver_medium = macroscopic
-macroscopic.epsilon = 1.5*8.8541878128e-12
-macroscopic.mu = 1.25663706212e-06
+macroscopic.epsilon = 1.5*epsilon0
+macroscopic.mu = mu0
 macroscopic.sigma = 0

Examples/Tests/embedded_boundary_rotated_cube/inputs_test_2d_embedded_boundary_rotated_cube

@@ -18,7 +18,6 @@ my_constants.xmin = -0.53
 my_constants.zmin = -0.53
 my_constants.xmax = 0.53
 my_constants.zmax = 0.53
-my_constants.pi = 3.141592653589793
 my_constants.theta = pi/8
 
 warpx.eb_implicit_function = "xr=x*cos(-theta)+z*sin(-theta); zr=-x*sin(-theta)+z*cos(-theta); max(max(xr+xmin,-(xr+xmax)), max(zr+zmin,-(zr+zmax)))"
@@ -29,12 +28,11 @@ my_constants.Lx = 1.06
 my_constants.Lz = 1.06
 my_constants.x_cent = 0.
 my_constants.z_cent = 0.
-my_constants.mu_0 = 1.25663706212e-06
 
 warpx.B_ext_grid_init_style = parse_B_ext_grid_function
 
 warpx.Bx_external_grid_function(x,y,z) = 0
-warpx.By_external_grid_function(x,y,z) = "mu_0 *
+warpx.By_external_grid_function(x,y,z) = "mu0 *
                                           cos(m * pi / Lx * (x*cos(-theta)+z*sin(-theta) - Lx / 2 - x_cent)) *
                                           cos(p * pi / Lz * (-x*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent))"
 warpx.Bz_external_grid_function(x,y,z) = 0

Examples/Tests/embedded_boundary_rotated_cube/inputs_test_3d_embedded_boundary_rotated_cube

@@ -20,7 +20,6 @@ my_constants.zmin = -0.5
 my_constants.xmax = 0.5
 my_constants.ymax = 0.5
 my_constants.zmax = 0.5
-my_constants.pi = 3.141592653589793
 my_constants.theta = pi/6
 
 warpx.eb_implicit_function = "max(max(max(x+xmin,-(x+xmax)), max(y*cos(-theta)-z*sin(-theta)+ymin,-(y*cos(-theta)-z*sin(-theta)+ymax))), max(y*sin(-theta)+z*cos(-theta)+zmin,-(y*sin(-theta)+z*cos(-theta)+zmax)))"
@@ -35,31 +34,30 @@ my_constants.x_cent = 0.
 my_constants.y_cent = 0.
 my_constants.z_cent = 0.
 my_constants.h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
-my_constants.mu_0 = 1.25663706212e-06
 
 warpx.B_ext_grid_init_style = parse_B_ext_grid_function
-warpx.Bx_external_grid_function(x,y,z) = "-2/h_2 * mu_0 * (m * pi / Lx) * (p * pi / Lz) *
+warpx.Bx_external_grid_function(x,y,z) = "-2/h_2 * mu0 * (m * pi / Lx) * (p * pi / Lz) *
                                           sin(m * pi / Lx * (x - Lx / 2 - x_cent)) *
                                           cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
                                           cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent))"
 
-warpx.By_external_grid_function(x,y,z) = "-2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) *
+warpx.By_external_grid_function(x,y,z) = "-2/h_2 * mu0 * (n * pi / Ly) * (p * pi / Lz) *
                                           cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
                                           sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
                                           cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
                                           cos(theta) -
-                                          mu_0 *
+                                          mu0 *
                                           cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
                                           cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
                                           sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
                                           sin(theta)"
 
-warpx.Bz_external_grid_function(x,y,z) = "mu_0 *
+warpx.Bz_external_grid_function(x,y,z) = "mu0 *
                                           cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
                                           cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
                                           sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *
                                           cos(theta) -
-                                          2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) *
+                                          2/h_2 * mu0 * (n * pi / Ly) * (p * pi / Lz) *
                                           cos(m * pi / Lx * (x - Lx / 2 - x_cent)) *
                                           sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) *
                                           cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) *

Examples/Tests/implicit/analysis_2d_psatd.py

@@ -22,7 +22,7 @@
 max_delta_E = np.abs(delta_E).max()
 
 # This case should have near machine precision conservation of energy
-tolerance_rel_energy = 2.1e-14
+tolerance_rel_energy = 2.4e-14
 
 print(f"max change in energy: {max_delta_E}")
 print(f"tolerance: {tolerance_rel_energy}")

Examples/Tests/initial_distribution/analysis.py

@@ -134,7 +134,7 @@
 x_rms = 0.25
 z_cut = 2.0
 q_tot = -1.0e-20
-q_e = -1.602176634e-19
+q_e = -scc.e
 npart = q_tot / q_e
 db = bin_value[1] - bin_value[0]
 

Examples/Tests/ionization_dsmc/inputs_base_3d

@@ -104,6 +104,8 @@ counts.intervals = 5
 counts.path = diags/
 counts.type = ParticleNumber
 
+warpx.synchronize_velocity_for_diagnostics = 1
+
 # Diagnostics
 diag1.intervals = 250
 diag1.diag_type = Full

Examples/Tests/load_external_field/analysis_3d.py

@@ -23,7 +23,7 @@
 
 tolerance = 1.0e-8
 x0 = 0.12238072
-y0 = 0.00965394
+y0 = 0.00965395
 z0 = 4.31754477
 
 filename = sys.argv[1]

Examples/Tests/nuclear_fusion/analysis_proton_boron_fusion.py

@@ -80,7 +80,7 @@
 
 ## Checks whether this is the 2D or the 3D test
 with open("./warpx_used_inputs") as warpx_used_inputs:
-    is_2D = re.search("geometry.dims\s*=\s*2", warpx_used_inputs.read())
+    is_2D = re.search(r"geometry.dims\s*=\s*2", warpx_used_inputs.read())
 warpx_used_inputs.close()
 
 ## Some numerical parameters for this test
@@ -697,8 +697,11 @@ def check_initial_energy2(data):
         ## Tolerance is quite high because we don't have a lot of alphas to produce good
         ## statistics and an event like alpha1 emitted exactly in direction of proton & alpha2
         ## emitted exactly in direction opposite to Beryllium is somewhat rare.
+        print(
+            f"Check energy max: {np.amax(energy_alpha2_simulation)} {max_energy_alpha23} {(np.amax(energy_alpha2_simulation) - max_energy_alpha23) / max_energy_alpha23}"
+        )
         assert is_close(
-            np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=5e-2
+            np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=5.1e-2
         )
         assert is_close(
             np.amin(energy_alpha2_simulation), min_energy_alpha23, atol=3.218e-14

Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py

@@ -82,7 +82,7 @@
     "helium4": 4.00260325413 * scc.m_u,
     "neutron": 1.0013784193052508 * scc.m_p,
 }
-m_reduced = np.product([mass[s] for s in reactant_species]) / np.sum(
+m_reduced = np.prod([mass[s] for s in reactant_species]) / np.sum(
     [mass[s] for s in reactant_species]
 )
 

Examples/Tests/ohm_solver_em_modes/analysis_rz.py

@@ -127,13 +127,13 @@ def process(it):
     omega_fast = sim.vA * np.sqrt(R2 + np.sqrt(R2**2 - P4))
     omega_slow = sim.vA * np.sqrt(R2 - np.sqrt(R2**2 - P4))
     # Upper right corner
-    ax.plot(k * sim.l_i, omega_fast / sim.w_ci, "w--", label="$\omega_{fast}$")
+    ax.plot(k * sim.l_i, omega_fast / sim.w_ci, "w--", label=r"$\omega_{fast}$")
     ax.plot(
         k * sim.l_i,
         omega_slow / sim.w_ci,
         color="white",
         linestyle="--",
-        label="$\omega_{slow}$",
+        label=r"$\omega_{slow}$",
     )
     # Thermal resonance
     thermal_res = sim.w_ci + 3 * sim.v_ti * k
@@ -142,7 +142,7 @@ def process(it):
         thermal_res / sim.w_ci,
         color="magenta",
         linestyle="--",
-        label="$\omega = \Omega_i + 3v_{th,i}k$",
+        label=r"$\omega = \Omega_i + 3v_{th,i}k$",
     )
     ax.plot(
         -k * sim.l_i, thermal_res / sim.w_ci, color="magenta", linestyle="--", label=""
@@ -153,7 +153,7 @@ def process(it):
         thermal_res / sim.w_ci,
         color="magenta",
         linestyle="--",
-        label="$\omega = \Omega_i + 3v_{th,i}k$",
+        label=r"$\omega = \Omega_i + 3v_{th,i}k$",
     )
     ax.plot(
         -k * sim.l_i, thermal_res / sim.w_ci, color="magenta", linestyle="--", label=""
@@ -164,10 +164,10 @@ def process(it):
     ax.set_xlim(-1.75, 1.75)
     ax.set_ylim(0, 1.6)
 
-axes[0, 0].set_ylabel("$\omega/\Omega_{ci}$")
-axes[1, 0].set_ylabel("$\omega/\Omega_{ci}$")
-axes[1, 0].set_xlabel("$k_zl_i$")
-axes[1, 1].set_xlabel("$k_zl_i$")
+axes[0, 0].set_ylabel(r"$\omega/\Omega_{ci}$")
+axes[1, 0].set_ylabel(r"$\omega/\Omega_{ci}$")
+axes[1, 0].set_xlabel(r"$k_zl_i$")
+axes[1, 1].set_xlabel(r"$k_zl_i$")
 
 plt.savefig("normal_modes_disp.png", dpi=600)
 if not sim.test:
@@ -179,5 +179,5 @@ def process(it):
     amps = np.abs(F_kw[2, 1, len(kz) // 2 - 2 : len(kz) // 2 + 2])
     print("Amplitude sample: ", amps)
     assert np.allclose(
-        amps, np.array([59.23850009, 19.26746169, 92.65794174, 10.83627164])
+        amps, np.array([79.98358457, 27.24783417, 213.15227656, 24.12800388])
     )

Examples/Tests/ohm_solver_ion_beam_instability/analysis.py

@@ -98,7 +98,7 @@
     )
 
     # The theoretical growth rates for the 4th, 5th and 6th Fourier modes of
-    # the By-field was obtained from Fig. 12a of Munoz et al.
+    # the By-field was obtained from Fig. 12a of Munoz et al., https://doi.org/10.1016/j.cpc.2017.10.012
     # Note the rates here are gamma / w_ci
     gamma4 = 0.1915611861780133
     gamma5 = 0.20087036355662818
@@ -227,6 +227,6 @@
     # compared to the theoretical ones to determine if the physics test passes.
     # At creation, the full test (3d) had the following errors (ran on 1 V100):
     # m4_rms_error = 3.329; m5_rms_error = 1.052; m6_rms_error = 2.583
-    assert np.isclose(m4_rms_error, 1.515, atol=0.01)
-    assert np.isclose(m5_rms_error, 0.718, atol=0.01)
-    assert np.isclose(m6_rms_error, 0.357, atol=0.01)
+    assert np.isclose(m4_rms_error, 1.546, atol=0.01)
+    assert np.isclose(m5_rms_error, 0.734, atol=0.01)
+    assert np.isclose(m6_rms_error, 0.367, atol=0.01)