Merge branch 'ECP-WarpX:development' into add_feature_projection_div_…

…cleaner

.github/workflows/clang_tidy.yml

@@ -13,6 +13,7 @@ jobs:
         dim: [1, 2, RZ, 3]
     name: clang-tidy-${{ matrix.dim }}D
     runs-on: ubuntu-22.04
+    timeout-minutes: 120
     if: github.event.pull_request.draft == false
     steps:
     - uses: actions/checkout@v4

.github/workflows/cuda.yml

@@ -115,7 +115,7 @@ jobs:
         which nvcc || echo "nvcc not in PATH!"
 
         git clone https://github.com/AMReX-Codes/amrex.git ../amrex
-        cd ../amrex && git checkout --detach dcb9cc0383dcc71e38dee9070574e325a812f8bf && cd -
+        cd ../amrex && git checkout --detach 20e6f2eadf0c297517588ba38973ec7c7084fa31 && cd -
         make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_FFT=TRUE USE_CCACHE=TRUE -j 4
 
         ccache -s

Docs/source/highlights.rst

@@ -14,6 +14,11 @@ Plasma-Based Acceleration
 
 Scientific works in laser-plasma and beam-plasma acceleration.
 
+#. Shrock JE, Rockafellow E, Miao B, Le M, Hollinger RC, Wang S, Gonsalves AJ, Picksley A, Rocca JJ, and Milchberg HM
+   **Guided Mode Evolution and Ionization Injection in Meter-Scale Multi-GeV Laser Wakefield Accelerators**.
+   Phys. Rev. Lett. **133**, 045002, 2024
+   `DOI:10.1103/PhysRevLett.133.045002 <https://link.aps.org/doi/10.1103/PhysRevLett.133.045002>`__
+
 #. Ross AJ, Chappell J, van de Wetering JJ, Cowley J, Archer E, Bourgeois N, Corner L, Emerson DR, Feder L, Gu XJ, Jakobsson O, Jones H, Picksley A, Reid L, Wang W, Walczak R, Hooker SM.
    **Resonant excitation of plasma waves in a plasma channel**.
    Phys. Rev. Research **6**, L022001, 2024

Examples/Tests/collision/analysis_collision_1d.py

