diff --git a/docs/api/fdfd.md b/docs/api/fdfd.md index 83328f0..20b126f 100644 --- a/docs/api/fdfd.md +++ b/docs/api/fdfd.md @@ -27,3 +27,10 @@ stretching can produce above-light-line PML artifacts. The cap is not a general spectral guarantee, and weaker absorption can increase boundary reflections. Check domain size, PML thickness, strength, mode localization, and residuals. The general scattering profile remains uncapped by default. + +`uniform_grid_scpml` builds unit widths by applying `stretch_with_scpml` to +both faces of each enabled axis. Both builders divide imaginary grading by +physical layer thickness and sample electric widths at cell centers and +magnetic widths at cell edges. This corrects the uniform builder's previous +thickness-independent strength and reversed E/H sample offsets; retune existing +uniform-builder profiles using domain/PML convergence checks. diff --git a/meanas/fdfd/scpml.py b/meanas/fdfd/scpml.py index acbcbdd..eaa7c5a 100644 --- a/meanas/fdfd/scpml.py +++ b/meanas/fdfd/scpml.py @@ -82,6 +82,8 @@ def uniform_grid_scpml( thicknesses; zero disables an axis. ``s_function`` specifies the imaginary grading, defaulting to prepare_s_function(). ``kappa_function`` optionally returns dimensionless real stretch values >= 1 at normalized depths. + Each face is built by stretch_with_scpml, including its physical-thickness + scaling and electric-center / magnetic-edge sampling convention. ``max_imag_ratio`` optionally bounds abs(Im(s))/Re(s) after frequency scaling. A value of 1 is a useful mode-solver starting point, not a proof @@ -101,15 +103,15 @@ def uniform_grid_scpml( if s_function is None: s_function = prepare_s_function() result = [[numpy.ones(n) for n in shape] for _ in range(2)] - correction = numpy.sqrt(epsilon_effective) * numpy.real(omega) - for axis, (n, thickness) in enumerate(zip(shape, thicknesses, strict=True)): + for axis, thickness in enumerate(thicknesses): if not thickness: continue - for family, offset in enumerate((0.0, 0.5)): - position = numpy.arange(n, dtype=float) + offset - depth = ((thickness - position).clip(0) + (position - (n - thickness)).clip(0)) / thickness - factors = _stretch_factors(depth, correction, s_function, kappa_function, max_imag_ratio) - result[family][axis] = numpy.where(depth > 0, factors, 1.0) + for polarity in (-1, 1): + stretch_with_scpml( + result, axis, polarity, omega, epsilon_effective, thickness, + s_function, kappa_function=kappa_function, + max_imag_ratio=max_imag_ratio, copy=False, + ) return result diff --git a/meanas/test/test_fdfd_algebra_helpers.py b/meanas/test/test_fdfd_algebra_helpers.py index bf863a9..bf6e973 100644 --- a/meanas/test/test_fdfd_algebra_helpers.py +++ b/meanas/test/test_fdfd_algebra_helpers.py @@ -51,10 +51,6 @@ def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray: return pm * (curl_term - OMEGA**2 * magnetic * masked_h) -def _normalized_distance(u: numpy.ndarray, size: int, thickness: int) -> numpy.ndarray: - return ((thickness - u).clip(0) + (u - (size - thickness)).clip(0)) / thickness - - def test_h_full_matches_dense_reference_with_and_without_mu() -> None: for mu in (None, MU): matrix_result = apply_fdfd_matrix( @@ -155,14 +151,16 @@ def test_prepare_s_function_matches_closed_form_polynomial() -> None: def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None: s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0) dxes = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0, epsilon_effective=4.0, s_function=s_function) - correction = numpy.sqrt(4.0) * 2.0 - for axis, size, thickness in ((0, 6, 2), (2, 3, 1)): - grid = numpy.arange(size, dtype=float) - expected_a = 1 + 1j * s_function(_normalized_distance(grid, size, thickness)) / correction - expected_b = 1 + 1j * s_function(_normalized_distance(grid + 0.5, size, thickness)) / correction - assert_close(dxes[0][axis], expected_a) - assert_close(dxes[1][axis], expected_b) + correction = thickness * numpy.sqrt(4.0) * 2.0 + for family, offset in enumerate((0.5, 0.0)): + positions = numpy.arange(size, dtype=float) + offset + expected = numpy.ones(size, dtype=complex) + low_depth = (thickness - positions[:thickness]) / thickness + high_depth = (positions[-thickness:] - (size - thickness)) / thickness + expected[:thickness] += 1j * s_function(low_depth) / correction + expected[-thickness:] += 1j * s_function(high_depth) / correction + assert_close(dxes[family][axis], expected) assert_close(dxes[0][1], 1.0) assert_close(dxes[1][1], 1.0) diff --git a/meanas/test/test_numerical_feedback.py b/meanas/test/test_numerical_feedback.py index b773115..7c15091 100644 --- a/meanas/test/test_numerical_feedback.py +++ b/meanas/test/test_numerical_feedback.py @@ -373,3 +373,30 @@ def test_derivatives_do_not_subtract_across_independent_rows(forward): derivative = (deriv_forward if forward else deriv_back)()[2] with np.errstate(over='raise', invalid='raise'): assert_allclose(derivative(fields), 0) + + +@pytest.mark.parametrize(('shape', 'thicknesses'), [ + ((6, 8, 4), (0, 0, 0)), + ((6, 8, 4), (1, 1, 1)), + ((6, 8, 4), (2, 3, 0)), + ((6, 8, 4), (3, 4, 2)), +]) +@pytest.mark.parametrize('kappa', [None, 'polynomial', 'constant']) +@pytest.mark.parametrize('cap', [None, .7, 1.0]) +def test_scpml_builders_agree_including_real_profiles_and_touching_faces(shape, thicknesses, kappa, cap): + profile = scpml.prepare_s_function(ln_R=-9, m=2.5) + real_profile = { + None: None, + 'polynomial': lambda depth: 1 + 2 * depth**2, + 'constant': lambda depth: np.full_like(depth, 2), + }[kappa] + options = dict(omega=.8, epsilon_effective=2.25, s_function=profile, + kappa_function=real_profile, max_imag_ratio=cap) + uniform = scpml.uniform_grid_scpml(shape, thicknesses, **options) + explicit = grid(shape) + for axis, thickness in enumerate(thicknesses): + for polarity in (-1, 1): + explicit = scpml.stretch_with_scpml(explicit, axis, polarity, thickness=thickness, **options) + for family in range(2): + for axis in range(3): + assert_allclose(uniform[family][axis], explicit[family][axis], rtol=1e-14, atol=1e-14)