Merge ECP-WarpX:development into EZoni:remove_tmp_particle_data

.github/workflows/cuda.yml

@@ -127,7 +127,7 @@ jobs:
         which nvcc || echo "nvcc not in PATH!"
 
         git clone https://github.com/AMReX-Codes/amrex.git ../amrex
-        cd ../amrex && git checkout --detach 69f1ac884c6aba4d9881260819ade3bb25ed8aad && cd -
+        cd ../amrex && git checkout --detach 78bdf0faabc4101d5333ebb421e553efcc7ec04e && 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

.pre-commit-config.yaml

@@ -69,7 +69,7 @@ repos:
 # Python: Ruff linter & formatter
 #   https://docs.astral.sh/ruff/
 - repo: https://github.com/astral-sh/ruff-pre-commit
-  rev: v0.9.3
+  rev: v0.9.4
   hooks:
     # Run the linter
     - id: ruff

CMakeLists.txt

@@ -1,7 +1,7 @@
 # Preamble ####################################################################
 #
 cmake_minimum_required(VERSION 3.24.0)
-project(WarpX VERSION 25.01)
+project(WarpX VERSION 25.02)
 
 include(${WarpX_SOURCE_DIR}/cmake/WarpXFunctions.cmake)
 

Docs/requirements.txt

@@ -27,5 +27,6 @@ sphinx-copybutton
 sphinx-design
 sphinx_rtd_theme>=1.1.1
 sphinxcontrib-bibtex
+sphinxcontrib-googleanalytics
 sphinxcontrib-napoleon
 yt  # for checksumAPI

Docs/source/conf.py

@@ -56,8 +56,13 @@
     "sphinx_design",
     "breathe",
     "sphinxcontrib.bibtex",
+    "sphinxcontrib.googleanalytics",
 ]
 
+# Google Analytics
+googleanalytics_id = "G-QZGY5060MZ"
+googleanalytics_enabled = True
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 
@@ -107,9 +112,9 @@ def __init__(self, *args, **kwargs):
 # built documents.
 #
 # The short X.Y version.
-version = "25.01"
+version = "25.02"
 # The full version, including alpha/beta/rc tags.
-release = "25.01"
+release = "25.02"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

Docs/source/refs.bib

@@ -35,6 +35,17 @@ @ARTICLE{Birdsall1991
 year = {1991}
 }
 