@@ -0,0 +1,126 @@
+#!/usr/bin/env python3
+
+# Copyright 2024 Justin Angus
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from the script `inputs_1d`.
+# run locally: python analysis_vandb_1d.py diags/diag1000600/
+#
+# This is a 1D intra-species Coulomb scattering relaxation test consisting
+# of a low-density population streaming into a higher density population at rest.
+# Both populations belong to the same carbon12 ion species.
+# See test T1b from JCP 413 (2020) by D. Higginson, et al.
+#
+import os
+import sys
+
+import numpy as np
+import yt
+from scipy.constants import e
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# this will be the name of the plot file
+fn = sys.argv[1]
+ds = yt.load(fn)
+data = ds.covering_grid(level = 0, left_edge = ds.domain_left_edge, dims = ds.domain_dimensions)
+
+# carbon 12 ion (mass = 12*amu - 6*me)
+mass = 1.992100316897910e-26
+
+# Separate macroparticles from group A (low weight) and group B (high weight)
+# by sorting based on weight
+sorted_indices = data['ions','particle_weight'].argsort()
+sorted_wp = data['ions', 'particle_weight'][sorted_indices].value
+sorted_px = data['ions', 'particle_momentum_x'][sorted_indices].value
+sorted_py = data['ions', 'particle_momentum_y'][sorted_indices].value
+sorted_pz = data['ions', 'particle_momentum_z'][sorted_indices].value
+
+# Find the index 'Npmin' that separates macroparticles from group A and group B
+Np = len(sorted_wp)
+wpmin = sorted_wp.min();
+wpmax = sorted_wp.max();
+for i in range(len(sorted_wp)):
+    if sorted_wp[i] > wpmin:
+        Npmin = i
+        break
+
+NpA = Npmin
+wpA = wpmin;
+NpB = Np - Npmin
+wpB = wpmax;
+NpAs = 0
+NpAe = Npmin
+NpBs = Npmin
+NpBe = Np
+
+#############
+
+sorted_px_sum = np.abs(sorted_px).sum();
+sorted_py_sum = np.abs(sorted_py).sum();
+sorted_pz_sum = np.abs(sorted_pz).sum();
+sorted_wp_sum = np.abs(sorted_wp).sum();
+
+# compute mean velocities
+wAtot = wpA*NpA
+wBtot = wpB*NpB
+
+uBx = uBy = uBz = 0.
+for i in range(NpBs,NpBe):
+    uBx += wpB*sorted_px[i]
+    uBy += wpB*sorted_py[i]
+    uBz += wpB*sorted_pz[i]
+uBx /= (mass*wBtot) # [m/s]
+uBy /= (mass*wBtot) # [m/s]
+uBz /= (mass*wBtot) # [m/s]
+
+uAx = uAy = uAz = 0.
+for i in range(NpAs,NpAe):
+    uAx += wpA*sorted_px[i]
+    uAy += wpA*sorted_py[i]
+    uAz += wpA*sorted_pz[i]
+uAx /= (mass*wAtot) # [m/s]
+uAy /= (mass*wAtot) # [m/s]
+uAz /= (mass*wAtot) # [m/s]
+
+# compute temperatures
+TBx = TBy = TBz = 0.
+for i in range(NpBs,NpBe):
+    TBx += wpB*(sorted_px[i]/mass - uBx)**2
+    TBy += wpB*(sorted_py[i]/mass - uBy)**2
+    TBz += wpB*(sorted_pz[i]/mass - uBz)**2
+TBx *= mass/(e*wBtot)
+TBy *= mass/(e*wBtot)
+TBz *= mass/(e*wBtot)
+
+TAx = TAy = TAz = 0.
+for i in range(NpAs,NpAe):
+    TAx += wpA*(sorted_px[i]/mass - uAx)**2
+    TAy += wpA*(sorted_py[i]/mass - uAy)**2
+    TAz += wpA*(sorted_pz[i]/mass - uAz)**2
+TAx *= mass/(e*wAtot)
+TAy *= mass/(e*wAtot)
+TAz *= mass/(e*wAtot)
+
+TApar = TAz
+TAperp = (TAx + TAy)/2.0
+TA = (TAx + TAy + TAz)/3.0
+
+TBpar = TBz
+TBperp = (TBx + TBy)/2.0
+TB = (TBx + TBy + TBz)/3.0
+
+TApar_30ps_soln = 6.15e3 # TA parallel solution at t = 30 ps
+error = np.abs(TApar-TApar_30ps_soln)/TApar_30ps_soln
+tolerance = 0.02
+print('TApar at 30ps error = ', error);
+print('tolerance = ', tolerance);
+assert error < tolerance
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn)

Examples/Tests/collision/analysis_collision_2d.py

@@ -114,4 +114,4 @@
                                           dim, species_name)
 
 test_name = os.path.split(os.getcwd())[1]
-checksumAPI.evaluate_checksum(test_name, fn, do_particles=False)
+checksumAPI.evaluate_checksum(test_name, fn)

Examples/Tests/collision/inputs_1d

