issue with psatd solver with pml BC and domain decomposition.

[gardner-chris]

Hello WarpX community. I wasn't sure whether this is a bug or if I'm doing something wrong. I'm trying to use the moving window with psatd solver and pml boundary conditions in the radial direction for some LWFA simulations in RZ.

As a test, I ran simulations with only a laser and no particles. If I use the domain decomposition, by setting amr.max_grid_size to be a 1/4 of the longitudinal cells for example, I get non-physical fields at the divisions between grids. The issue becomes less severe as I turn up the 'psatd.noz' parameter to higher values but doesn't go away.

So my question is, is this expected behavior? The documentation for fbpic says that a stencil order of 32 is usually a good tradeoff so it seems odd to me that I would need to turn it up to much larger values when I'm only simulating a laser in free space. I would appreciate any advice for making these RZ psatd simulations work better.

Expected behavior

I would expect all of these to look like the output that I got with the 'open' boundary condition.

How to reproduce

Here is the input deck that I used.

#################################
####### GENERAL PARAMETERS ######
#################################
max_step = 2250
amr.n_cell = 256 1024
amr.max_grid_size = 256
amr.blocking_factor = 256
geometry.dims = RZ
geometry.prob_lo = 0. -70.e-6
geometry.prob_hi = 40.e-6 10.e-6
amr.max_level = 0

warpx.n_rz_azimuthal_modes = 2

algo.maxwell_solver = psatd
#psatd.nx_guard = 32
#psatd.ny_guard = 32
#psatd.nz_guard = 256

#psatd.nox = 32
#psatd.noy = 32
#psatd.noz = 128

boundary.field_lo = none damped
boundary.field_hi = pml damped

warpx.pml_ncell = 32
warpx.do_pml_in_domain = 0

#################################
############ NUMERICS ###########
#################################
warpx.verbose = 1
warpx.do_dive_cleaning = 0
warpx.use_filter = 1
warpx.filter_npass_each_dir = 0 1
warpx.cfl = 0.99
warpx.do_moving_window = 1
warpx.moving_window_dir = z
warpx.moving_window_v = 1.0
warpx.start_moving_window_step = 728
algo.load_balance_intervals = -1
algo.particle_shape = 3

#################################
############ LASER ##############
#################################
lasers.names = laser0

laser0.profile = Gaussian

###if using gaussian fill out here
laser0.profile_duration = 2.97263e-14
laser0.profile_waist = 10e-6
laser0.profile_t_peak = 175.0e-15
laser0.profile_focal_distance = 225.0e-6

laser0.position = 0. 0. -25e-6
laser0.direction = 0. 0. 1.
laser0.polarization = 1. 0. 0.
laser0.wavelength = 0.8e-6
laser0.e_max = 1.13467e13

# Diagnostics
diagnostics.diags_names = diag
diag.diag_type = Full
diag.fields_to_plot = Er Et Ez 
diag.format = openpmd
diag.openpmd_backend = h5
diag.file_prefix = diag
diag.intervals = 750::250
diag.dump_last_timestep = 1

I figured it would also be good to include the top of the output file since it has more information:

Initializing AMReX (25.03)...
MPI initialized with 1 MPI processes
MPI initialized with thread support level 3
Initializing CUDA...
CUDA initialized with 1 device.
AMReX (25.03) initialized
PICSAR (25.01)
WarpX (25.03-29-g343dfd6c3596-dirty)

    __        __             __  __
    \ \      / /_ _ _ __ _ __\ \/ /
     \ \ /\ / / _` | '__| '_ \\  /
      \ V  V / (_| | |  | |_) /  \
       \_/\_/ \__,_|_|  | .__/_/\_\
                        |_|

Level 0: dt = 2.579909799e-16 ; dx = 1.5625e-07 ; dz = 7.8125e-08

Grids Summary:
  Level 0   4 grids  262144 cells  100 % of domain
            smallest grid: 256 x 256  biggest grid: 256 x 256

-------------------------------------------------------------------------------
--------------------------- MAIN EM PIC PARAMETERS ----------------------------
-------------------------------------------------------------------------------
Precision:            | DOUBLE
Particle precision:   | DOUBLE
Geometry:             | 2D (RZ)
                      | - n_rz_azimuthal_modes = 2
Operation mode:       | Electromagnetic
                      | - vacuum
-------------------------------------------------------------------------------
Current Deposition:   | direct
Particle Pusher:      | Boris
Charge Deposition:    | standard
Field Gathering:      | energy-conserving
Particle Shape Factor:| 3
-------------------------------------------------------------------------------
Maxwell Solver:       | PSATD
                      | - update with rho is ON
                      | - current correction is ON
                      | - collocated grid
Guard cells           | - ng_alloc_EB = (32,128)
 (allocated for E/B)  |
-------------------------------------------------------------------------------
Moving window:        |    ON
                      |  - moving_window_dir = z
                      |  - moving_window_v = 299792458
-------------------------------------------------------------------------------

System information

Please check all relevant boxes and provide details.

• Operating system (name and version):
  • Linux: Rocky 8.10 Green Obsidian
  • macOS: e.g., macOS Monterey 12.4
  • Windows: e.g., Windows 11 Pro
• Version of WarpX: 25.03
• Installation method:
  • Conda
  • Spack
  • PyPI
  • Brew
  • From source with CMake
  • Module system on an HPC cluster
• Other dependencies: I don't think so.
• Computational resources:
  • MPI: e.g., 2 MPI processes
  • OpenMP: e.g., 2 OpenMP threads
  • CPU: e.g., 2 CPUs
  • GPU: NVIDIA V100

I used the HPC3 cluster at UCI, and installed using the instructions on the WarpX documentation.

Steps taken so far

To test I ran simulations scanning psatd.noz values(with default psatd.nz_guard) and psatd.nz_guard values(with default psatd.noz) with the input deck below. I also ran a test of the 'open' boundary condition and didn't have any issues which makes me think that something is going wrong with the pml boundary.

In the image below I have the results of these tests. Each row is for the input parameter in the y-axis of each column, except for the last row which is for the open boundary. I extended the psatd.noz = 64 and 128 simulations to a longer time since it took longer for the issue to appear. It looks like increasing the noz value makes the problem less severe but the nz_guard cells don't make much difference. Let me know if there's anything else you think I should test.

[EZoni]

Based on the fact that it seems like you have scanned over the number of guard cells quite thoroughly with no evident improvements, this makes me think of a number of possible issues:

• It might be that the PSATD equations in the PML domain in RZ geometry, implemented in Implement PML for the outer RZ boundary with PSATD #2211, are only an "extrapolation" of the PSATD equations in the regular domain, instead of being derived rigorously. We had a similar situation in the past for the PSATD equations in the PML in Cartesian geometry. Those were fixed in Fix PSATD equations used with PML #1513, based on the analytical derivation outlined in the notebook added in Add SymPy notebook to derive PSATD equations in PML #4316.
• Around the same time, we also noticed that using divergence cleaning for the PSATD equations in the PML domain, through the parameters warpx.do_pml_dive_cleaning and warpx.do_pml_divb_cleaning, together with divergence cleaning in the regular domain, seemed important to mitigate edge instabilities in certain 3D simulations. The divergence cleaning option for the PSATD equations in the PML domain was implemented in Divergence cleaning for PSATD in PML #1600, based again on the analytical derivation outlined in the notebook added in Add SymPy notebook to derive PSATD equations in PML #4316. However, I believe that divergence cleaning for the PSATD equations in the PML domain was never implemented in RZ geometry, after Implement PML for the outer RZ boundary with PSATD #2211.

@RemiLehe @dpgrote
Do you think these above may possibly explain the issue reported here? If so, I'm not sure if we can provide a timeline for a fix. Should open BCs in the radial direction (so boundary.field_hi = none damped instead of boundary.field_hi = pml damped) be the recommended workaround in the meantime?

[RemiLehe]

This seems similar to this FBPIC issue:
fbpic/fbpic#534
However, the above FBPIC issue is with particles.

A few things that we discussed with @EZoni @dpgrote and that could be tried:

• Check whether the above issue also exists in FBPIC (by running the exact same case). If they both have this issue, we should dive into the equations, and maybe try to implement divergence cleaning (which is currently not implemented in the RZ PML)
• Checking whether something similar to Set guard cells in corners with mixed boundary conditions #3842 would fix the issue.

[RemiLehe]

@gardner-chris Would you be able to re-run the above test with

amr.max_grid_size = 1024
amr.blocking_factor = 1024

instead of

amr.max_grid_size = 256
amr.blocking_factor = 256

and let us know whether it fixed the issue?

(This change would avoid having multiple boxes that need to communicate guard cells.)

[gardner-chris]

Hello, thanks for the responses. I ran the suggested test with pml BC's, default guard cells but with the amr.max_grid_size and amr.blocking_factor set to have only one grid and there are no issues. The bug only appears with domain decomposition.