+@misc{Janssen2016
+author = {Janssen, J. F. J. and Pitchford L. C. and Hagelaar G. J. M. and van Dijk J.},
+doi = {10.1088/0963-0252/25/5/055026},
+journal = {Plasma Sources Science and Technology},
+number = {5},
+pages = {055026},
+title = {{Evaluation of angular scattering models for electron-neutral collisions in Monte Carlo simulations}},
+volume = {25},
+year = {2016}
+}
+
 @misc{Lim2007,
 author = {Lim, Chul-Hyun},
 issn = {0419-4217},

Docs/source/theory/boundary_conditions.rst

@@ -301,3 +301,23 @@ the right boundary is reflecting.
 
 .. bibliography::
    :keyprefix: bc-
+
+.. _theory-bc-pmc:
+
+Perfect Magnetic Conductor
+----------------------------
+
+This boundary can be used to model a symmetric surface, where charges and current are
+symmetric across the boundary.
+This is equivalent to the Neumann (zero-derivative) boundary condition.
+For the electromagnetic solve, at PMC, the tangential magnetic field and the normal electric
+field are odd across the boundary and set to 0 on the boundary.
+In the guard-cell region, those fields are set equal and
+opposite to the respective field component in the mirror location across the PMC boundary.
+The other components, the normal magnetic field and tangential electric field, are even
+and set equal to the field component in the mirror location in the domain across the PMC boundary.
+
+The PMC boundary condition also impacts the deposition of charge and current density.
+The charge and current densities deposited into the guard cells are reflected back into
+the domain, adding them to the mirror cells in the domain.
+This represents the charge and current from the virtual symmetric particles in the guard cells.

Docs/source/theory/multiphysics/collisions.rst

@@ -121,13 +121,23 @@ The particle velocity in the COM frame is then isotropically scattered using the
 Back scattering
 ^^^^^^^^^^^^^^^
 
-The process is the same as for elastic scattering above expect the scattering angle is fixed at :math:`\pi`, meaning the particle velocity in the COM frame is updated to :math:`-\vec{u}_c`.
+The process is the same as for elastic scattering above except the scattering angle is fixed at :math:`\pi`, meaning the particle velocity in the COM frame is updated to :math:`-\vec{u}_c`.
 
 Excitation
 ^^^^^^^^^^
 
 The process is also the same as for elastic scattering except the excitation energy cost is subtracted from the particle energy. This is done by reducing the velocity before a scattering angle is chosen.
 
+Forward scattering
+^^^^^^^^^^^^^^^^^^
+
+This process operates in two ways:
+
+1. If an excitation energy cost is provided, the energy cost is subtracted from the particle energy and no scattering is performed.
+2. If an excitation energy cost is not provided, the particle is not scattered and the velocity is unchanged (corresponding to a scattering angle of :math:`0` in the elastic scattering process above).
+
+See :cite:t:`b-Janssen2016` for a recommended use of this process.
+
 Benchmarks
 ----------
 

Docs/source/usage/parameters.rst

@@ -533,6 +533,8 @@ Domain Boundary Conditions
 
     * ``pec``: This option can be used to set a Perfect Electric Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pec>` for more details. Note that PEC boundary is invalid at `r=0` for the RZ solver. Please use ``none`` option. This boundary condition does not work with the spectral solver.
 
+    * ``pmc``: This option can be used to set a Perfect Magnetic Conductor at the simulation boundary. Please see the :ref:`PEC theory section <theory-bc-pmc>` for more details. This is equivalent to ``Neumann``. This boundary condition does not work with the spectral solver.
+
     * ``pec_insulator``: This option specifies a mixed perfect electric conductor and insulator boundary, where some part of the
       boundary is PEC and some is insulator. In the insulator portion, the normal fields are extrapolated and the tangential fields
       are either set to the specified value or extrapolated. The region that is insulator is specified using a spatially dependent expression with the insulator being in the area where the value of the expression is greater than zero.
@@ -2167,8 +2169,8 @@ Details about the collision models can be found in the :ref:`theory section <mul
 
 * ``<collision_name>.scattering_processes`` (`strings` separated by spaces)
     Only for ``dsmc`` and ``background_mcc``. The scattering processes that should be
-    included. Available options are ``elastic``, ``back`` & ``charge_exchange``
-    for ions and ``elastic``, ``excitationX`` & ``ionization`` for electrons.
+    included. Available options are ``elastic``, ``excitationX``, ``forward``, ``back``, and ``charge_exchange``
+    for ions and ``elastic``, ``excitationX``, ``ionization`` & ``forward`` for electrons.
     Multiple excitation events can be included for electrons corresponding to
     excitation to different levels, the ``X`` above can be changed to a unique
     identifier for each excitation process. For each scattering process specified
@@ -3182,6 +3184,12 @@ This shifts analysis from post-processing to runtime calculation of reduction op
         Note that the fields are averaged on the cell centers before their maximum values are
         computed.
 
+    * ``FieldPoyntingFlux``
+        Integrates the normal Poynting flux over each domain boundary surface and also integrates the flux over time.
+        This provides the power and total energy loss into or out of the simulation domain.
+        The output columns are the flux for each dimension on the lower boundaries, then the higher boundaries,
+        then the integrated energy loss for each dimension on the the lower and higher boundaries.
+
     * ``FieldProbe``
         This type computes the value of each component of the electric and magnetic fields
         and of the Poynting vector (a measure of electromagnetic flux) at points in the domain.

Examples/Tests/ohm_solver_em_modes/analysis_rz.py

@@ -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([61.02377286, 19.80026021, 100.47687017, 10.83331295])
+        amps, np.array([59.23850009, 19.26746169, 92.65794174, 10.83627164])
     )