@@ -0,0 +1,90 @@
+#################################
+########## CONSTANTS ############
+#################################
+
+my_constants.nA = 1.e25       # m^-3
+my_constants.NpA = 400        # m^-3
+my_constants.UA = 6.55e5      # m/s
+my_constants.TA = 500.        # eV
+#
+my_constants.nB = 1.e26       # m^-3
+my_constants.NpB = 400        # m^-3
+my_constants.UB = 0.          # m/s
+my_constants.TB = 500.        # eV
+#
+my_constants.q_c12 = 6.*q_e
+my_constants.m_c12 = 12.*m_u - 6.*m_e
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 600
+amr.n_cell = 180
+amr.max_level = 0
+amr.blocking_factor = 4
+geometry.dims = 1
+geometry.prob_lo = 0.
+geometry.prob_hi = 0.01
+
+#################################
+###### Boundary Condition #######
+#################################
+boundary.field_lo = periodic
+boundary.field_hi = periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.const_dt = 0.05e-12
+warpx.use_filter = 0
+
+# Do not evolve the E and B fields
+algo.maxwell_solver = none
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = ions
+
+ions.charge = q_c12
+ions.mass = m_c12
+ions.do_not_deposit = 1
+
+ions.injection_sources = groupA groupB
+
+ions.groupA.injection_style = "NUniformPerCell"
+ions.groupA.num_particles_per_cell_each_dim = NpA
+ions.groupA.profile = constant
+ions.groupA.density = nA  # number per m^3
+ions.groupA.momentum_distribution_type = "gaussian"
+ions.groupA.uz_m = UA/clight
+ions.groupA.ux_th = sqrt(TA*q_e/m_c12)/clight
+ions.groupA.uy_th = sqrt(TA*q_e/m_c12)/clight
+ions.groupA.uz_th = sqrt(TA*q_e/m_c12)/clight
+
+ions.groupB.injection_style = "NUniformPerCell"
+ions.groupB.num_particles_per_cell_each_dim = NpB
+ions.groupB.profile = constant
+ions.groupB.density = nB  # number per m^3
+ions.groupB.momentum_distribution_type = "gaussian"
+ions.groupB.uz_m = UB/clight
+ions.groupB.ux_th = sqrt(TB*q_e/m_c12)/clight
+ions.groupB.uy_th = sqrt(TB*q_e/m_c12)/clight
+ions.groupB.uz_th = sqrt(TB*q_e/m_c12)/clight
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = collision1
+collision1.species = ions ions
+collision1.CoulombLog = 10.0
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 600
+diag1.diag_type = Full

Examples/Tests/electrostatic_sphere_eb/analysis_rz_mr.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+# Copyright 2024 Olga Shapoval, Edoardo Zoni
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This script tests the embedded boundary in RZ.
+# A cylindrical surface (r=0.1) has a fixed potential 1 V.
+# The outer surface has 0 V fixed.
+# Thus the analytical solution has the form:
+# phi(r) = A+B*log(r), Er(r) = -B/r.
+
+import os
+import sys
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+tolerance = 0.004
+print(f'tolerance = {tolerance}')
+
+fn = sys.argv[1]
+
+def find_first_non_zero_from_bottom_left(matrix):
+    for i in range(matrix.shape[0]):
+        for j in range(matrix.shape[1]):
+            if (matrix[i][j] != 0) and (matrix[i][j] != np.nan):
+                return (i, j)
+    return i, j
+
+def find_first_non_zero_from_upper_right(matrix):
+    for i in range(matrix.shape[0]-1, -1, -1):
+        for j in range(matrix.shape[1]-1, -1, -1):
+            if (matrix[i][j] != 0) and (matrix[i][j] != np.nan):
+                return (i, j)
+    return i,j
+
+def get_fields(ts, level):
+    if level == 0:
+        Er, info = ts.get_field('E', 'r', iteration=0)
+        phi, info = ts.get_field('phi', iteration=0)
+    else:
+        Er, info = ts.get_field(f'E_lvl{level}', 'r', iteration=0)
+        phi, info = ts.get_field(f'phi_lvl{level}', iteration=0)
+    return Er, phi, info
+
+def get_error_per_lev(ts,level):
+    Er, phi, info = get_fields(ts, level)
+
+    nr_half = info.r.shape[0] // 2
+    dr = info.dr
+
+    Er_patch = Er[:,nr_half:]
+    phi_patch = phi[:,nr_half:]
+    r1 = info.r[nr_half:]
+    patch_left_lower_i, patch_left_lower_j = find_first_non_zero_from_bottom_left(Er_patch)
+    patch_right_upper_i, patch_right_upper_j = find_first_non_zero_from_upper_right(Er_patch)
+
+    # phi and Er field on the MR patch
+    phi_sim = phi_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1]
+    Er_sim =  Er_patch[patch_left_lower_i:patch_right_upper_i+1, patch_left_lower_j:patch_right_upper_j+1]
+    r =  r1[patch_left_lower_j:patch_right_upper_j+1]
+
+    B = 1.0/np.log(0.1/0.5)
+    A = -B*np.log(0.5)
+
+    # outside EB and last cutcell
+    rmin = np.min(np.argwhere(r >= (0.1+dr)))
+    rmax = -1
+    r = r[rmin:rmax]
+    phi_sim = phi_sim[:,rmin:rmax]
+    Er_sim = Er_sim[:,rmin:rmax]
+
+    phi_theory = A + B*np.log(r)
+    phi_theory = np.tile(phi_theory, (phi_sim.shape[0],1))
+    phi_error = np.max(np.abs(phi_theory-phi_sim) / np.abs(phi_theory))
+
+    Er_theory = -B/r
+    Er_theory = np.tile(Er_theory, (Er_sim.shape[0],1))
+    Er_error = np.max(np.abs(Er_theory-Er_sim) / np.abs(Er_theory))
+
+    print(f'max error of phi[lev={level}]: {phi_error}')
+    print(f'max error of  Er[lev={level}]: {Er_error}')
+    assert(phi_error < tolerance)
+    assert(Er_error < tolerance)
+
+ts = OpenPMDTimeSeries(fn)
+level_fields = [field for field in ts.avail_fields if 'lvl' in field]
+nlevels = 0 if level_fields == [] else int(level_fields[-1][-1])
+for level in range(nlevels+1):
+    get_error_per_lev(ts,level)
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn, output_format="openpmd")

