MLMG "Bottom solve failed" + spurious transverse E-field with EB electrodes and Neumann/periodic BCs

[n01r]

Description

An electrostatic MLMG solve with an embedded-boundary (EB) Dirichlet potential and all-Neumann/periodic domain BCs prints MLMG: Bottom solve failed. and leaves an unphysical transverse field (Ey, Ez ≠ 0 where they should be exactly 0). It is triggered by domain decomposition into more than one box — not by MPI rank count, resolution, or the EB implicit-function form (all ruled out below). On GPU the failure collapses convergence and the artifact reaches ~1e-3 V/m; on CPU the same failure is benign.

This issue was first flagged by @agargone.

Expected behavior

Two planar electrodes at x = ±0.052 in a box with Neumann-x/periodic-y,z give a 1-D field: Ex(x) only, Ey = Ez = 0 everywhere. The solve should converge and produce zero transverse field regardless of how the domain is split into boxes/ranks.

How to reproduce

Inputs:

• input_repro (saved as txt in attachment)
• reproduce.sh

Single rank, no mpirun needed — only amr.max_grid_size (i.e. the number of boxes) changes:

./warpx.3d input_repro amr.max_grid_size=528   # 1 box  -> 0 failures, max|Ey| = 0
./warpx.3d input_repro amr.max_grid_size=64    # 18 boxes -> "Bottom solve failed", max|Ey| != 0

Compare grep -c "Bottom solve failed" <log> and column max_Ey in diags/reducedfiles/fmax.txt. Key input settings: algo.maxwell_solver = Yee, boundary.field_lo/hi = neumann periodic periodic, warpx.eb_potential(x,y,z,t) = "-100000.0*(x<(-0.052+5.0e-3))", warpx.self_fields_verbosity = 2.

System information

• Operating system: [x] Linux — Ubuntu on Windows (WSL2), kernel 6.6.87.2 (originally observed on Perlmutter, SUSE SLES)
• Version of WarpX: 26.08-47-g9e32d6b19 (AMReX 26.08-26-g2cf4fbcde02a)
• Installation method: [x] From source with CMake
• Other dependencies: none required for the reproducer (EB on; no OpenPMD/QED/FFT needed)
• Computational resources:
  • [x] MPI — bug appears with any decomposition giving >1 box (even 1 rank + small max_grid_size)
  • [x] OpenMP — CPU build
  • [x] CPU — reproduces the warning (benign artifact)
  • [x] GPU — NVIDIA (Perlmutter A100, and locally RTX 3500 Ada); the severe ~1e-3 artifact

Steps taken so far

Isolated the trigger to multi-box decomposition (single box → 0 failures / Ey = 0). Ruled out MPI rank count, grid resolution, and the EB implicit-function form. See details below.

Additional information

Root cause and full evidence in the collapsible sections below.

Root cause (near-singular coarse operator + multi-box bottom solve)

Because warpx.eb_potential is set, WarpX does one electrostatic MLMG solve at t = 0 (operator amrex::MLEBNodeFDLaplacian, via ablastr::fields::computePhi).

• Neumann-x + periodic-y,z ⇒ the Laplacian has a constant null mode, pinned only by the EB Dirichlet electrodes.
• MLEBNodeFDLaplacian limits geometric coarsening so the EB stays represented ⇒ shallow MG hierarchy (2–3 levels) ⇒ the bottom solve runs on a still-sizable, EB-cut coarse grid.
• As one box, the default BiCGStab bottom solver handles it (0 failures, Ey = 0). Split into multiple boxes (smaller max_grid_size, or more MPI ranks via refine_grid_layout), the agglomerated coarse operator becomes ill-conditioned for
  BiCGStab and the bottom solve fails.
• With the coarse correction discarded, MG falls back to smoothing only. Smoothers reduce high-frequency error (residual still drops) but not the smooth transverse mode — which survives as the spurious Ey/Ez (residual ≠ error).

Evidence (box-count scan; CPU vs GPU)

max|Ex| ≈ 9.6e5 V/m in all runs; physical Ey/Ez should be 0.

Single rank, varying amr.max_grid_size (no mpirun), dx = 2e-4 (n_cell = 528 80 8):

max_grid_size	boxes	bottom-solve failures	max|Ey| (V/m)
528	1	0	0.0 (exact)
256	3	6	2.5e-7
128	5	7	2.8e-7
64	18	7	3.1e-7
32	51	7	4.1e-7

Same layout (4 ranks / 4 boxes), CPU vs GPU — same trigger, very different severity:

grid	backend	failures	MLMG iters	max|Ey|
2e-4	GPU	34	34	2.8e-3
2e-4	CPU	0	7	0.0
1e-4	GPU	42	42	7.2e-3
1e-4	CPU	7	7	4.5e-7

On CPU the failed bottom solve is benign (MLMG still converges in a few iterations); on GPU it collapses convergence and leaves a ~1e-3 V/m artifact. The GPU-specific cause is still under investigation.

What it is NOT (ruled out: rank count, resolution, implicit function)

• Not MPI rank count: one rank with a small max_grid_size (>1 box) reproduces it; rank count only amplifies severity.
• Not resolution: both dx = 2e-4 and 1e-4 fail once there is >1 box.
• Not the implicit-function form: replacing the step function
  2*(x<=-0.052)+2*(x>=0.052)-(x>-0.052)*(x<0.052) with the smooth signed distance abs(x)-0.052 (identical geometry) does not help — bottom-solve failures persist and the artifact is unchanged (the electrodes are node-aligned, so both give the same cut).

Possible fix / workarounds

The bottom solve on this near-singular, EB-cut, multi-box coarse operator needs a more robust solver than the default BiCGStab. ablastr/fields/PoissonSolver.H constructs amrex::MLMG with the default bottom solver and exposes no knob; adding e.g.
warpx.self_fields_bottom_solver = {bicgstab, cg, hypre, smoother} wired to MLMG::setBottomSolver (default unchanged; hypre needs -DAMReX_HYPRE=ON) would let users switch to hypre (BoomerAMG). Workaround without a code change: tighten warpx.self_fields_required_precision (reduces, does not remove, the artifact).

[n01r]

Here is an image showing the spurious fields in a comparison of the different decompositions / resolutions.

[n01r]

Trying out the algebraic solver from Hypre using AMReX-Codes/amrex#5649, to see if the bottom-solves will succeed.

[n01r]

Good news! It seems like doing the bottom solve with Hypre's algebraic solver avoids the failures and reduces the artifacts down to below the tolerance level at around ~1e-7 V/cm!
I just was not able to run it on GPU yet, due to Perlmutter being down at the moment, and my local machine having issues due to WSL+CUDA+Hypre where not all features are supported, yet.

[n01r]

Confirming: using Hypre as the bottom-solver works on GPU, avoids the failures and reduces the artifacts! 🥳

[RemiLehe]

Thanks a lot for the bug report!

I am not sure I understand this part though:

As one box, the default BiCGStab bottom solver handles it (0 failures, Ey = 0). Split into multiple boxes (smaller max_grid_size, or more MPI ranks via refine_grid_layout), the agglomerated coarse operator becomes ill-conditioned for BiCGStab.

The fact that it is one box vs multiple box should only change how the problem is distributed, but not its condition number no? i.e. the bottom solver is still mathematically solving the exact same linear system of equations (at the bottom level), albeit in a distributed way: is that correct?

Or is there maybe a bug in AMReX that makes the linear system (at the bottom level) mathematically different depending on whether it is distributed or on a single box?

[WeiqunZhang]

AMReX-Codes/amrex#5657 I think this fixes it.