Examples/Tests/pec/CMakeLists.txt

@@ -40,3 +40,33 @@ add_warpx_test(
     "analysis_default_regression.py --path diags/diag1000010"  # checksum
     OFF  # dependency
 )
+
+add_warpx_test(
+    test_3d_pmc_field  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_pmc_field  # inputs
+    "analysis_pec.py diags/diag1000134"  # analysis
+    "analysis_default_regression.py --path diags/diag1000134"  # checksum
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_2d_pec_field_insulator_implicit  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_pec_field_insulator_implicit  # inputs
+    "analysis_pec_insulator_implicit.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_2d_pec_field_insulator_implicit_restart  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_pec_field_insulator_implicit_restart  # inputs
+    "analysis_pec_insulator_implicit.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
+    test_2d_pec_field_insulator_implicit  # dependency
+)

Examples/Tests/pec/analysis_pec_insulator_implicit.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from
+# the scripts `inputs_test_2d_pec_field_insulator_implicit` and
+# `inputs_test_2d_pec_field_insulator_implicit_restart`.
+# The scripts model an insulator boundary condition on part of the
+# upper x boundary that pushes B field into the domain. The implicit
+# solver is used, converging to machine tolerance. The energy accounting
+# should be exact to machine precision, so that the energy is the system
+# should be the same as the amount of energy pushed in from the boundary.
+# This is checked using the FieldEnergy and FieldPoyntingFlux reduced
+# diagnostics.
+import sys
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+# this will be the name of the plot file
+fn = sys.argv[1]
+
+EE = np.loadtxt(f"{fn}/../reducedfiles/fieldenergy.txt", skiprows=1)
+SS = np.loadtxt(f"{fn}/../reducedfiles/poyntingflux.txt", skiprows=1)
+SSsum = SS[:, 2:6].sum(1)
+EEloss = SS[:, 7:].sum(1)
+
+dt = EE[1, 1]
+
+fig, ax = plt.subplots()
+ax.plot(EE[:, 0], EE[:, 2], label="field energy")
+ax.plot(SS[:, 0], -EEloss, label="-flux*dt")
+ax.legend()
+ax.set_xlabel("time (s)")
+ax.set_ylabel("energy (J)")
+fig.savefig("energy_history.png")
+
+fig, ax = plt.subplots()
+ax.plot(EE[:, 0], (EE[:, 2] + EEloss) / EE[:, 2].max())
+ax.set_xlabel("time (s)")
+ax.set_ylabel("energy difference/max energy (1)")
+fig.savefig("energy_difference.png")
+
+tolerance_rel = 1.0e-13
+
+energy_difference_fraction = np.abs((EE[:, 2] + EEloss) / EE[:, 2].max()).max()
+print(f"energy accounting error = {energy_difference_fraction}")
+print(f"tolerance_rel = {tolerance_rel}")
+
+assert energy_difference_fraction < tolerance_rel

Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit

