From 282ed2c38b913dd3e9c4a0b46fa29b2bf22e0e59 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 17:32:21 -0700 Subject: [PATCH] [fdmath] correct Gridlock Yee spacing order Use center-to-center intervals for forward E derivatives and base-cell widths for backward H derivatives in the Gridlock examples. Document the distinction from meshes defined directly by forward intervals and fix the random fixture's backward averaging direction. Add coordinate-derivative regression coverage for graded periodic and nonperiodic Gridlock meshes using both functional and sparse operators. Validation: 1689 tests passed, 586 skipped; no new Ruff findings in changed Python files. --- docs/api/fdmath.md | 13 +++++++++++++ examples/eme.py | 2 +- examples/eme_bend.py | 2 +- examples/fdfd0.py | 2 +- examples/fdfd1.py | 2 +- examples/fdtd.py | 2 +- examples/tcyl.py | 2 +- examples/waveguide.py | 4 ++-- meanas/test/conftest.py | 2 +- meanas/test/test_fdmath_functional.py | 28 +++++++++++++++++++++++++++ 10 files changed, 50 insertions(+), 9 deletions(-) diff --git a/docs/api/fdmath.md b/docs/api/fdmath.md index 3623207..b3011a6 100644 --- a/docs/api/fdmath.md +++ b/docs/api/fdmath.md @@ -1,5 +1,18 @@ # fdmath +For Gridlock's default E-field Yee shifts, construct solver spacings as +`dxes = [grid.autoshifted_dxyz(), grid.dxyz]`. The forward E curl needs +distances between adjacent base-cell centers; the backward H curl needs +distances between adjacent upper cell edges. Reversing these two width families +introduces an error on graded grids that uniform grids hide. + +If constructing the mesh directly from **forward derivative intervals** `d`, +with H samples at their midpoints, the backward intervals instead follow +`(d + numpy.roll(d, 1)) / 2` for periodic boundaries. These `d` arrays are +already center-to-center distances, whereas `grid.dxyz` contains base-cell +volume widths. Applying this averaging rule to `grid.dxyz` does not fix the +ordering error above. + ::: meanas.fdmath ## Functional and sparse operators diff --git a/examples/eme.py b/examples/eme.py index 3a26dc8..547a75e 100644 --- a/examples/eme.py +++ b/examples/eme.py @@ -109,7 +109,7 @@ def build_geometry( z=Extent(min=0, max=1e6), ) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] dxes_2d = [[d[0], d[1]] for d in dxes] return grid, epsilon, dxes_2d, omega diff --git a/examples/eme_bend.py b/examples/eme_bend.py index 75be3e1..0a4cf74 100644 --- a/examples/eme_bend.py +++ b/examples/eme_bend.py @@ -99,7 +99,7 @@ def build_geometry( z=Extent(min=-1e6, center=0), ) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] dxes_2d = [[d[0], d[1]] for d in dxes] return grid, epsilon, dxes_2d, omega diff --git a/examples/fdfd0.py b/examples/fdfd0.py index a6227c5..fbffc02 100644 --- a/examples/fdfd0.py +++ b/examples/fdfd0.py @@ -67,7 +67,7 @@ def test0(solver=generic_solver): num_points = 24, ) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] for a in (0, 1, 2): for p in (-1, 1): dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, diff --git a/examples/fdfd1.py b/examples/fdfd1.py index f983d37..39649c3 100644 --- a/examples/fdfd1.py +++ b/examples/fdfd1.py @@ -50,7 +50,7 @@ def test1(solver=generic_solver): epsilon = grid.allocate(n_air**2, dtype=numpy.float32) grid.draw_cuboid(epsilon, x=dict(center=0, span=8e3), y=dict(center=0, span=w), z=dict(center=0, span=th), foreground=n_wg**2) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] for a in (0, 1, 2): for p in (-1, 1): dxes = scpml.stretch_with_scpml(dxes,omega=omega, axis=a, polarity=p, diff --git a/examples/fdtd.py b/examples/fdtd.py index de7b33d..cfc82cd 100644 --- a/examples/fdtd.py +++ b/examples/fdtd.py @@ -141,7 +141,7 @@ def main() -> None: ee = numpy.zeros_like(epsilon, dtype=complex) hh = numpy.zeros_like(epsilon, dtype=complex) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] # PMLs in every direction pml_params = [ diff --git a/examples/tcyl.py b/examples/tcyl.py index e6dc15f..7b20d11 100644 --- a/examples/tcyl.py +++ b/examples/tcyl.py @@ -46,7 +46,7 @@ def test1(solver=generic_solver): epsilon = grid.allocate(n_air**2, dtype=numpy.float32) grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], foreground=n_wg**2) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] for a in (1, 2): for p in (-1, 1): dxes = scpml.stretch_with_scpml( diff --git a/examples/waveguide.py b/examples/waveguide.py index 7301c6c..bdf7e4c 100644 --- a/examples/waveguide.py +++ b/examples/waveguide.py @@ -119,7 +119,7 @@ def main( # Add SCPML stretching to the FDFD grid; this matches the CPML-backed FDTD # run below so the two solvers see the same absorbing boundary profile. - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] for a in (0, 1, 2): for p in (-1, 1): dxes = scpml.stretch_with_scpml(dxes, omega=omega, axis=a, polarity=p, thickness=pml_thickness) @@ -269,7 +269,7 @@ def main2(): ee = numpy.zeros_like(epsilon, dtype=complex) hh = numpy.zeros_like(epsilon, dtype=complex) - dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes = [grid.autoshifted_dxyz(), grid.dxyz] # PMLs in every direction pml_params = [ diff --git a/meanas/test/conftest.py b/meanas/test/conftest.py index 5b577c6..11d1d77 100644 --- a/meanas/test/conftest.py +++ b/meanas/test/conftest.py @@ -93,6 +93,6 @@ def dxes( dxes[eh][ax][dxes[eh][ax].size // 2] *= 1.1 elif request.param == 'random': dxe = [prng.uniform(low=1.0 * dx, high=1.1 * dx, size=s) for s in shape[1:]] - dxh = [(d + numpy.roll(d, -1)) / 2 for d in dxe] + dxh = [(d + numpy.roll(d, 1)) / 2 for d in dxe] dxes = [dxe, dxh] return dxes diff --git a/meanas/test/test_fdmath_functional.py b/meanas/test/test_fdmath_functional.py index 368eba3..f6d1491 100644 --- a/meanas/test/test_fdmath_functional.py +++ b/meanas/test/test_fdmath_functional.py @@ -1,4 +1,6 @@ import numpy +import pytest +from gridlock import Grid from ..fdmath import functional as fd_functional from ..fdmath import operators as fd_operators @@ -58,3 +60,29 @@ def test_curls_match_sparse_operators_on_nonuniform_grid() -> None: assert_fields_close(curl_forward, matrix_forward, atol=1e-12, rtol=1e-12) assert_fields_close(curl_back, matrix_back, atol=1e-12, rtol=1e-12) + + +@pytest.mark.parametrize('periodic', [True, False]) +def test_gridlock_spacings_differentiate_physical_coordinates(periodic) -> None: + grid = Grid([[0, 1, 3, 7, 15], [-2, 0, 5, 6], [3, 6, 10]], periodic=periodic) + # Transverse E samples lie at cell centers and H samples at upper cell edges. + dxes = [grid.autoshifted_dxyz(), grid.dxyz] + for forward in (True, False): + if forward: + coordinates = grid.xyz + functions = fd_functional.deriv_forward(dxes[0]) + matrices = fd_operators.deriv_forward(dxes[0]) + else: + coordinates = [edges[1:] for edges in grid.exyz] + functions = fd_functional.deriv_back(dxes[1]) + matrices = fd_operators.deriv_back(dxes[1]) + for axis in range(3): + shape = [1, 1, 1] + shape[axis] = len(coordinates[axis]) + field = numpy.broadcast_to(coordinates[axis].reshape(shape), grid.shape) + interior = [slice(None)] * 3 + # A linear coordinate is discontinuous at the periodic wrap. + interior[axis] = slice(None, -1) if forward else slice(1, None) + assert_close(functions[axis](field)[tuple(interior)], 1) + matrix_result = (matrices[axis] @ field.ravel()).reshape(grid.shape) + assert_close(matrix_result[tuple(interior)], 1)