Update memaster with latest upstream/development (#169)

* Moving Frame Field Probe Functionality (ECP-WarpX#2996)

* DRAFT for Moving Frame

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Adds m_last_compute_step

* Generalized SetParticlePosition, semicolons and formatting

* Update Source/Diagnostics/ReducedDiags/FieldProbe.H

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Moved particle push before calculations. Condensed stuff.

* Fixed Velocity, added start and stop time functionality

* Documentation

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Empty-Commit

* Fixed ParallelFor from review

* FP moving window in laser_acceleration 1d, 2d, 3d test

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Line detector instead of point

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Refactoring: More General Interfaces (ECP-WarpX#2991)

* Refactoring: More General Interfaces (J, rho)

* Refactoring: More General Interfaces (E, B)

* Refactoring: More General Interfaces (rho)

* Fix const Correctness

* Fix const Correctness (continued)

* Fix Bugs

* Refactoring: More General Interfaces (averaged E, averaged B)

* LoadBalanceCosts: num_cells & num_macro_particles (ECP-WarpX#3019)

* Examples: LBC in Laser-Ion Example

* LoadBalanceCosts: num_cells & num_macro_particles

* Fix Script: Wrong Header Counting

Spliced out the first to columns twice.

* Fix Script: Blocking Factor

The script did not work if a direction only had one block

* Fix Script: 2D/3D

There is no generic way from the data we write to find this out,
without breaking corner cases such as one-block in one direction.

* LoadBalanceCosts: Fix Unused Vars in Script

Leftover to ECP-WarpX#3019

* Semicoarsening in electrostatic solvers (ECP-WarpX#3024)

* Semicoarsening in electrostatic solvers

Enable semicoarsening in electrostatic solvers if the problem is anisotropic
due to relativistic beam and/or anisotropic cell spacing.

* Update relativistic_space_charge_initialization benchmark

* AMReX: Update to latest commit (ECP-WarpX#3021)

* Added PlasmaLens class to PICMI (ECP-WarpX#3025)

* Added PlasmaLens class to PICMI

* Added CI test

* Added PICMI input file

* Fixed the output dir for the CI test

* Add `_plt` to Output File Prefix

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Minor fixes removing warnings from MSVC C++ (ECP-WarpX#3022)

* Fix MPI Signal Handling on Perlmutter (ECP-WarpX#3029)

* Fix: `MPI_CXX_BOOL` -> `MPI_BYTE`

C99 types were aded in MPI-2.2, while Cray's MPICH fork in version
8.1.13 defines `MPI_CXX_BOOL` to `MPI_DATATYPE_NULL` on Perlmutter.

We could use `MPI_C_BOOL`, which is technically a `_Bool` from
[<cstdbool>](https://en.cppreference.com/w/cpp/header/cstdbool)
(deprecated: C++17; removed: C++20) - or we simply do a static
assert on `sizeof(bool)` and communicate as a `MPI_BYTE` or
`MPI_CHAR`.

* Signals: Do no MPI Comms if none is configured

* Docs: Link and Explain Allowed Signal Values

* Review Comments & Style

* Extend `do_pml_Lo/Hi` to MR Levels (ECP-WarpX#2890)

* Vay Deposition: Fix SyncCurrent, FFTs (ECP-WarpX#3012)

* Implemented the BCK fallback for the ECT solver (ECP-WarpX#3016)

* Implemented the BCK fallback for the ECT solver

* Ignored an unused variable

* Bug fix

* Fix issue with preprocessor directive

* Added comments

* Implemented suggestions from code review

* Fixed tilebox

* Improved warning

* Using WarpX::RecordWarning instead of amrex::Print()

* Ignoring RZ

* Improving comment

* Bug Fix

* Lowered warning priority

* Ignoring low priority warnings in ECT tests

Co-authored-by: lgiacome <lorenzo.giacome@cern.ch>

* No longer use a while loop in driving the simulation (#167)

* add persistent flag for whether a terminate signal was received that can be read by python

* handle checkpointing through sim_control

* change terminateFlag from type int to bool

* fix unintended changes to changelog

* specify terminateFlag return type

* code cleanup

* revert upstream changes

* remove use of while loop to run simulation

* changes requested during PR review

Co-authored-by: Tiberius Rheaume <35204125+TiberiusRheaume@users.noreply.github.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>
Co-authored-by: David Grote <grote1@llnl.gov>
Co-authored-by: Lorenzo Giacomel <47607756+lgiacome@users.noreply.github.com>
Co-authored-by: lgiacome <lorenzo.giacome@cern.ch>

.github/workflows/cuda.yml

@@ -106,7 +106,7 @@ jobs:
         which nvcc || echo "nvcc not in PATH!"
 
         git clone https://github.com/AMReX-Codes/amrex.git ../amrex
-        cd amrex && git checkout --detach 22.04 && cd -
+        cd amrex && git checkout --detach ec648a92a89ffc907a87d4b02ba4ddb8a2a932a5 && cd -
         make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 2
 
   build_nvhpc21-11-nvcc:

Docs/source/glossary.rst

@@ -15,6 +15,7 @@ Abbreviations
 * **BLAST:** `Beam, Plasma & Accelerator Simulation Toolkit <https://blast.lbl.gov>`__
 * **AMR:** adaptive mesh-refinement
 * **BC:** boundary condition (of a simulation)
+* **BCK:** `Benkler-Chavannes-Kuster <https://ieeexplore.ieee.org/document/1638381>`__ method, a stabilization technique for small cells in the electromagnetic solver
 * **BTD:** backtransformed diagnosics, a method to collect data for analysis from a *boosted frame* simulation
 * **CFL:** the Courant-Friedrichs-Lewy condition, a numerical parameter for the numerical convergence of PDE solvers
 * **CI:** continuous integration, automated tests that we perform before a proposed code-change is accepted; see PR

Docs/source/usage/parameters.rst

@@ -64,15 +64,6 @@ Overall simulation parameters
     It is mainly intended for debug purposes, and is best used with
     ``warpx.always_warn_immediately=1``.
 
-* ``warpx.break_signals`` (array of `string`, separated by spaces) optional
-    A list of signal names or numbers that the simulation should
-    handle by cleanly terminating at the next timestep
-
-* ``warpx.checkpoint_signals`` (array of `string`, separated by spaces) optional
-    A list of signal names or numbers that the simulation should
-    handle by outputting a checkpoint at the next timestep. A
-    diagnostic of type `checkpoint` must be configured.
-
 * ``warpx.random_seed`` (`string` or `int` > 0) optional
     If provided ``warpx.random_seed = random``, the random seed will be determined
     using `std::random_device` and `std::clock()`,
@@ -148,6 +139,37 @@ Overall simulation parameters
     Note that even with this set to ``1`` WarpX will not catch all out-of-memory events yet when operating close to maximum device memory.
     `Please also see the documentation in AMReX <https://amrex-codes.github.io/amrex/docs_html/GPU.html#inputs-parameters>`_.
 
+Signal Handling
+^^^^^^^^^^^^^^^
+
+WarpX can handle Unix (Linux/macOS) `process signals <https://en.wikipedia.org/wiki/Signal_(IPC)>`__.
+This can be useful to configure jobs on HPC and cloud systems to shut down cleanly when they are close to reaching their allocated walltime or to steer the simulation behavior interactively.
+
+Allowed signal names are documented in the `C++ standard <https://en.cppreference.com/w/cpp/utility/program/SIG_types>`__ and `POSIX <https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/signal.h.html>`__.
+We follow the same naming, but remove the ``SIG`` prefix, e.g., the WarpX signal configuration name for ``SIGINT`` is ``INT``.
+
+* ``warpx.break_signals`` (array of `string`, separated by spaces) optional
+    A list of signal names or numbers that the simulation should
+    handle by cleanly terminating at the next timestep
+
+* ``warpx.checkpoint_signals`` (array of `string`, separated by spaces) optional
+    A list of signal names or numbers that the simulation should
+    handle by outputting a checkpoint at the next timestep. A
+    diagnostic of type `checkpoint` must be configured.
+
+.. note::
+
+   Certain signals are only available on specific platforms, please see the links above for details.
+   Typically supported on Linux and macOS are ``HUP``, ``INT``, ``QUIT``, ``ABRT``, ``USR1``, ``USR2``, ``TERM``, ``TSTP``, ``URG``, and ``IO`` among others.
+
+   Signals to think about twice before overwriting in *interactive simulations*:
+   Note that ``INT`` (interupt) is the signal that ``Ctrl+C`` sends on the terminal, which most people use to abort a process; once overwritten you need to abort interactive jobs with, e.g., ``Ctrl+\`` (``QUIT``) or sending the ``KILL`` signal.
+   The ``TSTP`` (terminal stop) command is sent interactively from ``Ctrl+Z`` to temporarily send a process to sleep (until send in the background with commands such as ``bg`` or continued with ``fg``), overwriting it would thus disable that functionality.
+   The signals ``KILL`` and ``STOP`` cannot be used.
+
+   The ``FPE`` signal should not be overwritten in WarpX, as it is `controlled by AMReX <https://amrex-codes.github.io/amrex/docs_html/Debugging.html#breaking-into-debuggers>`__ for :ref:`debug workflows that catch invalid floating-point operations <debugging_warpx>`.
+
+
 .. _running-cpp-parameters-box:
 
 Setting up the field mesh
@@ -2335,6 +2357,7 @@ Reduced Diagnostics
         otherwise it is set to ``1``.
         Integrated electric and magnetic field components can instead be obtained by specifying
         ``<reduced_diags_name>.integrate == true``.
+        In a *moving window* simulation, the FieldProbe can be set to follow the moving frame by specifying ``<reduced_diags_name>.do_moving_window_FP = 1`` (default 0).
 
         .. warning::
 

Docs/source/usage/workflows/plot_distribution_mapping.rst

@@ -52,7 +52,8 @@ can be used to read the data:
     import plot_distribution_mapping as pdm
 
     sim_knapsack = pdm.SimData('LBC_knapsack.txt', # Data directory
-                               [2800]              # Files to process
+                               [2800],             # Files to process
+                               is_3D=False         # if this is a 2D sim
                               )
     sim_sfc = pdm.SimData('LBC_sfc.txt', [2800])
 

Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py

@@ -34,11 +34,11 @@
 m = 0
 n = 1
 p = 1
-Lx = 1.06
-Ly = 1.06
-Lz = 1.06
+Lx = 1
+Ly = 1
+Lz = 1
 h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2
-theta = np.pi/8
+theta = np.pi/6
 
 # Open the right plot file
 filename = sys.argv[1]

Examples/Modules/embedded_boundary_rotated_cube/inputs_3d

@@ -1,5 +1,5 @@
 stop_time = 5.303669113650618e-09
-amr.n_cell = 32 32 32
+amr.n_cell = 64 64 64
 amr.max_grid_size = 128
 amr.max_level = 0
 
@@ -13,23 +13,23 @@ boundary.field_hi = pec pec pec
 
 algo.maxwell_solver = ect
 
-my_constants.xmin = -0.53
-my_constants.ymin = -0.53
-my_constants.zmin = -0.53
-my_constants.xmax = 0.53
-my_constants.ymax = 0.53
-my_constants.zmax = 0.53
+my_constants.xmin = -0.5
+my_constants.ymin = -0.5
+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/8
+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)))"
 
 my_constants.m = 0
 my_constants.n = 1
 my_constants.p = 1
-my_constants.Lx = 1.06
-my_constants.Ly = 1.06
-my_constants.Lz = 1.06
+my_constants.Lx = 1
+my_constants.Ly = 1
+my_constants.Lz = 1
 my_constants.x_cent = 0.
 my_constants.y_cent = 0.
 my_constants.z_cent = 0.

Examples/Physics_applications/laser_acceleration/inputs_1d

@@ -70,3 +70,15 @@ diagnostics.diags_names = diag1
 diag1.intervals = 100
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+
+# Reduced Diagnostics
+warpx.reduced_diags_names               = FP
+
+FP.type = FieldProbe
+FP.intervals = 10
+FP.integrate = 0
+FP.probe_geometry = Line
+FP.z_probe = -56e-6
+FP.z1_probe = 12e-6
+FP.resolution = 100
+FP.do_moving_window_FP = 1

Examples/Physics_applications/laser_acceleration/inputs_2d

@@ -87,3 +87,17 @@ diagnostics.diags_names = diag1
 diag1.intervals = 200
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+
+# Reduced Diagnostics
+warpx.reduced_diags_names               = FP
+
+FP.type = FieldProbe
+FP.intervals = 10
+FP.integrate = 0
+FP.probe_geometry = Line
+FP.x_probe = 0
+FP.z_probe = -56e-6
+FP.x1_probe = 0
+FP.z1_probe = 12e-6
+FP.resolution = 300
+FP.do_moving_window_FP = 1

Examples/Physics_applications/laser_acceleration/inputs_3d

@@ -77,3 +77,19 @@ diagnostics.diags_names = diag1
 diag1.intervals = 100
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+
+# Reduced Diagnostics
+warpx.reduced_diags_names               = FP
+
+FP.type = FieldProbe
+FP.intervals = 10
+FP.integrate = 0
+FP.probe_geometry = Line
+FP.x_probe = 0
+FP.y_probe = 0
+FP.z_probe = -56e-6
+FP.x1_probe = 0
+FP.y1_probe = 0
+FP.z1_probe = 12e-6
+FP.resolution = 300
+FP.do_moving_window_FP = 1

Examples/Physics_applications/laser_ion/inputs

@@ -64,7 +64,7 @@ warpx.numprocs = 1 2   # 2 MPI ranks
 #   The grid & block parameters above are needed for load balancing:
 #   an average of ~10 grids per MPI rank (and device) are a good granularity
 #   to allow efficient load-balancing as the simulation evolves
-algo.load_balance_intervals = 10
+algo.load_balance_intervals = 100
 algo.load_balance_costs_update = Heuristic
 
 # particle bin-sorting on GPU (ideal defaults not investigated in 2D)
@@ -231,7 +231,7 @@ openPMDbw.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz<0)
 # histograms with 2.0 degree acceptance angle in fw direction
 # 2 deg * pi / 180 : 0.03490658503 rad
 # half-angle +/-   : 0.017453292515 rad
-warpx.reduced_diags_names                   = histuH histue histuzAll FieldProbe_Z FieldProbe_ScatPoint FieldProbe_ScatLine
+warpx.reduced_diags_names                   = histuH histue histuzAll FieldProbe_Z FieldProbe_ScatPoint FieldProbe_ScatLine LBC
 
 histuH.type                                 = ParticleHistogram
 histuH.intervals                            = 100
@@ -291,6 +291,10 @@ FieldProbe_ScatLine.x1_probe = 2.5e-6
 FieldProbe_ScatLine.z1_probe = 15e-6
 FieldProbe_ScatLine.resolution = 201
 
+# check computational load per box
+LBC.type = LoadBalanceCosts
+LBC.intervals = 100
+
 #################################
 # Physical Background
 #

Examples/Tests/Langmuir/analysis_langmuir_multi.py

@@ -149,14 +149,16 @@ def get_theoretical_field( field, t ):
 
 assert( error_rel < tolerance_rel )
 
-# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when
-# current correction (psatd.do_current_correction=1) is applied or when
-# Vay current deposition (algo.current_deposition=vay) is used
+# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E)
+# with current correction (and periodic single box option) or with Vay current deposition
+if current_correction:
+    tolerance = 1e-9
+elif vay_deposition:
+    tolerance = 1e-3
 if current_correction or vay_deposition:
     rho  = data[('boxlib','rho')].to_ndarray()
     divE = data[('boxlib','divE')].to_ndarray()
     error_rel = np.amax( np.abs( divE - rho/epsilon_0 ) ) / np.amax( np.abs( rho/epsilon_0 ) )
-    tolerance = 1.e-9
     print("Check charge conservation:")
     print("error_rel = {}".format(error_rel))
     print("tolerance = {}".format(tolerance))

Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py

@@ -121,14 +121,16 @@ def get_theoretical_field( field, t ):
 
 assert( error_rel < tolerance_rel )
 
-# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E) when
-# current correction (psatd.do_current_correction=1) is applied or when
-# Vay current deposition (algo.current_deposition=vay) is used
+# Check relative L-infinity spatial norm of rho/epsilon_0 - div(E)
+# with current correction (and periodic single box option) or with Vay current deposition
+if current_correction:
+    tolerance = 1e-9
+elif vay_deposition:
+    tolerance = 1e-3
 if current_correction or vay_deposition:
     rho  = data[('boxlib','rho')].to_ndarray()
     divE = data[('boxlib','divE')].to_ndarray()
     error_rel = np.amax( np.abs( divE - rho/epsilon_0 ) ) / np.amax( np.abs( rho/epsilon_0 ) )
-    tolerance = 1.e-9
     print("Check charge conservation:")
     print("error_rel = {}".format(error_rel))
     print("tolerance = {}".format(tolerance))

Examples/Tests/plasma_lens/PICMI_inputs_3d.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+
+from pywarpx import picmi
+
+# Physical constants
+c = picmi.constants.c
+q_e = picmi.constants.q_e
+
+# Number of time steps
+max_steps = 84
+
+# Number of cells
+nx = 16
+ny = 16
+nz = 16
+
+# Physical domain
+xmin = -1.
+xmax =  1.
+ymin = -1.
+ymax =  1.
+zmin =  0.
+zmax =  2.
+
+# Create grid
+grid = picmi.Cartesian3DGrid(number_of_cells = [nx, ny, nz],
+                             lower_bound = [xmin, ymin, zmin],
+                             upper_bound = [xmax, ymax, zmax],
+                             lower_boundary_conditions = ['dirichlet', 'dirichlet', 'dirichlet'],
+                             upper_boundary_conditions = ['dirichlet', 'dirichlet', 'dirichlet'],
+                             lower_boundary_conditions_particles = ['absorbing', 'absorbing', 'absorbing'],
+                             upper_boundary_conditions_particles = ['absorbing', 'absorbing', 'absorbing'])
+
+# Particles
+vel_z = 0.5*c
+multiparticles_distribution = picmi.ParticleListDistribution(x = [0.05, 0.],
+                                                             y = [0., 0.04],
+                                                             z = [0.05, 0.05],
+                                                             ux = [0., 0.],
+                                                             uy = [0., 0.],
+                                                             uz = [vel_z, vel_z],
+                                                             weight = [1., 1.])
+
+electrons = picmi.Species(particle_type = 'electron',
+                          name = 'electrons',
+                          initial_distribution = multiparticles_distribution)
+
+# Plasma lenses
+plasma_lenses = picmi.PlasmaLens(period = 0.5,
+                                 starts = [0.1, 0.11, 0.12, 0.13],
+                                 lengths = [0.1, 0.11, 0.12, 0.13],
+                                 strengths_E = [600000., 800000., 600000., 200000.],
+                                 strengths_B = [0.0, 0.0, 0.0, 0.0])
+
+# Electromagnetic solver
+solver = picmi.ElectromagneticSolver(grid = grid,
+                                     method = 'Yee',
+                                     cfl = 0.7)
+
+# Diagnostics
+part_diag1 = picmi.ParticleDiagnostic(name = 'diag1',
+                                      period = max_steps,
+                                      species = [electrons],
+                                      data_list = ['ux', 'uy', 'uz'],
+                                      write_dir = '.',
+                                      warpx_file_prefix = 'Python_plasma_lens_plt')
+
+# Set up simulation
+sim = picmi.Simulation(solver = solver,
+                       max_steps = max_steps,
+                       verbose = 1,
+                       particle_shape = 'linear',
+                       warpx_serialize_initial_conditions = 1,
+                       warpx_do_dynamic_scheduling = 0)
+
+# Add plasma electrons
+sim.add_species(electrons, layout = None)
+
+# Add the plasma lenses
+sim.add_applied_field(plasma_lenses)
+
+# Add diagnostics
+sim.add_diagnostic(part_diag1)
+
+# Write input file that can be used to run with the compiled version
+#sim.write_input_file(file_name = 'inputs_3d_picmi')
+
+# Advance simulation until last time step
+sim.step(max_steps)

Examples/Tests/plasma_lens/analysis.py

@@ -60,7 +60,11 @@ def applylens(x0, vx0, vz0, gamma, lens_length, lens_strength):
     return x1, vx1
 
 clight = c
-vel_z = eval(ds.parameters.get('my_constants.vel_z'))
+try:
+    vel_z = eval(ds.parameters.get('my_constants.vel_z'))
+except TypeError:
+    # vel_z is not saved in my_constants with the PICMI version
+    vel_z = 0.5*c
 
 plasma_lens_period = float(ds.parameters.get('particles.repeated_plasma_lens_period'))
 plasma_lens_starts = [float(x) for x in ds.parameters.get('particles.repeated_plasma_lens_starts').split()]

Examples/Tests/plasma_lens/inputs_3d

@@ -18,6 +18,7 @@ boundary.particle_hi = absorbing absorbing absorbing
 
 # Algorithms
 algo.particle_shape = 1
+warpx.cfl = 0.7
 
 my_constants.vel_z = 0.5*clight
 

Examples/Tests/plasma_lens/inputs_boosted_3d

@@ -20,6 +20,7 @@ boundary.particle_hi = absorbing absorbing absorbing
 algo.particle_shape = 1
 warpx.gamma_boost = 2.
 warpx.boost_direction = z
+warpx.cfl = 0.7
 
 my_constants.vel_z = 0.5*clight