249 lines
8.7 KiB
Python
249 lines
8.7 KiB
Python
from typing import Any
|
|
|
|
import numpy
|
|
import pytest
|
|
|
|
from .. import fdtd
|
|
from ..fdtd.base import maxwell_e, maxwell_h
|
|
from ..fdtd.pml import cpml_params, updates_with_cpml
|
|
from .utils import assert_close
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
('axis', 'polarity', 'thickness', 'epsilon_eff'),
|
|
[(3, 1, 4, 1.0), (0, 0, 4, 1.0), (0, 1, 2, 1.0), (0, 1, 4, 0.0)],
|
|
)
|
|
def test_cpml_params_reject_invalid_arguments(axis: int, polarity: int, thickness: int, epsilon_eff: float) -> None:
|
|
with pytest.raises(ValueError, match='Invalid axis|Invalid polarity|wise to have a pml|epsilon_eff must be positive'):
|
|
cpml_params(axis=axis, polarity=polarity, dt=0.1, thickness=thickness, epsilon_eff=epsilon_eff)
|
|
|
|
|
|
def test_cpml_params_shapes_and_region() -> None:
|
|
params = cpml_params(axis=1, polarity=1, dt=0.1, thickness=3)
|
|
p0e, p1e, p2e = params['param_e']
|
|
p0h, p1h, p2h = params['param_h']
|
|
|
|
assert p0e.shape == (1, 3, 1)
|
|
assert p1e.shape == (1, 3, 1)
|
|
assert p2e.shape == (1, 3, 1)
|
|
assert p0h.shape == (1, 3, 1)
|
|
assert p1h.shape == (1, 3, 1)
|
|
assert p2h.shape == (1, 3, 1)
|
|
assert params['region'][1] == slice(-3, None)
|
|
|
|
|
|
def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
|
|
shape = (3, 4, 4, 4)
|
|
epsilon = numpy.ones(shape, dtype=float)
|
|
e = numpy.zeros(shape, dtype=float)
|
|
h = numpy.zeros(shape, dtype=float)
|
|
dxes = [[numpy.ones(4), numpy.ones(4), numpy.ones(4)] for _ in range(2)]
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=3)
|
|
|
|
update_e, update_h = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
|
|
update_e(e, h, epsilon)
|
|
update_h(e, h)
|
|
|
|
assert not e.any()
|
|
assert not h.any()
|
|
|
|
|
|
def _unit_dxes(shape: tuple[int, int, int, int]) -> list[list[numpy.ndarray]]:
|
|
return [[numpy.ones(n, dtype=float) for n in shape[1:]] for _ in range(2)]
|
|
|
|
|
|
def _real_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray:
|
|
total = numpy.prod(shape, dtype=int)
|
|
return numpy.arange(start, start + total, dtype=float).reshape(shape) / total
|
|
|
|
|
|
def _complex_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray:
|
|
real = _real_field(shape, start)
|
|
imag = _real_field(shape, start + real.size)
|
|
return real + 1j * imag
|
|
|
|
|
|
def test_updates_with_cpml_matches_base_updates_when_all_faces_disabled() -> None:
|
|
shape = (3, 4, 5, 6)
|
|
epsilon = _real_field(shape, 1.0) + 2.0
|
|
mu = _real_field(shape, 4.0) + 1.5
|
|
e = _real_field(shape, 10.0)
|
|
h = _real_field(shape, 100.0)
|
|
dxes = _unit_dxes(shape)
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
|
|
update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
|
|
update_e_base = maxwell_e(dt=0.1, dxes=dxes)
|
|
update_h_base = maxwell_h(dt=0.1, dxes=dxes)
|
|
|
|
e_cpml = e.copy()
|
|
h_cpml = h.copy()
|
|
e_base = e.copy()
|
|
h_base = h.copy()
|
|
|
|
update_e_cpml(e_cpml, h_cpml, epsilon)
|
|
update_e_base(e_base, h_base, epsilon)
|
|
update_h_cpml(e_cpml, h_cpml, mu)
|
|
update_h_base(e_base, h_base, mu)
|
|
|
|
assert_close(e_cpml, e_base)
|
|
assert_close(h_cpml, h_base)
|
|
|
|
|
|
def test_updates_with_cpml_matches_base_updates_with_complex_dtype_when_all_faces_disabled() -> None:
|
|
shape = (3, 4, 5, 6)
|
|
epsilon = _real_field(shape, 1.0) + 2.0
|
|
mu = _real_field(shape, 4.0) + 1.5
|
|
e = _complex_field(shape, 10.0)
|
|
h = _complex_field(shape, 100.0)
|
|
dxes = _unit_dxes(shape)
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
|
|
update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon, dtype=complex)
|
|
update_e_base = maxwell_e(dt=0.1, dxes=dxes)
|
|
update_h_base = maxwell_h(dt=0.1, dxes=dxes)
|
|
|
|
e_cpml = e.copy()
|
|
h_cpml = h.copy()
|
|
e_base = e.copy()
|
|
h_base = h.copy()
|
|
|
|
update_e_cpml(e_cpml, h_cpml, epsilon)
|
|
update_e_base(e_base, h_base, epsilon)
|
|
update_h_cpml(e_cpml, h_cpml, mu)
|
|
update_h_base(e_base, h_base, mu)
|
|
|
|
assert_close(e_cpml, e_base)
|
|
assert_close(h_cpml, h_base)
|
|
|
|
|
|
def test_updates_with_cpml_only_changes_the_configured_face_region() -> None:
|
|
shape = (3, 6, 6, 6)
|
|
epsilon = numpy.ones(shape, dtype=float)
|
|
mu = numpy.ones(shape, dtype=float)
|
|
e = _real_field(shape, 1.0)
|
|
h = _real_field(shape, 100.0)
|
|
dxes = _unit_dxes(shape)
|
|
thickness = 3
|
|
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=thickness)
|
|
|
|
update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
|
|
update_e_base = maxwell_e(dt=0.1, dxes=dxes)
|
|
update_h_base = maxwell_h(dt=0.1, dxes=dxes)
|
|
|
|
e_cpml = e.copy()
|
|
h_cpml = h.copy()
|
|
e_base = e.copy()
|
|
h_base = h.copy()
|
|
|
|
update_e_cpml(e_cpml, h_cpml, epsilon)
|
|
update_e_base(e_base, h_base, epsilon)
|
|
update_h_cpml(e_cpml, h_cpml, mu)
|
|
update_h_base(e_base, h_base, mu)
|
|
|
|
e_untouched = slice(thickness, None)
|
|
h_untouched = slice(thickness, -1)
|
|
assert_close(e_cpml[:, e_untouched, :, :], e_base[:, e_untouched, :, :])
|
|
assert_close(h_cpml[:, h_untouched, :, :], h_base[:, h_untouched, :, :])
|
|
|
|
changed_e = numpy.any(numpy.abs(e_cpml[:, :thickness, :, :] - e_base[:, :thickness, :, :]) > 1e-12)
|
|
changed_h = numpy.any(numpy.abs(h_cpml[:, :thickness, :, :] - h_base[:, :thickness, :, :]) > 1e-12)
|
|
assert changed_e
|
|
assert changed_h
|
|
|
|
|
|
def test_cpml_plane_wave_phasor_decays_monotonically_through_outgoing_pml() -> None:
|
|
dt = 0.4
|
|
period_steps = 24
|
|
omega = 2 * numpy.pi / (period_steps * dt)
|
|
shape = (3, 80, 1, 1)
|
|
thickness = 8
|
|
source_x = 16
|
|
warmup_periods = 10
|
|
accumulation_periods = 6
|
|
total_steps = period_steps * (warmup_periods + accumulation_periods)
|
|
|
|
epsilon = numpy.ones(shape, dtype=float)
|
|
dxes = _unit_dxes(shape)
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
for polarity_index, polarity in enumerate((-1, 1)):
|
|
params[0][polarity_index] = cpml_params(axis=0, polarity=polarity, dt=dt, thickness=thickness)
|
|
|
|
update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon)
|
|
|
|
e = numpy.zeros(shape, dtype=float)
|
|
h = numpy.zeros(shape, dtype=float)
|
|
e_accumulator = numpy.zeros((1, *shape), dtype=complex)
|
|
|
|
for step in range(total_steps):
|
|
update_e(e, h, epsilon)
|
|
|
|
source = numpy.cos(omega * (step + 0.5) * dt)
|
|
e[1, source_x, 0, 0] -= dt * source
|
|
|
|
if step >= period_steps * warmup_periods:
|
|
fdtd.accumulate_phasor_e(e_accumulator, omega, dt, e, step + 1)
|
|
|
|
update_h(e, h)
|
|
|
|
profile = numpy.abs(e_accumulator[0, 1, :, 0, 0])
|
|
right_pml = profile[-thickness:]
|
|
interior = profile[-thickness - 6:-thickness]
|
|
interior_level = interior.mean()
|
|
|
|
assert interior_level > 1.0
|
|
assert right_pml[-1] < interior_level / 100
|
|
assert profile[0] < interior_level / 100
|
|
assert numpy.all(numpy.diff(right_pml) <= interior_level * 1e-3)
|
|
|
|
|
|
@pytest.mark.complete
|
|
def test_cpml_point_source_total_energy_reaches_late_time_plateau() -> None:
|
|
dt = 0.3
|
|
period_steps = 24
|
|
omega = 2 * numpy.pi / (period_steps * dt)
|
|
cycles = 1000
|
|
sample_every_cycles = 50
|
|
sample_stride = period_steps * sample_every_cycles
|
|
shape = (3, 9, 9, 9)
|
|
thickness = 3
|
|
center = shape[1] // 2
|
|
|
|
epsilon = numpy.ones(shape, dtype=float)
|
|
dxes = _unit_dxes(shape)
|
|
params: list[list[dict[str, Any] | None]] = [[None, None] for _ in range(3)]
|
|
for axis in range(3):
|
|
for polarity_index, polarity in enumerate((-1, 1)):
|
|
params[axis][polarity_index] = cpml_params(axis=axis, polarity=polarity, dt=dt, thickness=thickness)
|
|
|
|
update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon)
|
|
|
|
e = numpy.zeros(shape, dtype=float)
|
|
h = numpy.zeros(shape, dtype=float)
|
|
sampled_energies: list[float] = []
|
|
|
|
for step in range(period_steps * cycles):
|
|
h_before = h.copy()
|
|
update_e(e, h, epsilon)
|
|
|
|
source = numpy.cos(omega * (step + 0.5) * dt)
|
|
e[1, center, center, center] -= dt * source
|
|
|
|
update_h(e, h)
|
|
|
|
if (step + 1) % sample_stride == 0:
|
|
total_energy = fdtd.energy_estep(h0=h_before, e1=e, h2=h, epsilon=epsilon, dxes=dxes).sum().real
|
|
sampled_energies.append(total_energy)
|
|
|
|
energies = numpy.asarray(sampled_energies)
|
|
late_window = energies[-5:]
|
|
previous_window = energies[-10:-5]
|
|
late_mean = late_window.mean()
|
|
|
|
assert energies.size == cycles // sample_every_cycles
|
|
assert late_mean > 0.1
|
|
assert (late_window.max() - late_window.min()) / late_mean < 1e-4
|
|
assert abs(late_mean - previous_window.mean()) / late_mean < 1e-4
|