Merge branch 'development' into picmi_interface_for_implicit

.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 0c3273f5e591815909180f8ffaf5b793cabbf9bc && 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

Docs/source/usage/python.rst

@@ -53,6 +53,10 @@ Field solvers define the updates of electric and magnetic fields.
 
 .. autoclass:: pywarpx.picmi.ElectrostaticSolver
 
+Object that allows smoothing of fields.
+
+.. autoclass:: pywarpx.picmi.BinomialSmoother
+
 Evolve Schemes
 --------------
 

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/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