MPI-dependent deposition for 2D point and zero-height line sources on chunk boundaries

[meNOwanTstudy]

MPI-dependent deposition for 2D point/line sources on chunk boundaries

Summary

Ordinary mp.Source deposition for 2D zero-size point
sources and zero-height line sources appears to depend on the MPI chunk
decomposition when the discrete source support crosses an internal chunk
boundary.

The physical setup is unchanged: same cell, source, resolution, and
geometry_center. Only the MPI process count changes. I expected
Simulation.get_source(mp.Ez) to produce the same discrete support and total
source weight up to roundoff.

I have only checked this in 2D; I have not checked the corresponding 1D or 3D cases.

Minimal reproducer

import meep as mp
import numpy as np


resolution = 40
cell_size = mp.Vector3(9.85, 18.75)
geometry_center = mp.Vector3(0.0125, 0.0125)

cases = {
    "point_xy_boundary": (mp.Vector3(0.0125, 4.6875), mp.Vector3()),
    "hline_y_boundary": (mp.Vector3(0.0125, 4.6875), mp.Vector3(7.85, 0)),
    "hline_y_off_boundary_control": (mp.Vector3(0.0125, 4.0), mp.Vector3(7.85, 0)),
}

for name, (center, size) in cases.items():
    sim = mp.Simulation(
        cell_size=cell_size,
        boundary_layers=[mp.PML(1.0)],
        sources=[
            mp.Source(
                mp.GaussianSource(frequency=1.0, fwidth=0.5),
                component=mp.Ez,
                center=center,
                size=size,
            )
        ],
        resolution=resolution,
        geometry_center=geometry_center,
    )
    sim.init_sim()

    source = np.asarray(sim.get_source(mp.Ez))
    x_coords, y_coords, _, _ = sim.get_array_metadata()
    nz = np.argwhere(np.abs(source) > 1e-14)
    x_support = {round(float(x_coords[i]), 12) for i, _ in nz}
    y_support = {round(float(y_coords[j]), 12) for _, j in nz}

    if mp.am_master():
        print(
            f"nproc={mp.count_processors()} case={name} "
            f"nnz={len(nz)} x_support={len(x_support)} "
            f"y_support={len(y_support)} l1={np.sum(np.abs(source)):.12g}"
        )
    sim.reset_meep()

Run:

python repro_2d_zero_size_source_mpi.py
mpirun -np 4 python repro_2d_zero_size_source_mpi.py
mpirun -np 8 python repro_2d_zero_size_source_mpi.py

Observed output

case=point_xy_boundary
nproc=1  nnz=4    x_support=2    y_support=2    l1=1600
nproc=4  nnz=2    x_support=1    y_support=2    l1=800
nproc=8  nnz=1    x_support=1    y_support=1    l1=400

case=hline_y_boundary
nproc=1  nnz=630  x_support=315  y_support=2    l1=12560
nproc=4  nnz=630  x_support=315  y_support=2    l1=12560
nproc=8  nnz=315  x_support=315  y_support=1    l1=6280

case=hline_y_off_boundary_control
nproc=1  nnz=315  x_support=315  y_support=1    l1=12560
nproc=4  nnz=315  x_support=315  y_support=1    l1=12560
nproc=8  nnz=315  x_support=315  y_support=1    l1=12560

For the point source at (0.0125, 4.6875), the serial run deposits onto four
Yee-grid samples. With MPI, support points are dropped when internal chunk
boundaries cross that support: nproc=4 halves the total source weight, and
nproc=8 reduces it to one quarter.

For the zero-height horizontal line source at the same y, the serial support
has two y rows. At nproc=8, one row is omitted. Moving the same line source
away from the chunk boundary keeps the deposition invariant.

Additional observations

The same issue is not limited to exactly zero width in x. For a dx-wide
zero-height source at the same center, the deposited l1 also depends on the
MPI decomposition, but not monotonically: l1=40,40,60,30 for nproc=1,2,4,8.

This can affect adjoint workflows as well. In one fixed
mpa.OptimizationProblem([x]) evaluation with the same boundary point source,
the objective changes from 3.1461 (nproc=1) to 0.6659 (nproc=4) and
0.1682 (nproc=8). Moving the point source away from the boundary keeps the
objective essentially invariant.

Question

I noticed that the Python IndexedSource path calls fields.add_srcdata with a
needs_boundary_fix argument, while ordinary mp.Source uses
fields.add_volume_source. Is that boundary-fix logic intended only for
IndexedSource, or should ordinary mp.Source deposition have analogous
handling when its discrete support spans multiple MPI chunks?