@@ -0,0 +1,73 @@
+# Maximum number of time steps
+max_step = 20
+
+# number of grid points
+amr.n_cell = 32 32
+amr.blocking_factor = 16
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo = 0.     2.e-2 # physical domain
+geometry.prob_hi = 1.e-2  3.e-2
+
+# Boundary condition
+boundary.field_lo = neumann      periodic
+boundary.field_hi = pec_insulator periodic
+
+insulator.area_x_hi(y,z) = (2.25e-2 <= z and z <= 2.75e-2)
+insulator.By_x_hi(y,z,t) = min(t/1.0e-12,1)*1.e1*3.3e-4
+
+warpx.serialize_initial_conditions = 1
+
+# Implicit setup
+# Note that this is the CFL step size for the explicit simulation, over 2.
+# This value allows quick convergence of the Picard solver.
+warpx.const_dt = 7.37079480234276e-13/2.
+
+algo.maxwell_solver = Yee
+algo.evolve_scheme = "theta_implicit_em"
+#algo.evolve_scheme = "semi_implicit_em"
+
+implicit_evolve.theta = 0.5
+#implicit_evolve.max_particle_iterations = 21
+#implicit_evolve.particle_tolerance = 1.0e-12
+
+implicit_evolve.nonlinear_solver = "picard"
+picard.verbose = true
+picard.max_iterations = 25
+picard.relative_tolerance = 0.0
+picard.absolute_tolerance = 0.0
+picard.require_convergence = false
+
+#implicit_evolve.nonlinear_solver = "newton"
+#newton.verbose = true
+#newton.max_iterations = 20
+#newton.relative_tolerance = 1.0e-20
+#newton.absolute_tolerance = 0.0
+#newton.require_convergence = false
+
+#gmres.verbose_int = 2
+#gmres.max_iterations = 1000
+#gmres.relative_tolerance = 1.0e-20
+#gmres.absolute_tolerance = 0.0
+
+# Verbosity
+warpx.verbose = 1
+
+# Diagnostics
+diagnostics.diags_names = diag1 chk
+diag1.intervals = 20
+diag1.diag_type = Full
+
+chk.intervals = 10
+chk.diag_type = Full
+chk.format = checkpoint
+
+warpx.reduced_diags_names = fieldenergy poyntingflux
+poyntingflux.type = FieldPoyntingFlux
+poyntingflux.intervals = 1
+fieldenergy.type = FieldEnergy
+fieldenergy.intervals = 1

Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart

@@ -0,0 +1,5 @@
+# base input parameters
+FILE = inputs_test_2d_pec_field_insulator_implicit
+
+# test input parameters
+amr.restart = "../test_2d_pec_field_insulator_implicit/diags/chk000010"

Examples/Tests/pec/inputs_test_3d_pmc_field

@@ -0,0 +1,54 @@
+# Set-up to test the PMC Boundary condition for the fields
+# Constructive interference between the incident and reflected wave result in a
+# standing wave.
+
+# max step
+max_step = 134
+
+# number of grid points
+amr.n_cell = 32 32 256
+
+# Maximum allowable size of each subdomain
+amr.max_grid_size = 1024
+amr.blocking_factor = 32
+
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 3
+geometry.prob_lo     =  -8.e-6 -8.e-6 -4.e-6
+geometry.prob_hi     =   8.e-6  8.e-6  4.e-6
+
+# Boundary condition
+boundary.field_lo = periodic periodic pmc
+boundary.field_hi = periodic periodic pmc
+
+warpx.serialize_initial_conditions = 1
+
+# Verbosity
+warpx.verbose = 1
+
+# Algorithms
+algo.current_deposition = esirkepov
+# CFL
+warpx.cfl = 0.9
+
+
+my_constants.z1 = -2.e-6
+my_constants.z2 = 2.e-6
+my_constants.wavelength = 1.e-6
+warpx.E_ext_grid_init_style = parse_E_ext_grid_function
+warpx.Ez_external_grid_function(x,y,z) = "0."
+warpx.Ex_external_grid_function(x,y,z) = "0."
+warpx.Ey_external_grid_function(x,y,z) = "((1.e5*sin(2*pi*(z)/wavelength)) * (z<z2) * (z>z1))"
+
+warpx.B_ext_grid_init_style = parse_B_ext_grid_function
+warpx.Bx_external_grid_function(x,y,z)= "(((-1.e5*sin(2*pi*(z)/wavelength))/clight))*(z<z2) * (z>z1) "
+warpx.By_external_grid_function(x,y,z)= "0."
+warpx.Bz_external_grid_function(x,y,z) = "0."
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 134
+diag1.diag_type = Full
+diag1.fields_to_plot = Ey Bx