Examples/Tests/electrostatic_sphere_eb/inputs_rz_mr

@@ -30,3 +30,4 @@ diagnostics.diags_names = diag1
 diag1.intervals = 1
 diag1.diag_type = Full
 diag1.fields_to_plot = Er phi
+diag1.format = openpmd

Python/pywarpx/Bucket.py

@@ -59,8 +59,6 @@ def attrlist(self):
         for attr, value in self.argvattrs.items():
             if value is None:
                 continue
-            # --- repr is applied to value so that for floats, all of the digits are included.
-            # --- The strip of "'" is then needed when value is a string.
             if isinstance(value, str):
                 if value.find('=') > -1:
                     # --- Expressions with temporary variables need to be inside quotes
@@ -73,11 +71,11 @@ def attrlist(self):
                     continue
                 # --- For lists, tuples, and arrays make a space delimited string of the values.
                 # --- The lambda is needed in case this is a list of strings.
-                rhs = ' '.join(map(lambda s : repr(s).strip("'"), value))
+                rhs = ' '.join(map(lambda s : f'{s}', value))
             elif isinstance(value, bool):
                 rhs = 1 if value else 0
             else:
                 rhs = value
-            attrstring = '{0}.{1} = {2}'.format(self.instancename, attr, repr(rhs).strip("'"))
+            attrstring = f'{self.instancename}.{attr} = {rhs}'
             result += [attrstring]
         return result

Python/pywarpx/picmi.py

@@ -1352,6 +1352,9 @@ def init(self, kw):
 
     def solver_initialize_inputs(self):
 
+        # Open BC means FieldBoundaryType::Open for electrostatic sims, rather than perfectly-matched layer
+        BC_map['open'] = 'open'
+
         self.grid.grid_initialize_inputs()
 
         if self.relativistic:
@@ -1372,6 +1375,8 @@ def solver_initialize_inputs(self):
             pywarpx.boundary.potential_hi_y = self.grid.potential_ymax
             pywarpx.boundary.potential_hi_z = self.grid.potential_zmax
 
+        pywarpx.warpx.poisson_solver = self.method
+
 
 class GaussianLaser(picmistandard.PICMI_GaussianLaser):
     def laser_initialize_inputs(self):