Compare commits
25 changed files with 342 additions and 938 deletions
|
|
@ -14,47 +14,6 @@
|
|||
|
||||
::: meanas.fdfd.farfield
|
||||
|
||||
## Bloch magnetic materials
|
||||
|
||||
Bloch state vectors store the transverse Fourier coefficients of magnetic flux
|
||||
`B = mu * H`. They equal `H` when permeability is omitted. Pass the same `mu`
|
||||
to `bloch.hmn_2_hxyz(...)` and `bloch.hmn_2_exyz(...)` when converting solved
|
||||
vectors to fields. The eigensolver uses the projected `mu**-1` energy metric
|
||||
for finite positive real permeability, including anisotropic and spatially
|
||||
varying profiles. Its preconditioner remains finite at the Gamma point.
|
||||
|
||||
## Far-field normalization
|
||||
|
||||
The near-to-far transform includes the cell area `dx * dy`, so normalized
|
||||
far-field amplitudes describe a Fourier integral and converge under spatial
|
||||
refinement. This changes amplitudes previously returned as unweighted FFT sums.
|
||||
The inverse applies the matching Fourier quadrature and returns spacings for
|
||||
the requested output shape. Both transforms center their origins at index
|
||||
`size // 2`; padding and cropping preserve that sample. The inverse does not
|
||||
modify input spectra.
|
||||
|
||||
Transforms use the library's `exp(+i omega t)` time convention: outgoing plane
|
||||
waves vary as `exp(-i k dot r)`, so the angular projection uses a positive
|
||||
Fourier phase. The inverse uses a negative phase. Far-field vector pairs are
|
||||
ordered as polar and azimuthal components; near-field pairs are Cartesian.
|
||||
|
||||
## Solver and port calibration
|
||||
|
||||
The default QMR wrapper raises `RuntimeError` on breakdown or failure to
|
||||
converge. PEC rows impose zero electric fields with unit diagonal entries;
|
||||
`solvers.generic` zeroes the corresponding source and initial-guess entries.
|
||||
H-only and EH operators impose analogous constraints for PMC locations.
|
||||
|
||||
FDFD port calibration uses the complete incoming-amplitude matrix from all
|
||||
excitations, solving `S @ incoming_amplitudes = raw_outgoing_amplitudes`.
|
||||
Diagnostics include these matrices and the incoming/outgoing flux Gram
|
||||
matrices. Guided powers retain cross terms between nonorthogonal modes.
|
||||
|
||||
Lossless modes at or above the propagation-grid band edge require a finer
|
||||
mesh for power normalization. Lorentz-normalized evanescent modes use the
|
||||
decaying spatial branch. Expansion on graded grids uses each electric
|
||||
component's own Yee intervals.
|
||||
|
||||
## SCPML conventions
|
||||
|
||||
`stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the
|
||||
|
|
|
|||
|
|
@ -1,44 +1,97 @@
|
|||
"""
|
||||
Bloch eigenmode solvers and operators in a spatial Fourier basis.
|
||||
Bloch eigenmode solver/operators
|
||||
|
||||
The two stored coefficients at each wavevector represent transverse magnetic
|
||||
flux B = mu H. They coincide with H when mu is omitted. Using B ensures
|
||||
k dot B = 0 even in anisotropic or spatially varying magnetic materials; H
|
||||
itself need not be transverse. The historical h_mn function names are retained.
|
||||
Pass the same mu to hmn_2_hxyz and hmn_2_exyz when reconstructing these modes.
|
||||
This module contains functions for generating and solving the
|
||||
3D Bloch eigenproblem. The approach is to transform the problem
|
||||
into the (spatial) fourier domain, transforming the equation
|
||||
|
||||
For each Fourier wavevector, generate_kmn constructs unit vectors m and n with
|
||||
k cross m = |k| n. The Cartesian Fourier field is B_k = b_m m + b_n n.
|
||||
The Maxwell operator applies, in order:
|
||||
1/mu * curl(1/eps * curl(H_eigenmode)) = (w/c)^2 H_eigenmode
|
||||
|
||||
H_k = fftn(ifftn(B_k) / mu)
|
||||
E_k = fftn(ifftn(k cross H_k) / epsilon)
|
||||
A(B_k) = k cross E_k
|
||||
into
|
||||
|
||||
Material division is performed in real space and curls in Fourier space.
|
||||
The explicit factors of i are omitted, giving the historical negative squared
|
||||
frequency eigenvalues. Wavevectors and frequencies use cycles per length:
|
||||
frequency = sqrt(-real(eigenvalue)). Reconstructed electric fields retain the
|
||||
historical overall frequency scaling and are not power-normalized.
|
||||
conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
|
||||
|
||||
For finite positive real permeability, the flux operator is self-adjoint in
|
||||
the projected magnetic energy metric M = P mu^-1 P, where P projects onto the
|
||||
transverse Fourier subspace. The block variational solve uses that metric in
|
||||
its Rayleigh quotients, gradients, line searches, and final Ritz extraction.
|
||||
The inverse-curl preconditioner regularizes zero wavevectors so static modes
|
||||
remain accessible at Gamma. Nonzero-frequency modes have no static Fourier
|
||||
flux component at Gamma; the two stored zero-wavevector basis vectors choose
|
||||
two representatives of the static sector.
|
||||
where:
|
||||
|
||||
Typical usage:
|
||||
- the `_k` subscript denotes a 3D fourier transformed field
|
||||
- each component of `H_k` corresponds to a plane wave with wavevector `k`
|
||||
- `x` is the cross product
|
||||
- `conv()` denotes convolution
|
||||
|
||||
values, vectors = eigsolve(5, k0, reciprocal_lattice, epsilon, mu=mu)
|
||||
frequency = numpy.sqrt(-numpy.real(values[0]))
|
||||
electric = hmn_2_exyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0])
|
||||
magnetic = hmn_2_hxyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0])
|
||||
Since `k` and `H` are orthogonal for each plane wave, we can use each
|
||||
`k` to create an orthogonal basis (k, m, n), with `k x m = n`, and
|
||||
`|m| = |n| = 1`. The cross products are then simplified as follows:
|
||||
|
||||
- `h` is shorthand for `H_k`
|
||||
- `(...)_xyz` denotes the `(x, y, z)` basis
|
||||
- `(...)_kmn` denotes the `(k, m, n)` basis
|
||||
- `hm` is the component of `h` in the `m` direction, etc.
|
||||
|
||||
We know
|
||||
|
||||
k @ h = kx hx + ky hy + kz hz = 0 = hk
|
||||
h = hk + hm + hn = hm + hn
|
||||
k = kk + km + kn = kk = |k|
|
||||
|
||||
We can write
|
||||
|
||||
k x h = (ky hz - kz hy,
|
||||
kz hx - kx hz,
|
||||
kx hy - ky hx)_xyz
|
||||
= ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
|
||||
= (0, (m x k) @ h, (n x k) @ h)_kmn # triple product ordering
|
||||
= (0, kk (-n @ h), kk (m @ h))_kmn # (m x k) = -|k| n, etc.
|
||||
= |k| (0, -h @ n, h @ m)_kmn
|
||||
|
||||
which gives us a straightforward way to perform the cross product
|
||||
while simultaneously transforming into the `_kmn` basis.
|
||||
We can also write
|
||||
|
||||
k x h = (km hn - kn hm,
|
||||
kn hk - kk hn,
|
||||
kk hm - km hk)_kmn
|
||||
= (0, -kk hn, kk hm)_kmn
|
||||
= (-kk hn)(mx, my, mz)_xyz + (kk hm)(nx, ny, nz)_xyz
|
||||
= |k| (hm * (nx, ny, nz)_xyz
|
||||
- hn * (mx, my, mz)_xyz)
|
||||
|
||||
which gives us a way to perform the cross product while simultaneously
|
||||
trasnforming back into the `_xyz` basis.
|
||||
|
||||
We can also simplify `conv(X_k, Y_k)` as `fftn(X * ifftn(Y_k))`.
|
||||
|
||||
Using these results and storing `H_k` as `h = (hm, hn)`, we have
|
||||
|
||||
e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
|
||||
b_mn = |k| (-e_xyz @ n, e_xyz @ m)
|
||||
h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
|
||||
|
||||
which forms the operator from the left side of the equation.
|
||||
|
||||
We can then use a preconditioned block Rayleigh iteration algorithm, as in
|
||||
SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods
|
||||
for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001)
|
||||
(similar to that used in MPB) to find the eigenvectors for this operator.
|
||||
|
||||
===
|
||||
|
||||
Typically you will want to do something like
|
||||
|
||||
recip_lattice = numpy.diag(1/numpy.array(epsilon[0].shape * dx))
|
||||
n, v = bloch.eigsolve(5, k0, recip_lattice, epsilon)
|
||||
f = numpy.sqrt(-numpy.real(n[0]))
|
||||
n_eff = norm(recip_lattice @ k0) / f
|
||||
|
||||
v2e = bloch.hmn_2_exyz(k0, recip_lattice, epsilon)
|
||||
e_field = v2e(v[0])
|
||||
|
||||
k, f = find_k(frequency=1/1550,
|
||||
tolerance=(1/1550 - 1/1551),
|
||||
direction=[1, 0, 0],
|
||||
G_matrix=recip_lattice,
|
||||
epsilon=epsilon,
|
||||
band=0)
|
||||
|
||||
The block iteration follows the approach of S. G. Johnson and J. D.
|
||||
Joannopoulos, Optics Express 8, 173-190 (2001).
|
||||
"""
|
||||
|
||||
from typing import Any, cast
|
||||
|
|
@ -121,7 +174,7 @@ def generate_kmn(
|
|||
m = numpy.broadcast_to([0, 1, 0], tuple(shape[:3]) + (3,)).astype(float)
|
||||
n = numpy.broadcast_to([0, 0, 1], tuple(shape[:3]) + (3,)).astype(float)
|
||||
|
||||
xy_non0 = numpy.any(k_xyz[:, :, :, :2] != 0, axis=3)
|
||||
xy_non0 = numpy.any(k_xyz[:, :, :, 0:1] != 0, axis=3)
|
||||
if numpy.any(xy_non0):
|
||||
u = numpy.cross(k_xyz[xy_non0], [0, 0, 1])
|
||||
m[xy_non0, :] = u / norm(u, axis=1)[:, None]
|
||||
|
|
@ -142,20 +195,13 @@ def maxwell_operator(
|
|||
mu: fdfield | None = None
|
||||
) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]:
|
||||
"""
|
||||
Generate the Maxwell operator on transverse magnetic-flux coefficients.
|
||||
Generate the Maxwell operator
|
||||
|
||||
k x conv(1/eps_k, k x conv(1/mu_k, ___))
|
||||
conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
|
||||
|
||||
The corresponding real-space curl-curl problem is
|
||||
which is the spatial-frequency-space representation of
|
||||
|
||||
curl(1/eps * curl(1/mu * ___))
|
||||
|
||||
with the overall factor ``-(2*pi)**2`` omitted for the module's
|
||||
cycles-per-length frequency convention.
|
||||
|
||||
The state represents B = mu H, whose Fourier components are transverse.
|
||||
It equals H when mu is omitted. Eigenvalues retain the historical negative
|
||||
squared-frequency convention (the two explicit factors of i are omitted).
|
||||
1/mu * curl(1/eps * curl(___))
|
||||
|
||||
The operator is a function that acts on a vector h_mn of size `2 * epsilon[0].size`
|
||||
|
||||
|
|
@ -184,28 +230,23 @@ def maxwell_operator(
|
|||
"""
|
||||
Maxwell operator for Bloch eigenmode simulation.
|
||||
|
||||
h contains the two transverse magnetic-flux coefficients.
|
||||
h is complex 2-field in (m, n) basis, vectorized
|
||||
|
||||
Args:
|
||||
h: Raveled b_mn; size `2 * epsilon[0].size`. Not modified.
|
||||
h: Raveled h_mn; size `2 * epsilon[0].size`.
|
||||
Altered in-place.
|
||||
|
||||
Returns:
|
||||
Raveled transverse flux image under the Maxwell operator.
|
||||
Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)), returned
|
||||
and overwritten in-place of `h`.
|
||||
"""
|
||||
hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
|
||||
|
||||
#{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
|
||||
|
||||
# Convert transverse B into H before taking either curl. H itself need
|
||||
# not be transverse for anisotropic or spatially varying permeability.
|
||||
if mu is None:
|
||||
h_m, h_n = hin_m, hin_n
|
||||
else:
|
||||
b_xyz = m * hin_m + n * hin_n
|
||||
h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
|
||||
h_m = numpy.sum(h_xyz * m, axis=3, keepdims=True)
|
||||
h_n = numpy.sum(h_xyz * n, axis=3, keepdims=True)
|
||||
d_xyz = (n * h_m - m * h_n) * k_mag
|
||||
# cross product and transform into xyz basis
|
||||
d_xyz = (n * hin_m
|
||||
- m * hin_n) * k_mag
|
||||
|
||||
# divide by epsilon
|
||||
temp = ifftn(d_xyz, axes=range(3)) # reuses d_xyz if using pyfftw
|
||||
|
|
@ -216,7 +257,23 @@ def maxwell_operator(
|
|||
b_m = numpy.sum(e_xyz * n, axis=3, keepdims=True) * -k_mag
|
||||
b_n = numpy.sum(e_xyz * m, axis=3, keepdims=True) * +k_mag
|
||||
|
||||
return _assemble_hmn_vector(b_m, b_n)
|
||||
if mu is None:
|
||||
h_m, h_n = b_m, b_n
|
||||
else:
|
||||
# transform from mn to xyz
|
||||
b_xyz = (m * b_m
|
||||
+ n * b_n) # noqa
|
||||
|
||||
# divide by mu
|
||||
temp = ifftn(b_xyz, axes=range(3))
|
||||
temp /= mu
|
||||
h_xyz = fftn(temp, axes=range(3))
|
||||
|
||||
# transform back to mn
|
||||
h_m = numpy.sum(h_xyz * m, axis=3)
|
||||
h_n = numpy.sum(h_xyz * n, axis=3)
|
||||
|
||||
return _assemble_hmn_vector(h_m, h_n)
|
||||
|
||||
return operator
|
||||
|
||||
|
|
@ -225,12 +282,12 @@ def hmn_2_exyz(
|
|||
k0: ArrayLike,
|
||||
G_matrix: ArrayLike,
|
||||
epsilon: fdfield,
|
||||
mu: fdfield | None = None,
|
||||
) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]:
|
||||
"""
|
||||
Convert transverse Fourier flux coefficients into an unnormalized E field.
|
||||
Generate an operator which converts a vectorized spatial-frequency-space
|
||||
`h_mn` into an E-field distribution, i.e.
|
||||
|
||||
ifft(conv(1/eps_k, k x conv(1/mu_k, b_mn)))
|
||||
ifft(conv(1/eps_k, ik x h_mn))
|
||||
|
||||
The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
|
||||
|
||||
|
|
@ -241,26 +298,17 @@ def hmn_2_exyz(
|
|||
G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
||||
epsilon: Dielectric constant distribution for the simulation.
|
||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
||||
mu: Permeability used by the solve. The input coefficients describe B;
|
||||
divide by mu to recover H before computing the electric field.
|
||||
|
||||
Returns:
|
||||
Function for converting `h_mn` into `E_xyz`
|
||||
"""
|
||||
shape = epsilon[0].shape + (1,)
|
||||
epsilon = numpy.moveaxis(epsilon, 0, -1)
|
||||
if mu is not None:
|
||||
mu = numpy.moveaxis(mu, 0, -1)
|
||||
|
||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
||||
|
||||
def operator(h: NDArray[numpy.complex128]) -> cfdfield_t:
|
||||
hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
|
||||
if mu is not None:
|
||||
b_xyz = m * hin_m + n * hin_n
|
||||
h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
|
||||
hin_m = numpy.sum(h_xyz * m, axis=3, keepdims=True)
|
||||
hin_n = numpy.sum(h_xyz * n, axis=3, keepdims=True)
|
||||
d_xyz = (n * hin_m
|
||||
- m * hin_n) * k_mag
|
||||
|
||||
|
|
@ -275,12 +323,12 @@ def hmn_2_hxyz(
|
|||
k0: ArrayLike,
|
||||
G_matrix: ArrayLike,
|
||||
epsilon: fdfield,
|
||||
mu: fdfield | None = None,
|
||||
) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]:
|
||||
"""
|
||||
Convert transverse Fourier flux coefficients into a physical H field.
|
||||
Generate an operator which converts a vectorized spatial-frequency-space
|
||||
`h_mn` into an H-field distribution, i.e.
|
||||
|
||||
ifft(b_mn) / mu
|
||||
ifft(h_mn)
|
||||
|
||||
The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`.
|
||||
|
||||
|
|
@ -291,8 +339,6 @@ def hmn_2_hxyz(
|
|||
G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
||||
epsilon: Dielectric constant distribution for the simulation.
|
||||
Only `epsilon[0].shape` is used.
|
||||
mu: Permeability used by the solve. Input coefficients describe B;
|
||||
the returned physical H is B / mu. Default is unit permeability.
|
||||
|
||||
Returns:
|
||||
Function for converting `h_mn` into `H_xyz`
|
||||
|
|
@ -304,10 +350,7 @@ def hmn_2_hxyz(
|
|||
hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
|
||||
h_xyz = (m * hin_m
|
||||
+ n * hin_n)
|
||||
field = numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)])
|
||||
if mu is not None:
|
||||
field /= mu
|
||||
return cfdfield_t(field)
|
||||
return cfdfield_t(numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)]))
|
||||
|
||||
return operator
|
||||
|
||||
|
|
@ -319,16 +362,11 @@ def inverse_maxwell_operator_approx(
|
|||
mu: fdfield | None = None,
|
||||
) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]:
|
||||
"""
|
||||
Generate a positive approximate inverse of minus the flux Maxwell operator,
|
||||
Generate an approximate inverse of the Maxwell operator,
|
||||
|
||||
project_mn(conv(mu_k, cross_inv(conv(eps_k, cross_inv(___)))))
|
||||
ik x conv(eps_k, ik x conv(mu_k, ___))
|
||||
|
||||
At k = 0, inverse crosses use the smallest nonzero grid wavevector scale
|
||||
(or one on a single-point grid). This finite regularization keeps the
|
||||
preconditioner positive on static directions, which must remain accessible
|
||||
to the eigensolver at Gamma.
|
||||
|
||||
It is used for the preconditioned block variational solve.
|
||||
which can be used to improve the speed of ARPACK in shift-invert mode.
|
||||
|
||||
See the `meanas.fdfd.bloch` docstring for more information.
|
||||
|
||||
|
|
@ -349,9 +387,6 @@ def inverse_maxwell_operator_approx(
|
|||
mu = numpy.moveaxis(mu, 0, -1)
|
||||
|
||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
||||
nonzero = k_mag[k_mag != 0]
|
||||
reference = float(nonzero.min()) if nonzero.size else 1.0
|
||||
k_inverse = numpy.divide(1, k_mag, out=numpy.full_like(k_mag, 1 / reference), where=k_mag != 0)
|
||||
|
||||
def operator(h: NDArray[numpy.complex128]) -> NDArray[numpy.complex128]:
|
||||
"""
|
||||
|
|
@ -363,14 +398,31 @@ def inverse_maxwell_operator_approx(
|
|||
h: Raveled h_mn; size `2 * epsilon[0].size`.
|
||||
|
||||
Returns:
|
||||
Raveled, regularized approximate inverse flux image.
|
||||
Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
|
||||
"""
|
||||
hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
|
||||
|
||||
#{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
|
||||
|
||||
if mu is None:
|
||||
b_m, b_n = hin_m, hin_n
|
||||
else:
|
||||
# transform from mn to xyz
|
||||
h_xyz = (m * hin_m
|
||||
+ n * hin_n) # noqa
|
||||
|
||||
# multiply by mu
|
||||
temp = ifftn(h_xyz, axes=range(3))
|
||||
temp *= mu
|
||||
b_xyz = fftn(temp, axes=range(3))
|
||||
|
||||
# transform back to mn
|
||||
b_m = numpy.sum(b_xyz * m, axis=3, keepdims=True)
|
||||
b_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
|
||||
|
||||
# cross product and transform into xyz basis
|
||||
e_xyz = (n * hin_m - m * hin_n) * k_inverse
|
||||
e_xyz = (n * b_m
|
||||
- m * b_n) / k_mag
|
||||
|
||||
# multiply by epsilon
|
||||
temp = ifftn(e_xyz, axes=range(3))
|
||||
|
|
@ -378,41 +430,14 @@ def inverse_maxwell_operator_approx(
|
|||
d_xyz = fftn(temp, axes=range(3))
|
||||
|
||||
# cross product and transform into mn basis crossinv_t2c
|
||||
h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) * k_inverse
|
||||
h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) * -k_inverse
|
||||
if mu is not None:
|
||||
h_xyz = m * h_m + n * h_n
|
||||
b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3))
|
||||
h_m = numpy.sum(b_xyz * m, axis=3, keepdims=True)
|
||||
h_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
|
||||
h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) / +k_mag
|
||||
h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) / -k_mag
|
||||
|
||||
return _assemble_hmn_vector(h_m, h_n)
|
||||
|
||||
return operator
|
||||
|
||||
|
||||
def _transverse_material_operator(
|
||||
k0: ArrayLike,
|
||||
G_matrix: ArrayLike,
|
||||
material: fdfield,
|
||||
*,
|
||||
inverse: bool,
|
||||
) -> Callable[[NDArray], NDArray]:
|
||||
"""Apply P material^(+/-1) P on the transverse Fourier subspace."""
|
||||
shape = material[0].shape + (1,)
|
||||
_k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
||||
values = numpy.moveaxis(material, 0, -1)
|
||||
|
||||
def apply(h: NDArray) -> NDArray:
|
||||
hm, hn = (part.reshape(shape) for part in numpy.split(h, 2))
|
||||
field = ifftn(m * hm + n * hn, axes=range(3))
|
||||
field = field / values if inverse else field * values
|
||||
transformed = fftn(field, axes=range(3))
|
||||
return _assemble_hmn_vector(numpy.sum(transformed * m, axis=3), numpy.sum(transformed * n, axis=3))
|
||||
|
||||
return apply
|
||||
|
||||
|
||||
def find_k(
|
||||
frequency: float,
|
||||
tolerance: float,
|
||||
|
|
@ -503,20 +528,18 @@ def eigsolve(
|
|||
G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
||||
epsilon: Dielectric constant distribution for the simulation.
|
||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
||||
mu: Finite positive real diagonal permeability. Default is one.
|
||||
tolerance: Relative projected-gradient tolerance in the magnetic energy
|
||||
metric. Flat objective traces trigger a conjugate-direction restart.
|
||||
max_iters: Maximum iteration count; returns Ritz pairs from the final
|
||||
subspace if the limit is reached.
|
||||
reset_iters: Periodic conjugate-direction restart interval.
|
||||
callback: Optional no-argument callback after each subspace update.
|
||||
y0: Initial flux coefficients with shape (num_modes, 2 * grid_size).
|
||||
Copied before use.
|
||||
mu: Magnetic permability distribution for the simulation.
|
||||
Default `None` (1 everywhere).
|
||||
tolerance: Solver stops when fractional change in the objective
|
||||
`trace(Z.H @ A @ Z @ inv(Z Z.H))` is smaller than the tolerance
|
||||
max_iters: TODO
|
||||
reset_iters: TODO
|
||||
callback: TODO
|
||||
y0: TODO, initial guess
|
||||
|
||||
Returns:
|
||||
`(eigenvalues, eigenvectors)` where `eigenvalues[i]` corresponds to the
|
||||
transverse flux vector `eigenvectors[i, :]`. Pass mu to the field
|
||||
converters when it was supplied to this solve.
|
||||
vector `eigenvectors[i, :]`
|
||||
"""
|
||||
k0 = numpy.asarray(k0)
|
||||
|
||||
|
|
@ -533,24 +556,6 @@ def eigsolve(
|
|||
scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
|
||||
scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
|
||||
|
||||
metric = inverse_metric = None
|
||||
if mu is not None:
|
||||
if numpy.iscomplexobj(mu) or not numpy.isfinite(mu).all() or numpy.any(mu <= 0):
|
||||
raise ValueError('Bloch variational solving requires finite positive real permeability')
|
||||
metric = spalg.LinearOperator(
|
||||
shape=(h_size, h_size), dtype=complex,
|
||||
matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=True),
|
||||
)
|
||||
inverse_metric = spalg.LinearOperator(
|
||||
shape=(h_size, h_size), dtype=complex,
|
||||
matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=False),
|
||||
)
|
||||
|
||||
def mass(vectors: NDArray) -> NDArray:
|
||||
# P mu^-1 P makes the flux operator self-adjoint in the magnetic
|
||||
# energy inner product, even though it is not Euclidean-Hermitian.
|
||||
return vectors if metric is None else metric @ vectors
|
||||
|
||||
y_shape = (h_size, num_modes)
|
||||
|
||||
prev_E = 0.0
|
||||
|
|
@ -564,11 +569,11 @@ def eigsolve(
|
|||
if y0 is None:
|
||||
Z = rng.random(y_shape) + 1j * rng.random(y_shape)
|
||||
else:
|
||||
Z = numpy.array(y0, dtype=complex, copy=True).T
|
||||
Z = numpy.asarray(y0).T
|
||||
|
||||
while True:
|
||||
Z *= num_modes / norm(Z)
|
||||
ZtZ = Z.conj().T @ mass(Z)
|
||||
ZtZ = Z.conj().T @ Z
|
||||
try:
|
||||
U = numpy.linalg.inv(ZtZ)
|
||||
except numpy.linalg.LinAlgError:
|
||||
|
|
@ -587,28 +592,27 @@ def eigsolve(
|
|||
|
||||
for i in range(max_iters):
|
||||
Zt = numpy.conj(Z.T, out=Zt)
|
||||
ZtZ = Zt @ mass(Z)
|
||||
ZtZ = Zt @ Z
|
||||
U = numpy.linalg.inv(ZtZ)
|
||||
AZ = scipy_op @ Z.copy()
|
||||
ZtAZ = Zt @ mass(AZ)
|
||||
ZtAZ = Zt @ AZ
|
||||
ZtAZU = ZtAZ @ U
|
||||
E_signed = real(trace(ZtAZU))
|
||||
sgn = numpy.sign(E_signed)
|
||||
E = numpy.abs(E_signed)
|
||||
|
||||
# G = AZU projected onto the space orthonormal to Z via (1 - ZUZt)
|
||||
G = mass(AZ @ U - Z @ U @ ZtAZU) * sgn
|
||||
G = (AZ @ U - Z @ U @ ZtAZU) * sgn
|
||||
|
||||
if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7):
|
||||
# A flat line search alone does not imply a stationary subspace.
|
||||
# Restart with steepest descent if conjugacy has been lost.
|
||||
prev_traceGtKG = 0
|
||||
|
||||
if norm(G) <= tolerance * max(norm(mass(AZ) @ U), 1e-7):
|
||||
logger.info(
|
||||
'Optimization succeded: '
|
||||
f'[change in trace] {abs(E - prev_E)} - 5e-8 '
|
||||
f'< {tolerance} [tolerance] * {(E + prev_E) / 2} [value of trace]'
|
||||
)
|
||||
break
|
||||
# Approximate the inverse energy Hessian by a symmetric sandwich of
|
||||
# P mu P around the inverse-curl / epsilon preconditioner.
|
||||
KG = scipy_iop @ (G if inverse_metric is None else inverse_metric @ G)
|
||||
|
||||
KG = scipy_iop @ G # Preconditioned steepest descent direction
|
||||
traceGtKG = _rtrace_AtB(G, KG)
|
||||
|
||||
if prev_traceGtKG == 0 or i % reset_iters == 0:
|
||||
|
|
@ -638,11 +642,11 @@ def eigsolve(
|
|||
# Qi = inv(Q) = U'
|
||||
|
||||
AD = scipy_op @ D.copy()
|
||||
DtD = D.conj().T @ mass(D)
|
||||
DtAD = D.conj().T @ mass(AD)
|
||||
DtD = D.conj().T @ D
|
||||
DtAD = D.conj().T @ AD
|
||||
|
||||
symZtD = _symmetrize(Zt @ mass(D))
|
||||
symZtAD = _symmetrize(Zt @ mass(AD))
|
||||
symZtD = _symmetrize(Zt @ D)
|
||||
symZtAD = _symmetrize(Zt @ AD)
|
||||
|
||||
Qi_memo: list[float | None] = [None, None]
|
||||
|
||||
|
|
@ -733,10 +737,8 @@ def eigsolve(
|
|||
result = scipy.optimize.minimize_scalar(
|
||||
trace_func,
|
||||
method='bounded',
|
||||
# Center the period on zero so the line search can resolve small
|
||||
# steps without the relative angular floor near theta = pi.
|
||||
bounds=(-pi / 2, pi / 2),
|
||||
options={'xatol': min(tolerance, 1e-12)},
|
||||
bounds=(0, pi),
|
||||
options={'xatol': tolerance},
|
||||
)
|
||||
new_E = result.fun
|
||||
theta = result.x
|
||||
|
|
@ -755,9 +757,9 @@ def eigsolve(
|
|||
#
|
||||
# Recover eigenvectors from Z
|
||||
#
|
||||
U = numpy.linalg.inv(Z.conj().T @ mass(Z))
|
||||
U = numpy.linalg.inv(ZtZ)
|
||||
Y = Z @ scipy.linalg.sqrtm(U).astype(numpy.complex128)
|
||||
W = Y.conj().T @ mass(scipy_op @ Y.copy())
|
||||
W = Y.conj().T @ (scipy_op @ Y.copy())
|
||||
|
||||
eigvals, W_eigvecs = numpy.linalg.eig(W)
|
||||
eigvecs = Y @ W_eigvecs
|
||||
|
|
@ -767,10 +769,10 @@ def eigsolve(
|
|||
n = eigvals[i]
|
||||
v /= norm(v)
|
||||
Av = numpy.asarray(scipy_op @ v.copy()).reshape(-1)
|
||||
eigness = norm(Av - n * v)
|
||||
f = numpy.sqrt(max(0.0, -numpy.real(n)))
|
||||
df = numpy.sqrt(max(0.0, -numpy.real(n) + eigness))
|
||||
neff_err = kmag * (1 / df - 1 / f) if f > 0 and df > 0 else numpy.inf
|
||||
eigness = norm(Av - (v.conj() @ Av) * v)
|
||||
f = numpy.sqrt(-numpy.real(n))
|
||||
df = numpy.sqrt(-numpy.real(n) + eigness)
|
||||
neff_err = kmag * (1 / df - 1 / f)
|
||||
logger.info(f'eigness {i}: {eigness}\n neff_err: {neff_err}')
|
||||
|
||||
order = numpy.argsort(numpy.abs(eigvals))
|
||||
|
|
|
|||
|
|
@ -1,9 +1,8 @@
|
|||
"""
|
||||
Functions for performing near-to-farfield transformation (and the reverse).
|
||||
"""
|
||||
from typing import Any
|
||||
from typing import Any, cast
|
||||
from collections.abc import Sequence
|
||||
import operator
|
||||
import numpy
|
||||
from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
|
||||
from numpy import pi
|
||||
|
|
@ -14,34 +13,6 @@ type farfield_slice = NDArray[complexfloating]
|
|||
type transverse_slice_pair = Sequence[farfield_slice]
|
||||
|
||||
|
||||
def _padded_shape(shape: tuple[int, ...], padded_size: Sequence[int] | int | None) -> tuple[int, int]:
|
||||
if len(shape) != 2 or any(n == 0 for n in shape):
|
||||
raise ValueError('fields must be nonempty two-dimensional arrays')
|
||||
if padded_size is None:
|
||||
sizes = tuple(1 << (n - 1).bit_length() for n in shape)
|
||||
else:
|
||||
raw_sizes = (padded_size, padded_size) if numpy.isscalar(padded_size) else padded_size
|
||||
sizes = tuple(operator.index(n) for n in raw_sizes)
|
||||
if len(sizes) != 2 or any(n <= 0 for n in sizes):
|
||||
raise ValueError('padded_size must contain two positive integer sizes')
|
||||
return sizes[0], sizes[1]
|
||||
|
||||
|
||||
def _resize_centered(field: farfield_slice, shape: tuple[int, int]) -> farfield_slice:
|
||||
"""Pad or crop while preserving the sample at index size // 2 as the origin."""
|
||||
if field.shape == shape:
|
||||
return field
|
||||
result = numpy.zeros(shape, dtype=field.dtype)
|
||||
source, target = [], []
|
||||
for old, new in zip(field.shape, shape, strict=True):
|
||||
start_old, start_new = max(0, old // 2 - new // 2), max(0, new // 2 - old // 2)
|
||||
count = min(old - start_old, new - start_new)
|
||||
source.append(slice(start_old, start_old + count))
|
||||
target.append(slice(start_new, start_new + count))
|
||||
result[tuple(target)] = field[tuple(source)]
|
||||
return result
|
||||
|
||||
|
||||
def near_to_farfield(
|
||||
E_near: transverse_slice_pair,
|
||||
H_near: transverse_slice_pair,
|
||||
|
|
@ -53,11 +24,7 @@ def near_to_farfield(
|
|||
Compute the farfield, i.e. the distribution of the fields after propagation
|
||||
through several wavelengths of uniform medium.
|
||||
|
||||
The input fields should be complex phasors. The spatial origin is the
|
||||
sample at ``(nx // 2, ny // 2)``. Padding/cropping preserves that origin.
|
||||
Fourier integrals include the cell area ``dx * dy``.
|
||||
Phasors use ``exp(+i omega t)``; an outgoing plane wave has spatial
|
||||
dependence ``exp(-i k dot r)`` and is projected with a positive Fourier phase.
|
||||
The input fields should be complex phasors.
|
||||
|
||||
Args:
|
||||
E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
||||
|
|
@ -73,9 +40,9 @@ def near_to_farfield(
|
|||
Returns:
|
||||
Dict with keys
|
||||
|
||||
- `E`: Normalized polar and azimuthal E-field components; multiply by
|
||||
- `E_far`: Normalized E-field farfield; multiply by
|
||||
(i k exp(-i k r) / (4 pi r)) to get the actual field value.
|
||||
- `H`: Normalized polar and azimuthal H-field components; multiply by
|
||||
- `H_far`: Normalized H-field farfield; multiply by
|
||||
(i k exp(-i k r) / (4 pi r)) to get the actual field value.
|
||||
- `kx`, `ky`: Wavevector values corresponding to the x- and y- axes in E_far and H_far,
|
||||
normalized to wavelength (dimensionless).
|
||||
|
|
@ -95,10 +62,14 @@ def near_to_farfield(
|
|||
if not all(s == f.shape for f in [*E_near, *H_near]):
|
||||
raise Exception('All fields must be the same shape!')
|
||||
|
||||
padded_shape = _padded_shape(s, padded_size)
|
||||
if padded_size is None:
|
||||
padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
|
||||
if not hasattr(padded_size, '__len__'):
|
||||
padded_size = (padded_size, padded_size) # type: ignore # checked if sequence
|
||||
padded_shape = cast('Sequence[int]', padded_size)
|
||||
|
||||
En_fft = [fftshift(ifft2(ifftshift(_resize_centered(Eni, padded_shape)), norm='forward')) * dx * dy for Eni in E_near]
|
||||
Hn_fft = [fftshift(ifft2(ifftshift(_resize_centered(Hni, padded_shape)), norm='forward')) * dx * dy for Hni in H_near]
|
||||
En_fft = [fftshift(fft2(fftshift(Eni), s=padded_shape)) for Eni in E_near]
|
||||
Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_shape)) for Hni in H_near]
|
||||
|
||||
# Propagation vectors kx, ky
|
||||
k = 2 * pi
|
||||
|
|
@ -141,8 +112,8 @@ def near_to_farfield(
|
|||
outputs = {
|
||||
'E': E_far,
|
||||
'H': H_far,
|
||||
'dkx': float(2 * pi / (padded_shape[0] * dx)),
|
||||
'dky': float(2 * pi / (padded_shape[1] * dy)),
|
||||
'dkx': float(kxx[1] - kxx[0]),
|
||||
'dky': float(kyy[1] - kyy[0]),
|
||||
'kx': kx,
|
||||
'ky': ky,
|
||||
'theta': theta,
|
||||
|
|
@ -160,17 +131,18 @@ def far_to_nearfield(
|
|||
padded_size: Sequence[int] | int | None = None
|
||||
) -> dict[str, Any]:
|
||||
"""
|
||||
Reconstruct the propagating near-field spectrum from a normalized farfield.
|
||||
Compute the farfield, i.e. the distribution of the fields after propagation
|
||||
through several wavelengths of uniform medium.
|
||||
|
||||
Inputs are not modified. Spectral padding/cropping preserves zero frequency
|
||||
at index ``size // 2``. The output spacings follow its resized shape, and
|
||||
inverse Fourier quadrature restores the cell-area normalization.
|
||||
The input fields should be complex phasors.
|
||||
|
||||
Args:
|
||||
E_far: List of 2 ndarrays containing polar and azimuthal E components.
|
||||
E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
||||
E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
|
||||
Fields should be normalized so that
|
||||
E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
|
||||
H_far: List of 2 ndarrays containing polar and azimuthal H components.
|
||||
H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
||||
H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
|
||||
Fields should be normalized so that
|
||||
H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
|
||||
dkx: kx discretization, in units of wavelength.
|
||||
|
|
@ -196,7 +168,11 @@ def far_to_nearfield(
|
|||
if not all(s == f.shape for f in [*E_far, *H_far]):
|
||||
raise Exception('All fields must be the same shape!')
|
||||
|
||||
padded_shape = _padded_shape(s, padded_size)
|
||||
if padded_size is None:
|
||||
padded_size = (2 ** numpy.ceil(numpy.log2(s))).astype(int)
|
||||
if not hasattr(padded_size, '__len__'):
|
||||
padded_size = (padded_size, padded_size) # type: ignore # checked if sequence
|
||||
padded_shape = cast('Sequence[int]', padded_size)
|
||||
|
||||
k = 2 * pi
|
||||
kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
|
||||
|
|
@ -221,7 +197,9 @@ def far_to_nearfield(
|
|||
invalid_ind = kxy2 >= k * k
|
||||
theta[invalid_ind] = 0
|
||||
phi[invalid_ind] = 0
|
||||
E_far = [numpy.where(invalid_ind, 0, field) for field in E_far]
|
||||
for i in range(2):
|
||||
E_far[i][invalid_ind] = 0
|
||||
H_far[i][invalid_ind] = 0
|
||||
|
||||
# Normalized vector potentials N, L
|
||||
L = [0.5 * E_far[1],
|
||||
|
|
@ -259,11 +237,11 @@ def far_to_nearfield(
|
|||
),
|
||||
]
|
||||
|
||||
dx = 2 * pi / (padded_shape[0] * dkx)
|
||||
dy = 2 * pi / (padded_shape[1] * dky)
|
||||
E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_shape)) for Ei in En_fft]
|
||||
H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft]
|
||||
|
||||
E_near = [fftshift(fft2(ifftshift(_resize_centered(Ei, padded_shape)), norm='forward')) / (dx * dy) for Ei in En_fft]
|
||||
H_near = [fftshift(fft2(ifftshift(_resize_centered(Hi, padded_shape)), norm='forward')) / (dx * dy) for Hi in Hn_fft]
|
||||
dx = 2 * pi / (s[0] * dkx)
|
||||
dy = 2 * pi / (s[1] * dky)
|
||||
|
||||
outputs = {
|
||||
'E': E_near,
|
||||
|
|
|
|||
|
|
@ -226,8 +226,6 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi
|
|||
For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`.
|
||||
"""
|
||||
def exh(e: cfdfield, h: cfdfield) -> cfdfield_t:
|
||||
dtype = numpy.result_type(e, h, *dxes[0], *dxes[1])
|
||||
e, h = numpy.asarray(e, dtype=dtype), numpy.asarray(h, dtype=dtype)
|
||||
s = numpy.empty_like(e)
|
||||
ex = e[0] * dxes[0][0][:, None, None]
|
||||
ey = e[1] * dxes[0][1][None, :, None]
|
||||
|
|
|
|||
|
|
@ -57,8 +57,6 @@ def e_full(
|
|||
(del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
|
||||
|
||||
To make this matrix symmetric, use the preconditioners from `e_full_preconditioners()`.
|
||||
PEC rows impose E = 0 with unit diagonal entries; the corresponding
|
||||
right-hand-side entries must be zero (as done by solvers.generic).
|
||||
|
||||
Args:
|
||||
omega: Angular frequency of the simulation
|
||||
|
|
@ -95,8 +93,6 @@ def e_full(
|
|||
m_div = sparse.diags_array(1 / mu)
|
||||
|
||||
op = pe @ (ch @ pm @ m_div @ ce - omega**2 * e) @ pe
|
||||
if pec is not None:
|
||||
op += sparse.diags_array(numpy.where(pec, 1, 0))
|
||||
return op
|
||||
|
||||
|
||||
|
|
@ -141,9 +137,6 @@ def h_full(
|
|||
|
||||
del x (1/epsilon * del x) - omega**2 * mu
|
||||
|
||||
PMC rows impose H = 0 with unit diagonal entries. Their right-hand-side
|
||||
entries must be zero.
|
||||
|
||||
for use with the H-field, with wave equation
|
||||
$$ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M $$
|
||||
|
||||
|
|
@ -184,8 +177,6 @@ def h_full(
|
|||
m = sparse.diags_array(mu)
|
||||
|
||||
A = pm @ (ce @ pe @ e_div @ ch - omega**2 * m) @ pm
|
||||
if pmc is not None:
|
||||
A += sparse.diags_array(numpy.where(pmc, 1, 0))
|
||||
return A
|
||||
|
||||
|
||||
|
|
@ -208,9 +199,6 @@ def eh_full(
|
|||
[[-i * omega * epsilon, del x ],
|
||||
[del x, i * omega * mu]]
|
||||
|
||||
PEC electric rows and PMC magnetic rows impose zero fields with unit
|
||||
diagonal entries. The corresponding source entries must be zero.
|
||||
|
||||
for use with a field vector of the form `cat(vec(E), vec(H))`:
|
||||
$$ \begin{bmatrix}
|
||||
-\imath \omega \epsilon & \nabla \times \\
|
||||
|
|
@ -259,10 +247,8 @@ def eh_full(
|
|||
A1 = pe @ curl_back(dxes[1]) @ pm
|
||||
A2 = pm @ curl_forward(dxes[0]) @ pe
|
||||
|
||||
electric = -iwe if pec is None else -iwe + sparse.diags_array(numpy.where(pec, 1, 0))
|
||||
magnetic = iwm if pmc is None else iwm + sparse.diags_array(numpy.where(pmc, 1, 0))
|
||||
A = sparse.block_array([[electric, A1],
|
||||
[A2, magnetic]])
|
||||
A = sparse.block_array([[-iwe, A1],
|
||||
[A2, iwm]])
|
||||
return A
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ def _scipy_qmr(
|
|||
kwargs: Passed as **kwargs to the wrapped function
|
||||
|
||||
Returns:
|
||||
Converged solution. Raises RuntimeError on breakdown or nonconvergence.
|
||||
Guess for solution (returned even if didn't converge)
|
||||
"""
|
||||
|
||||
#
|
||||
|
|
@ -61,10 +61,7 @@ def _scipy_qmr(
|
|||
#
|
||||
# Run the actual solve
|
||||
#
|
||||
x, info = scipy.sparse.linalg.qmr(A, b, **kwargs)
|
||||
if info != 0:
|
||||
residual = norm(A @ x - b)
|
||||
raise RuntimeError(f'QMR failed to converge ({info=}, residual={residual:.6g})')
|
||||
x, _ = scipy.sparse.linalg.qmr(A, b, **kwargs)
|
||||
return x
|
||||
|
||||
|
||||
|
|
@ -83,7 +80,7 @@ def generic(
|
|||
E_guess: vcfdfield | None = None,
|
||||
) -> vcfdfield_t:
|
||||
"""
|
||||
Preconditioned FDFD solver using CSR sparse matrices.
|
||||
Conjugate gradient FDFD solver using CSR sparse matrices.
|
||||
|
||||
All ndarray arguments should be 1D arrays, as returned by `meanas.fdmath.vectorization.vec()`.
|
||||
|
||||
|
|
@ -95,8 +92,7 @@ def generic(
|
|||
epsilon: Dielectric constant distribution (at E-field locations)
|
||||
mu: Magnetic permeability distribution (at H-field locations)
|
||||
pec: Perfect electric conductor distribution
|
||||
(at E-field locations; non-zero value indicates PEC is present).
|
||||
Sources and initial guesses are zeroed at PEC locations.
|
||||
(at E-field locations; non-zero value indicates PEC is present)
|
||||
pmc: Perfect magnetic conductor distribution
|
||||
(at H-field locations; non-zero value indicates PMC is present)
|
||||
adjoint: If true, solves the adjoint problem.
|
||||
|
|
@ -105,8 +101,8 @@ def generic(
|
|||
`b`: `ArrayLike`;
|
||||
`x`: `ArrayLike`;
|
||||
Default is a wrapped version of `scipy.sparse.linalg.qmr()`
|
||||
which raises on nonconvergence and logs the residual every
|
||||
100 iterations.
|
||||
which doesn't return convergence info and logs the residual
|
||||
every 100 iterations.
|
||||
matrix_solver_opts: Passed as kwargs to `matrix_solver(...)`
|
||||
E_guess: Guess at the solution E-field. `matrix_solver` must accept an
|
||||
`x0` argument with the same purpose.
|
||||
|
|
@ -115,11 +111,10 @@ def generic(
|
|||
E-field which solves the system.
|
||||
"""
|
||||
|
||||
matrix_solver_opts = {} if matrix_solver_opts is None else dict(matrix_solver_opts)
|
||||
if matrix_solver_opts is None:
|
||||
matrix_solver_opts = dict()
|
||||
|
||||
b0 = -1j * omega * J
|
||||
if pec is not None:
|
||||
b0 = numpy.where(pec, 0, b0)
|
||||
A0 = operators.e_full(omega, dxes, epsilon=epsilon, mu=mu, pec=pec, pmc=pmc)
|
||||
|
||||
Pl, Pr = operators.e_full_preconditioners(dxes)
|
||||
|
|
@ -132,23 +127,17 @@ def generic(
|
|||
b = Pl @ b0
|
||||
|
||||
if E_guess is not None:
|
||||
if pec is not None:
|
||||
E_guess = numpy.where(pec, 0, E_guess)
|
||||
if adjoint:
|
||||
x0 = Pr.T.conjugate() @ E_guess
|
||||
else:
|
||||
x0 = Pl @ E_guess
|
||||
matrix_solver_opts['x0'] = x0
|
||||
|
||||
x = numpy.asarray(matrix_solver(A.tocsr(), b, **matrix_solver_opts))
|
||||
if x.shape != b.shape or not numpy.isfinite(x).all():
|
||||
raise RuntimeError('matrix_solver must return a finite solution vector matching the source shape')
|
||||
x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
|
||||
|
||||
if adjoint:
|
||||
x0 = Pl.T.conjugate() @ x
|
||||
else:
|
||||
x0 = Pr @ x
|
||||
|
||||
if pec is not None:
|
||||
x0 = numpy.where(pec, 0, x0)
|
||||
return vcfdfield_t(x0)
|
||||
|
|
|
|||
|
|
@ -801,40 +801,74 @@ def sensitivity(
|
|||
wavenumber: complex,
|
||||
omega: complex,
|
||||
dxes: dx_lists2_t,
|
||||
epsilon: vfdslice | vcfdslice,
|
||||
mu: vfdslice | vcfdslice | None = None,
|
||||
epsilon: vfdslice,
|
||||
mu: vfdslice | None = None,
|
||||
) -> vcfdslice_t:
|
||||
r"""
|
||||
Calculate the analytic sensitivity of an isolated mode's beta to epsilon.
|
||||
Given a waveguide structure (`dxes`, `epsilon`, `mu`) and mode fields
|
||||
(`e_norm`, `h_norm`, `wavenumber`, `omega`), calculates the sensitivity of the wavenumber
|
||||
$\beta$ to changes in the dielectric structure $\epsilon$.
|
||||
|
||||
For the reciprocal diagonal waveguide problem ``A_E E_xy = beta**2 E_xy``,
|
||||
the Yee-weighted Lorentz left eigenvector is the unconjugated row
|
||||
The output is a vector of the same size as `vec(epsilon)`, with each element specifying the
|
||||
sensitivity of `wavenumber` to changes in the corresponding element in `vec(epsilon)`, i.e.
|
||||
|
||||
$$ L^T = [H_y \Delta x_e \Delta y_h,\; -H_x \Delta x_h \Delta y_e]. $$
|
||||
$$sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}$$
|
||||
|
||||
Differentiating the eigenproblem and projecting with this left vector gives
|
||||
An adjoint approach is used to calculate the sensitivity; the derivation is provided here:
|
||||
|
||||
$$ \frac{\partial\beta}{\partial\epsilon_i}
|
||||
= \frac{L^T (\partial_{\epsilon_i} A_E) E_{xy}}
|
||||
{2\beta L^T E_{xy}}. $$
|
||||
Starting with the eigenvalue equation
|
||||
|
||||
$$\beta^2 E_{xy} = A_E E_{xy}$$
|
||||
|
||||
where $A_E$ is the waveguide operator from `operator_e()`, and $E_{xy} = \begin{bmatrix} E_x \\
|
||||
E_y \end{bmatrix}$,
|
||||
we can differentiate with respect to one of the $\epsilon$ elements (i.e. at one Yee grid point), $\epsilon_i$:
|
||||
|
||||
$$
|
||||
(2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy}
|
||||
= \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy}
|
||||
$$
|
||||
|
||||
We then multiply by $H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}$ from the left:
|
||||
|
||||
$$
|
||||
(2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
|
||||
= H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy}
|
||||
$$
|
||||
|
||||
However, $H_{yx}^\star$ is actually a left-eigenvector of $A_E$. This can be verified by inspecting
|
||||
the form of `operator_h` ($A_H$) and comparing its conjugate transpose to `operator_e` ($A_E$). Also, note
|
||||
$H_{yx}^\star \cdot E_{xy} = H^\star \times E$ recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201
|
||||
for a similar approach. Therefore,
|
||||
|
||||
$$
|
||||
H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
|
||||
$$
|
||||
|
||||
and we can simplify to
|
||||
|
||||
$$
|
||||
\partial_{\epsilon_i}(\beta)
|
||||
= \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}}
|
||||
$$
|
||||
|
||||
This expression can be quickly calculated for all $i$ by writing out the various terms of
|
||||
$\partial_{\epsilon_i} A_E$ and recognizing that the vector-matrix-vector products (i.e. scalars)
|
||||
$sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}$, indexed by $i$, can be expressed as
|
||||
elementwise multiplications $\vec{sens} = \vec{v}_{left} \star \vec{v}_{right}$
|
||||
|
||||
This pairing is bilinear: neither H nor the coordinate weights are
|
||||
conjugated. It therefore applies to complex PML stretches and reciprocal
|
||||
lossy materials as well as lossless grids. Quadrature weights multiply the
|
||||
left projection before derivative operators act. Overall field scaling
|
||||
cancels between numerator and denominator.
|
||||
|
||||
Args:
|
||||
e_norm: Vectorized E_xyz eigenfield for the mode.
|
||||
h_norm: Matching vectorized H_xyz eigenfield at the same phase reference.
|
||||
wavenumber: Nonzero beta of the continuous-z eigenproblem.
|
||||
omega: Angular frequency of the system.
|
||||
dxes: Two-dimensional Yee widths, including any complex PML stretches.
|
||||
epsilon: Vectorized diagonal permittivity.
|
||||
mu: Optional vectorized diagonal permeability; defaults to one.
|
||||
e_norm: Normalized, vectorized E_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
|
||||
h_norm: Normalized, vectorized H_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
|
||||
wavenumber: Propagation constant for the mode. The z-axis is assumed to be continuous (i.e. without numerical dispersion).
|
||||
omega: The angular frequency of the system.
|
||||
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
||||
epsilon: Vectorized dielectric constant grid
|
||||
mu: Vectorized magnetic permeability grid (default 1 everywhere)
|
||||
|
||||
Returns:
|
||||
Vector ``d beta / d epsilon_i`` in the same component order as epsilon.
|
||||
Sparse matrix representation of the operator.
|
||||
"""
|
||||
if mu is None:
|
||||
mu = numpy.ones_like(epsilon)
|
||||
|
|
@ -849,17 +883,16 @@ def sensitivity(
|
|||
mu_x, mu_y, _mu_z = numpy.split(mu, 3)
|
||||
mu_yx = sparse.diags_array(numpy.hstack((mu_y, mu_x)))
|
||||
|
||||
da_exxhyy = vec(dxes[0][0][:, None] * dxes[1][1][None, :])
|
||||
da_eyyhxx = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
|
||||
ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2])
|
||||
hx, hy, _hz = numpy.split(h_norm, 3)
|
||||
# The Yee-weighted H projection is the left eigenvector on a graded grid.
|
||||
hv_yx = numpy.concatenate([hy, -hx]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
|
||||
da_exxhyy = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
|
||||
da_eyyhxx = vec(dxes[1][1][None, :] * dxes[0][0][:, None])
|
||||
ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
|
||||
hx, hy, hz = numpy.split(h_norm, 3)
|
||||
hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx]))
|
||||
|
||||
sens_xy1 = (hv_yx @ (omega * omega * mu_yx)) * ev_xy
|
||||
sens_xy2 = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
|
||||
sens_z = (hv_yx @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy)
|
||||
norm = hv_yx @ ev_xy
|
||||
sens_xy1 = (hv_yx_conj @ (omega * omega * mu_yx)) * ev_xy
|
||||
sens_xy2 = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
|
||||
sens_z = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy)
|
||||
norm = hv_yx_conj @ ev_xy
|
||||
|
||||
sens_tot = numpy.concatenate([sens_xy1 + sens_xy2, sens_z]) / (2 * wavenumber * norm)
|
||||
return vcfdslice_t(sens_tot)
|
||||
|
|
|
|||
|
|
@ -143,9 +143,6 @@ def mode_from_2d(
|
|||
|
||||
``polarity`` is +1 toward increasing indices and -1 toward decreasing
|
||||
indices. ``normalization`` follows waveguide_2d.normalized_fields_e.
|
||||
Lossless modes at or beyond the longitudinal grid band edge cannot be
|
||||
power-normalized. Lorentz normalization can represent these modes
|
||||
beyond that edge and selects the decaying branch for real beta.
|
||||
"""
|
||||
if polarity not in (-1, 1):
|
||||
raise ValueError('polarity must be -1 or +1')
|
||||
|
|
@ -157,21 +154,15 @@ def mode_from_2d(
|
|||
raise ValueError('the 2D eigenvector and wavenumber must be finite and match the cross-section')
|
||||
wavenumber_2d = complex(wavenumber_2d)
|
||||
# Correct wavenumber to account for numerical dispersion.
|
||||
argument = wavenumber_2d * dx_prop / 2
|
||||
if numpy.isreal(argument) and abs(argument.real) >= 1:
|
||||
if normalization == 'power':
|
||||
raise ValueError('mode is at or beyond the longitudinal grid band edge; refine the propagation grid')
|
||||
argument = complex(argument.real, -0.0)
|
||||
wavenumber = 2 / dx_prop * numpy.arcsin(argument)
|
||||
wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2)
|
||||
|
||||
shape = [d.size for d in args_2d['dxes'][0]]
|
||||
ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, prop_phase=dx_prop * wavenumber, normalization=normalization, **args_2d)
|
||||
e = unvec(ve, shape)
|
||||
h = unvec(vh, shape)
|
||||
|
||||
# Reversing beta changes Ez and transverse H, while Hz is unchanged.
|
||||
e[2] *= polarity
|
||||
h[:2] *= polarity
|
||||
# Adjust for propagation direction
|
||||
h *= polarity # type: ignore # mypy issue with numpy
|
||||
|
||||
# Apply phase shift to H-field
|
||||
h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
|
||||
|
|
@ -380,18 +371,12 @@ def expand_e(
|
|||
$$
|
||||
|
||||
to each copied slice.
|
||||
|
||||
Transverse electric samples use forward intervals dx_e[i] between
|
||||
indices i and i+1. The longitudinal electric component uses backward
|
||||
intervals dx_h[i] between indices i-1 and i on its staggered grid.
|
||||
"""
|
||||
slices = tuple(slices)
|
||||
|
||||
# Determine phase factors for parallel slices
|
||||
a_shape = numpy.roll([1, -1, 1, 1], axis)
|
||||
# dx_e[i] spans transverse electric samples i -> i+1. The longitudinal
|
||||
# component lies on the complementary grid, with backward intervals dx_h.
|
||||
a_E = numpy.concatenate(([0.0], numpy.real(dxes[0][axis][:-1]).cumsum()))
|
||||
a_E = numpy.real(dxes[0][axis]).cumsum()
|
||||
r_E = a_E - a_E[slices[axis]]
|
||||
iphi = polarity * -1j * wavenumber
|
||||
phase_E = numpy.exp(iphi * r_E).reshape(a_shape)
|
||||
|
|
@ -406,9 +391,6 @@ def expand_e(
|
|||
slices_in = (slice(None), *slices)
|
||||
|
||||
E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
|
||||
a_long = numpy.real(dxes[1][axis]).cumsum()
|
||||
phase_long = numpy.exp(iphi * (a_long - a_long[slices[axis]])).reshape(a_shape[1:])
|
||||
E_expanded[(axis, *slices_exp_l)] = phase_long * E[(axis, *slices)]
|
||||
return cfdfield_t(E_expanded)
|
||||
|
||||
|
||||
|
|
@ -567,41 +549,10 @@ def _decompose_port_field(
|
|||
matrix, overlaps, mode_count = decomposition
|
||||
field_vec = vec(field)
|
||||
measured = numpy.asarray([field_vec @ vec(overlap).conj() for overlap in overlaps])
|
||||
coeffs, _resid, rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
|
||||
if rank != matrix.shape[1]:
|
||||
raise numpy.linalg.LinAlgError('monitor modal basis is rank deficient')
|
||||
coeffs, _resid, _rank, _sing = numpy.linalg.lstsq(matrix, measured, rcond=None)
|
||||
return coeffs[:mode_count], coeffs[mode_count:]
|
||||
|
||||
|
||||
def _port_power_gram(
|
||||
modes: Sequence[Waveguide3DMode],
|
||||
dxes: dx_lists_t,
|
||||
axis: int,
|
||||
polarity: int,
|
||||
slices: Sequence[slice],
|
||||
) -> NDArray[numpy.complex128]:
|
||||
"""Hermitian flux metric, including interference between nonorthogonal modes."""
|
||||
u, v = (axis + 1) % 3, (axis + 2) % 3
|
||||
|
||||
def width(family: int, component: int) -> NDArray:
|
||||
values = numpy.real(dxes[family][component][slices[component]])
|
||||
shape = [1, 1, 1]
|
||||
shape[component] = values.size
|
||||
return values.reshape(shape)
|
||||
|
||||
area_uv = width(0, u) * width(1, v)
|
||||
area_vu = width(1, u) * width(0, v)
|
||||
dx = float(numpy.real(dxes[0][axis][slices[axis]]).item())
|
||||
gram = numpy.empty((len(modes), len(modes)), dtype=complex)
|
||||
for row, test in enumerate(modes):
|
||||
h = test['H'][(slice(None), *slices)].conj()
|
||||
for column, mode in enumerate(modes):
|
||||
e = mode['E'][(slice(None), *slices)]
|
||||
phase = numpy.exp(-1j * polarity * mode['wavenumber'] * dx)
|
||||
gram[row, column] = 0.5 * polarity * phase * numpy.sum(e[u] * h[v] * area_uv - e[v] * h[u] * area_vu)
|
||||
return 0.5 * (gram + gram.conj().T)
|
||||
|
||||
|
||||
def _copy_cross_section(
|
||||
target: NDArray,
|
||||
cross_section: NDArray,
|
||||
|
|
@ -685,11 +636,8 @@ def get_s_fdfd(
|
|||
Extract a two-port guided-mode scattering matrix from local 3D FDFD solves.
|
||||
|
||||
The returned S-matrix is ordered as `[[R_L, T_R_to_L], [T_L_to_R, R_R]]`, with
|
||||
reference planes at the supplied monitor slices. Calibration solves
|
||||
``S @ incoming_amplitudes = raw_outgoing_amplitudes`` using every measured
|
||||
incoming coefficient on both ports, including mixed or degenerate modes.
|
||||
Diagnostic guided powers use the full input/output flux Gram matrices,
|
||||
retaining interference between nonorthogonal modal basis vectors.
|
||||
reference planes at the supplied monitor slices. Each column is normalized
|
||||
by the incident modal amplitude measured on the corresponding input monitor.
|
||||
"""
|
||||
if incident_amplitude_atol <= 0:
|
||||
raise ValueError('incident_amplitude_atol must be positive')
|
||||
|
|
@ -761,8 +709,8 @@ def get_s_fdfd(
|
|||
slices=right_monitor,
|
||||
)
|
||||
|
||||
raw_outgoing = numpy.zeros((n_total, n_total), dtype=complex)
|
||||
incoming_amplitudes = numpy.zeros_like(raw_outgoing)
|
||||
s = numpy.zeros((n_total, n_total), dtype=complex)
|
||||
raw_outgoing = numpy.zeros_like(s)
|
||||
incident_amplitudes = numpy.zeros((n_total,), dtype=complex)
|
||||
electric_fields: list[cfdfield_t] = []
|
||||
magnetic_fields: list[cfdfield_t] = []
|
||||
|
|
@ -798,39 +746,26 @@ def get_s_fdfd(
|
|||
e_field = unvec(solvers.generic(**solve_kwargs), shape)
|
||||
left_in, left_out = _decompose_port_field(left_decomposition, e_field)
|
||||
right_in, right_out = _decompose_port_field(right_decomposition, e_field)
|
||||
incoming_column = numpy.concatenate((left_in, right_in))
|
||||
if numpy.linalg.norm(incoming_column) <= incident_amplitude_atol:
|
||||
if column < n_left:
|
||||
incident = left_in[column]
|
||||
else:
|
||||
incident = right_in[column - n_left]
|
||||
if abs(incident) <= incident_amplitude_atol:
|
||||
raise ValueError(f'incident modal amplitude for column {column} is zero')
|
||||
|
||||
raw_column = numpy.concatenate((left_out, right_out))
|
||||
|
||||
incident_amplitudes[column] = incoming_column[column]
|
||||
incoming_amplitudes[:, column] = incoming_column
|
||||
incident_amplitudes[column] = incident
|
||||
raw_outgoing[:, column] = raw_column
|
||||
s[:, column] = raw_column / incident
|
||||
if return_diagnostics:
|
||||
electric_fields.append(cfdfield_t(e_field))
|
||||
magnetic_fields.append(cfdfield_t(functional.e2h(omega, solve_dxes, mu=mu)(e_field)))
|
||||
|
||||
if not numpy.isfinite(incoming_amplitudes).all() or not numpy.isfinite(raw_outgoing).all():
|
||||
raise ValueError('modal calibration requires finite incoming and outgoing amplitudes')
|
||||
try:
|
||||
s = numpy.linalg.solve(incoming_amplitudes.T, raw_outgoing.T).T
|
||||
except numpy.linalg.LinAlgError as error:
|
||||
raise numpy.linalg.LinAlgError('incoming excitation matrix is singular') from error
|
||||
|
||||
if not return_diagnostics:
|
||||
return s
|
||||
|
||||
zero = numpy.zeros((n_left, n_right), dtype=complex)
|
||||
input_power_gram = numpy.block([
|
||||
[_port_power_gram(left_in_modes, dxes, axis, +1, left_monitor), zero],
|
||||
[zero.T, _port_power_gram(right_in_modes, dxes, axis, -1, right_monitor)],
|
||||
])
|
||||
output_power_gram = numpy.block([
|
||||
[_port_power_gram(left_out_modes, dxes, axis, -1, left_monitor), zero],
|
||||
[zero.T, _port_power_gram(right_out_modes, dxes, axis, +1, right_monitor)],
|
||||
])
|
||||
guided_powers = numpy.diag(s.conj().T @ output_power_gram @ s).real / numpy.diag(input_power_gram).real
|
||||
guided_powers = numpy.sum(numpy.abs(s) ** 2, axis=0).real
|
||||
diagnostics: dict[str, Any] = {
|
||||
'axis': axis,
|
||||
'left_mode_numbers': left_modes_i,
|
||||
|
|
@ -845,9 +780,6 @@ def get_s_fdfd(
|
|||
'electric_fields': tuple(electric_fields),
|
||||
'magnetic_fields': tuple(magnetic_fields),
|
||||
'incident_amplitudes': incident_amplitudes,
|
||||
'incoming_amplitudes': incoming_amplitudes,
|
||||
'incoming_power_gram': input_power_gram,
|
||||
'outgoing_power_gram': output_power_gram,
|
||||
'raw_outgoing_amplitudes': raw_outgoing,
|
||||
'guided_powers': guided_powers,
|
||||
'unresolved_powers': 1.0 - guided_powers,
|
||||
|
|
|
|||
|
|
@ -86,7 +86,7 @@ large-radius-compatible identities
|
|||
|
||||
$$
|
||||
E_z =
|
||||
\frac{1}{\imath \beta} \epsilon_{zz}^{-1}
|
||||
\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1}
|
||||
\begin{bmatrix}
|
||||
\hat{\partial}_r T_b \epsilon_{rr} &
|
||||
\hat{\partial}_y T_a \epsilon_{yy}
|
||||
|
|
@ -304,11 +304,6 @@ def linear_wavenumbers(
|
|||
Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad)
|
||||
and the mode's energy distribution.
|
||||
|
||||
The effective radius is weighted by transverse electric energy, using the
|
||||
physical Yee cell areas for each component. Radial E uses dx_e * dy_h;
|
||||
transverse E uses dx_h * dy_e. Real coordinate widths define these areas
|
||||
and radii when complex PML stretching is present.
|
||||
|
||||
Args:
|
||||
e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y)
|
||||
angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of
|
||||
|
|
@ -326,19 +321,14 @@ def linear_wavenumbers(
|
|||
|
||||
shape2d = (len(dxes[0][0]), len(dxes[0][1]))
|
||||
epsilon2d = unvec(epsilon, shape2d)[:2]
|
||||
physical_widths = [[numpy.real(dx) for dx in family] for family in dxes]
|
||||
ra = rmin + numpy.cumsum(physical_widths[0][0])
|
||||
rb = rmin + physical_widths[0][0] / 2.0 + numpy.concatenate((
|
||||
numpy.zeros(1),
|
||||
numpy.cumsum(physical_widths[1][0][:-1]),
|
||||
))
|
||||
areas = numpy.stack((
|
||||
numpy.outer(physical_widths[0][0], physical_widths[1][1]),
|
||||
numpy.outer(physical_widths[1][0], physical_widths[0][1]),
|
||||
ra = rmin + numpy.cumsum(dxes[0][0])
|
||||
rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate((
|
||||
numpy.zeros(1, dtype=dxes[1][0].dtype),
|
||||
numpy.cumsum(dxes[1][0][:-1]),
|
||||
))
|
||||
for ii in range(angular_wavenumbers.size):
|
||||
efield = unvec(e_xys[ii], shape2d, 2)
|
||||
energy = numpy.real((efield * efield.conj()) * epsilon2d) * areas
|
||||
energy = numpy.real((efield * efield.conj()) * epsilon2d)
|
||||
er_energy_vs_r = energy[0].sum(axis=1)
|
||||
ey_energy_vs_r = energy[1].sum(axis=1)
|
||||
energy_vs_r = er_energy_vs_r + ey_energy_vs_r
|
||||
|
|
@ -415,7 +405,7 @@ def exy2e(
|
|||
|
||||
n_pts = dxes[0][0].size * dxes[0][1].size
|
||||
exy2ez = (
|
||||
epsilon_z_inv
|
||||
Ta @ epsilon_z_inv
|
||||
@ sparse.hstack((Dbx @ Tb @ epsilon_x,
|
||||
Dby @ Ta @ epsilon_y))
|
||||
/ (1j * wavenumber)
|
||||
|
|
|
|||
|
|
@ -71,12 +71,6 @@ def deriv_back(
|
|||
return tuple(_derivative(axis, False, None if dx_h is None else dx_h[axis]) for axis in range(3))
|
||||
|
||||
|
||||
def _promote_field(field: NDArray, widths: Sequence[NDArray] | None) -> NDArray:
|
||||
"""Use one calculation dtype for all components of a curl."""
|
||||
dtype = field.dtype if widths is None else numpy.result_type(field, *widths)
|
||||
return numpy.asarray(field, dtype=dtype)
|
||||
|
||||
|
||||
def curl_forward(
|
||||
dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
|
||||
) -> Callable[[NDArray[floating | complexfloating]], NDArray[floating | complexfloating]]:
|
||||
|
|
@ -95,8 +89,8 @@ def curl_forward(
|
|||
Dx, Dy, Dz = deriv_forward(dx_e)
|
||||
|
||||
def ce_fun(e: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]:
|
||||
e = _promote_field(e, dx_e)
|
||||
output = numpy.empty_like(e)
|
||||
dtype = e.dtype if dx_e is None else numpy.result_type(e, *dx_e)
|
||||
output = numpy.empty_like(e, dtype=dtype)
|
||||
output[0] = Dy(e[2])
|
||||
output[1] = Dz(e[0])
|
||||
output[2] = Dx(e[1])
|
||||
|
|
@ -126,8 +120,8 @@ def curl_back(
|
|||
Dx, Dy, Dz = deriv_back(dx_h)
|
||||
|
||||
def ch_fun(h: NDArray[floating | complexfloating]) -> NDArray[floating | complexfloating]:
|
||||
h = _promote_field(h, dx_h)
|
||||
output = numpy.empty_like(h)
|
||||
dtype = h.dtype if dx_h is None else numpy.result_type(h, *dx_h)
|
||||
output = numpy.empty_like(h, dtype=dtype)
|
||||
output[0] = Dy(h[2])
|
||||
output[1] = Dz(h[0])
|
||||
output[2] = Dx(h[1])
|
||||
|
|
@ -145,7 +139,6 @@ def curl_forward_parts(
|
|||
Dx, Dy, Dz = deriv_forward(dx_e)
|
||||
|
||||
def mkparts_fwd(e: fdfield) -> tuple[tuple[fdfield, fdfield], ...]:
|
||||
e = _promote_field(e, dx_e)
|
||||
return ((-Dz(e[1]), Dy(e[2])),
|
||||
( Dz(e[0]), -Dx(e[2])),
|
||||
(-Dy(e[0]), Dx(e[1])))
|
||||
|
|
@ -159,7 +152,6 @@ def curl_back_parts(
|
|||
Dx, Dy, Dz = deriv_back(dx_h)
|
||||
|
||||
def mkparts_back(h: fdfield) -> tuple[tuple[fdfield, fdfield], ...]:
|
||||
h = _promote_field(h, dx_h)
|
||||
return ((-Dz(h[1]), Dy(h[2])),
|
||||
( Dz(h[0]), -Dx(h[2])),
|
||||
(-Dy(h[0]), Dx(h[1])))
|
||||
|
|
|
|||
|
|
@ -120,7 +120,7 @@ def deriv_forward(
|
|||
def deriv(axis: int) -> sparse.sparray:
|
||||
return shift_circ(axis, shape, 1) - sparse.eye_array(n)
|
||||
|
||||
Ds = [sparse.diags_array(numpy.reciprocal(dx.ravel(order='C'), dtype=numpy.result_type(dx, numpy.float64))) @ deriv(a)
|
||||
Ds = [sparse.diags_array(+1 / dx.ravel(order='C')) @ deriv(a)
|
||||
for a, dx in enumerate(dx_e_expanded)]
|
||||
|
||||
return Ds
|
||||
|
|
@ -147,7 +147,7 @@ def deriv_back(
|
|||
def deriv(axis: int) -> sparse.sparray:
|
||||
return shift_circ(axis, shape, -1) - sparse.eye_array(n)
|
||||
|
||||
Ds = [sparse.diags_array(-numpy.reciprocal(dx.ravel(order='C'), dtype=numpy.result_type(dx, numpy.float64))) @ deriv(a)
|
||||
Ds = [sparse.diags_array(-1 / dx.ravel(order='C')) @ deriv(a)
|
||||
for a, dx in enumerate(dx_h_expanded)]
|
||||
|
||||
return Ds
|
||||
|
|
|
|||
|
|
@ -1,61 +1,10 @@
|
|||
import numpy
|
||||
import pytest
|
||||
from numpy.linalg import norm
|
||||
|
||||
from ..fdfd import bloch
|
||||
from ._bloch_case import EPSILON, G_MATRIX, H_MN, K0_AXIAL, K0_GENERAL, MU, SHAPE, ZERO_H_MN
|
||||
from .utils import assert_close
|
||||
|
||||
|
||||
@pytest.mark.parametrize('k0', [[0.2, 0, 0], [0, 0.2, 0], [0, 0, 0.2], [0, 0.2, 0.1], [0, 0, 0]])
|
||||
def test_kmn_basis_is_transverse_on_axes_and_at_gamma(k0) -> None:
|
||||
k0 = numpy.asarray(k0)
|
||||
shape = (2, 3, 2)
|
||||
_magnitude, m, n = bloch.generate_kmn(k0, numpy.eye(3), shape)
|
||||
reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid(
|
||||
*(numpy.fft.fftfreq(size, 1 / size) for size in shape), indexing='ij',
|
||||
)), 0, -1)
|
||||
k = k0 - reciprocal_indices
|
||||
assert_close(numpy.sum(k * m, axis=-1), 0, atol=1e-14)
|
||||
assert_close(numpy.sum(k * n, axis=-1), 0, atol=1e-14)
|
||||
assert_close(numpy.sum(m * n, axis=-1), 0, atol=1e-14)
|
||||
assert_close(norm(m, axis=-1), 1)
|
||||
assert_close(norm(n, axis=-1), 1)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('mu', [None, MU])
|
||||
def test_inverse_preconditioner_is_finite_at_gamma(mu) -> None:
|
||||
operator = bloch.inverse_maxwell_operator_approx(numpy.zeros(3), G_MATRIX, EPSILON, mu)
|
||||
with numpy.errstate(divide='raise', invalid='raise', over='raise'):
|
||||
result = operator(H_MN)
|
||||
zero = operator(ZERO_H_MN)
|
||||
assert numpy.isfinite(result).all()
|
||||
assert_close(zero, 0)
|
||||
|
||||
|
||||
def test_anisotropic_magnetic_bloch_modes_obey_full_fourier_maxwell_equations() -> None:
|
||||
shape = EPSILON.shape[1:]
|
||||
size = 2 * numpy.prod(shape)
|
||||
operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
|
||||
matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(size, dtype=complex)])
|
||||
values, vectors = numpy.linalg.eig(matrix)
|
||||
reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid(
|
||||
*(numpy.fft.fftfreq(n, 1 / n) for n in shape), indexing='ij',
|
||||
)), 0, -1)
|
||||
k = K0_GENERAL - reciprocal_indices
|
||||
epsilon, mu = numpy.moveaxis(EPSILON, 0, -1), numpy.moveaxis(MU, 0, -1)
|
||||
h_converter = bloch.hmn_2_hxyz(K0_GENERAL, G_MATRIX, EPSILON, MU)
|
||||
e_converter = bloch.hmn_2_exyz(K0_GENERAL, G_MATRIX, EPSILON, MU)
|
||||
for value, vector in zip(values, vectors.T, strict=True):
|
||||
h = numpy.moveaxis(h_converter(vector), 0, -1)
|
||||
b = mu * h
|
||||
h_fourier = numpy.fft.fftn(h, axes=(0, 1, 2))
|
||||
electric = numpy.fft.ifftn(numpy.cross(k, h_fourier), axes=(0, 1, 2)) / epsilon
|
||||
curl_electric = numpy.fft.ifftn(numpy.cross(k, numpy.fft.fftn(electric, axes=(0, 1, 2))), axes=(0, 1, 2))
|
||||
assert_close(curl_electric, value * b, atol=1e-13, rtol=1e-11)
|
||||
assert_close(numpy.sum(k * numpy.fft.fftn(b, axes=(0, 1, 2)), axis=-1), 0, atol=1e-13)
|
||||
assert_close(numpy.moveaxis(e_converter(vector), 0, -1), electric, atol=1e-13)
|
||||
|
||||
def test_generate_kmn_general_case_returns_orthonormal_basis() -> None:
|
||||
k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_GENERAL, G_MATRIX, SHAPE)
|
||||
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ from numpy.testing import assert_allclose
|
|||
from types import SimpleNamespace
|
||||
|
||||
from ..fdfd import bloch
|
||||
from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, K0_GENERAL, MU, Y0, Y0_TWO_MODE, build_overlap_fixture
|
||||
from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, Y0, Y0_TWO_MODE, build_overlap_fixture
|
||||
from .utils import assert_close
|
||||
|
||||
|
||||
|
|
@ -81,42 +81,6 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
|
|||
assert callback_count > 0
|
||||
|
||||
|
||||
@pytest.mark.parametrize('count', [1, 2])
|
||||
def test_iteration_limit_returns_consistent_ritz_pairs_without_mutating_guess(count) -> None:
|
||||
rng = numpy.random.default_rng(11)
|
||||
guess = rng.normal(size=(count, H_SIZE)) + 1j * rng.normal(size=(count, H_SIZE))
|
||||
original = guess.copy()
|
||||
values, vectors = bloch.eigsolve(count, K0_X, G_MATRIX, EPSILON, y0=guess, max_iters=1)
|
||||
operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
|
||||
for value, vector in zip(values, vectors, strict=True):
|
||||
rayleigh = numpy.vdot(vector, operator(vector).ravel()) / numpy.vdot(vector, vector)
|
||||
assert_close(value, rayleigh, atol=1e-13)
|
||||
assert_allclose(vectors.conj() @ vectors.T, numpy.eye(count), atol=1e-12)
|
||||
numpy.testing.assert_array_equal(guess, original)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('mu', [None, MU])
|
||||
def test_eigsolve_at_gamma_returns_a_static_mode_with_small_residual(mu) -> None:
|
||||
with numpy.errstate(divide='raise', invalid='raise', over='raise'):
|
||||
values, vectors = bloch.eigsolve(1, numpy.zeros(3), G_MATRIX, EPSILON, mu, y0=Y0, max_iters=200, tolerance=1e-9)
|
||||
residual = bloch.maxwell_operator(numpy.zeros(3), G_MATRIX, EPSILON, mu)(vectors[0]).ravel() - values[0] * vectors[0]
|
||||
assert abs(values[0]) < 1e-12
|
||||
assert numpy.linalg.norm(residual) < 1e-9
|
||||
|
||||
|
||||
def test_anisotropic_magnetic_eigensolve_matches_dense_reference() -> None:
|
||||
operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
|
||||
matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(H_SIZE, dtype=complex)])
|
||||
reference = numpy.linalg.eigvals(matrix)
|
||||
reference = reference[numpy.argsort(numpy.abs(reference))[:3]]
|
||||
rng = numpy.random.default_rng(9)
|
||||
guess = rng.normal(size=(3, H_SIZE)) + 1j * rng.normal(size=(3, H_SIZE))
|
||||
values, vectors = bloch.eigsolve(3, K0_GENERAL, G_MATRIX, EPSILON, MU, y0=guess, max_iters=300, tolerance=1e-10)
|
||||
assert_allclose(values, reference, atol=1e-11, rtol=1e-9)
|
||||
for value, vector in zip(values, vectors, strict=True):
|
||||
assert numpy.linalg.norm(matrix @ vector - value * vector) < 1e-9
|
||||
|
||||
|
||||
def test_eigsolve_without_initial_guess_returns_finite_modes() -> None:
|
||||
eigvals, eigvecs = bloch.eigsolve(
|
||||
1,
|
||||
|
|
|
|||
|
|
@ -2,7 +2,6 @@ from typing import cast
|
|||
|
||||
import numpy
|
||||
import pytest
|
||||
from scipy.sparse.linalg import spsolve
|
||||
|
||||
from ..fdmath import vec
|
||||
from ..fdfd import eme, waveguide_2d, waveguide_3d, waveguide_cyl
|
||||
|
|
@ -548,50 +547,6 @@ def test_get_s_fdfd_from_2d_ports_identical_guides_are_nearly_transmissive() ->
|
|||
assert abs(abs(ss[0, 1]) - 1.0) < 0.03
|
||||
|
||||
|
||||
def test_multimode_fdfd_calibrates_the_full_incoming_matrix(monkeypatch) -> None:
|
||||
# Exactly degenerate polarizations have different source/monitor bases.
|
||||
# A uniform guide must still have zero interface reflection and conserve flux.
|
||||
original_solve = waveguide_2d.solve_mode
|
||||
seed = 0
|
||||
|
||||
def seeded_solve(*args, **kwargs):
|
||||
nonlocal seed
|
||||
seed += 1
|
||||
return original_solve(*args, **kwargs, rng=seed)
|
||||
|
||||
decompositions = []
|
||||
original_decomposition = waveguide_3d._port_decomposition
|
||||
|
||||
def capture_decomposition(*args, **kwargs):
|
||||
result = original_decomposition(*args, **kwargs)
|
||||
decompositions.append(result)
|
||||
return result
|
||||
|
||||
monkeypatch.setattr(waveguide_2d, 'solve_mode', seeded_solve)
|
||||
monkeypatch.setattr(waveguide_3d, '_port_decomposition', capture_decomposition)
|
||||
epsilon = numpy.ones((3, 2, 2))
|
||||
widths = [[numpy.full(2, 0.5), numpy.full(2, 0.5)] for _ in range(2)]
|
||||
ss, diagnostics = waveguide_3d.get_s_fdfd_from_2d_ports(
|
||||
epsilon, epsilon, omega=0.5, dxes_2d=widths, mode_numbers=[0, 1],
|
||||
pml_cells=4, pml_buffer_cells=1, source_monitor_cells=2, monitor_interface_cells=2,
|
||||
source_scales=[2, 0.3, 1j, 0.7j], matrix_solver=spsolve, return_diagnostics=True,
|
||||
)
|
||||
incoming, outgoing = [], []
|
||||
for field in diagnostics['electric_fields']:
|
||||
li, lo = waveguide_3d._decompose_port_field(decompositions[0], field)
|
||||
ri, ro = waveguide_3d._decompose_port_field(decompositions[1], field)
|
||||
incoming.append(numpy.concatenate((li, ri)))
|
||||
outgoing.append(numpy.concatenate((lo, ro)))
|
||||
incoming, outgoing = numpy.column_stack(incoming), numpy.column_stack(outgoing)
|
||||
assert numpy.linalg.norm(incoming - numpy.diag(numpy.diag(incoming))) > 0.1
|
||||
assert_close(ss @ incoming, outgoing, atol=1e-11, rtol=1e-11)
|
||||
assert_close(diagnostics['incoming_amplitudes'], incoming, atol=1e-12)
|
||||
assert_close(ss[:2, :2], 0, atol=1e-11)
|
||||
assert_close(ss[2:, 2:], 0, atol=1e-11)
|
||||
assert_close(ss.conj().T @ diagnostics['outgoing_power_gram'] @ ss, diagnostics['incoming_power_gram'], atol=1e-11)
|
||||
assert_close(diagnostics['guided_powers'], 1, atol=1e-11)
|
||||
|
||||
|
||||
def test_get_s_fdfd_from_2d_ports_mismatch_returns_finite_diagnostics() -> None:
|
||||
left = _fdfd_cross_section(2.25)
|
||||
right = _fdfd_cross_section(1.5)
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@ import dataclasses
|
|||
import pytest # type: ignore
|
||||
import numpy
|
||||
from numpy.typing import NDArray
|
||||
from scipy.sparse.linalg import spsolve
|
||||
#from numpy.testing import assert_allclose, assert_array_equal
|
||||
|
||||
from .. import fdfd
|
||||
|
|
@ -146,9 +145,7 @@ def sim(
|
|||
omega = omega,
|
||||
dxes = dxes,
|
||||
epsilon = eps_vec,
|
||||
# These small, nearly static systems are ill-conditioned. Test the
|
||||
# Maxwell residual and flux identities independently of QMR convergence.
|
||||
matrix_solver = spsolve,
|
||||
matrix_solver_opts = dict(atol=1e-15, rtol=1e-11),
|
||||
)
|
||||
e = unvec(e_vec, shape[1:])
|
||||
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
|
|||
curl_term = curl_term / mu
|
||||
curl_term = pm * curl_term
|
||||
curl_term = ch(curl_term)
|
||||
return pe * (curl_term - OMEGA**2 * EPSILON * masked_e) + numpy.where(PEC, E_FIELD, 0)
|
||||
return pe * (curl_term - OMEGA**2 * EPSILON * masked_e)
|
||||
|
||||
|
||||
def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
|
||||
|
|
@ -48,7 +48,7 @@ def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
|
|||
curl_term = ch(masked_h)
|
||||
curl_term = pe * (curl_term / EPSILON)
|
||||
curl_term = ce(curl_term)
|
||||
return pm * (curl_term - OMEGA**2 * magnetic * masked_h) + numpy.where(PMC, H_FIELD, 0)
|
||||
return pm * (curl_term - OMEGA**2 * magnetic * masked_h)
|
||||
|
||||
|
||||
def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
|
||||
|
|
@ -95,8 +95,6 @@ def test_eh_full_matches_manual_block_operator_with_masks() -> None:
|
|||
|
||||
dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD))
|
||||
dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD))
|
||||
dense_e += numpy.where(PEC, E_FIELD, 0)
|
||||
dense_h += numpy.where(PMC, H_FIELD, 0)
|
||||
|
||||
assert_fields_close(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
|
||||
assert_fields_close(matrix_h, dense_h, atol=1e-10, rtol=1e-10)
|
||||
|
|
|
|||
|
|
@ -98,76 +98,3 @@ def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
|
|||
assert restored['dy'] == pytest.approx(0.3)
|
||||
assert numpy.isfinite(restored['E'][0]).all()
|
||||
assert numpy.isfinite(restored['H'][0]).all()
|
||||
|
||||
|
||||
@pytest.mark.parametrize('shape', [(4, 4), (4, 5), (5, 4), (5, 7)])
|
||||
def test_farfield_padding_preserves_aperture_origin(shape) -> None:
|
||||
rng = numpy.random.default_rng(7)
|
||||
fields = [rng.normal(size=shape) + 1j * rng.normal(size=shape) for _ in range(4)]
|
||||
target = (8, 10)
|
||||
padding = [(new // 2 - old // 2, new - old - (new // 2 - old // 2)) for old, new in zip(shape, target, strict=True)]
|
||||
padded = [numpy.pad(field, padding) for field in fields]
|
||||
automatic = farfield.near_to_farfield(fields[:2], fields[2:], 0.5, 0.7, padded_size=target)
|
||||
manual = farfield.near_to_farfield(padded[:2], padded[2:], 0.5, 0.7, padded_size=target)
|
||||
for key in ('E', 'H'):
|
||||
numpy.testing.assert_allclose(automatic[key], manual[key], atol=1e-13)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('size', [4, 8])
|
||||
def test_farfield_amplitude_is_independent_of_aperture_sampling(size) -> None:
|
||||
# A constant 2 by 2 wavelength aperture has integral 4 on either mesh.
|
||||
ones = numpy.ones((size, size), dtype=complex)
|
||||
zeros = numpy.zeros_like(ones)
|
||||
result = farfield.near_to_farfield([ones, zeros], [zeros, ones], 2 / size, 2 / size, padded_size=size)
|
||||
assert result['E'][0][size // 2, size // 2] == pytest.approx(8)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('shape', [(4, 4), (16, 12), (7, 9)])
|
||||
def test_inverse_resizing_matches_plane_wave_quadrature(shape) -> None:
|
||||
spectrum_shape = (8, 6)
|
||||
dkx, dky = 0.4, 0.3
|
||||
amplitude = 1.7 + 0.2j
|
||||
e = [numpy.zeros(spectrum_shape, dtype=complex) for _ in range(2)]
|
||||
e[0][spectrum_shape[0] // 2 + 1, spectrum_shape[1] // 2] = 2 * amplitude
|
||||
h = [-e[1], e[0]]
|
||||
result = farfield.far_to_nearfield(e, h, dkx, dky, padded_size=shape)
|
||||
dx, dy = 2 * numpy.pi / (shape[0] * dkx), 2 * numpy.pi / (shape[1] * dky)
|
||||
assert result['dx'] == pytest.approx(dx)
|
||||
assert result['dy'] == pytest.approx(dy)
|
||||
x = (numpy.arange(shape[0]) - shape[0] // 2) * dx
|
||||
expected = numpy.broadcast_to(amplitude * dkx * dky / (2 * numpy.pi) ** 2 * numpy.exp(-1j * dkx * x[:, None]), shape)
|
||||
numpy.testing.assert_allclose(result['E'][0], expected, atol=1e-14)
|
||||
numpy.testing.assert_allclose(result['H'][1], expected / numpy.sqrt(1 - (dkx / (2 * numpy.pi)) ** 2), atol=1e-14)
|
||||
|
||||
|
||||
def test_inverse_accepts_readonly_spectra_without_modifying_them() -> None:
|
||||
rng = numpy.random.default_rng(9)
|
||||
fields = [rng.normal(size=(8, 8)).astype(complex) for _ in range(4)]
|
||||
originals = [field.copy() for field in fields]
|
||||
for field in fields:
|
||||
field.flags.writeable = False
|
||||
result = farfield.far_to_nearfield(fields[:2], fields[2:], dkx=3, dky=3)
|
||||
assert numpy.isfinite(result['E']).all()
|
||||
for field, original in zip(fields, originals, strict=True):
|
||||
numpy.testing.assert_array_equal(field, original)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('direction', [-1, 1])
|
||||
def test_farfield_direction_and_roundtrip_match_outgoing_phasors(direction) -> None:
|
||||
shape = (16, 8)
|
||||
dx = dy = 0.25
|
||||
k = 2 * numpy.pi
|
||||
kx = direction * 2 * numpy.pi / (shape[0] * dx)
|
||||
kz = numpy.sqrt(k * k - kx * kx)
|
||||
x = (numpy.arange(shape[0]) - shape[0] // 2) * dx
|
||||
phase = numpy.broadcast_to(numpy.exp(-1j * kx * x[:, None]), shape)
|
||||
zeros = numpy.zeros_like(phase)
|
||||
electric, magnetic = [kz / k * phase, zeros], [zeros, phase]
|
||||
result = farfield.near_to_farfield(electric, magnetic, dx, dy, padded_size=shape)
|
||||
intensity = numpy.abs(result['E'][0]) ** 2 + numpy.abs(result['E'][1]) ** 2
|
||||
peak = numpy.unravel_index(numpy.argmax(intensity), shape)
|
||||
assert result['kx'][peak] == pytest.approx(kx)
|
||||
assert result['phi'][peak] == pytest.approx(numpy.arccos(kz / k))
|
||||
restored = farfield.far_to_nearfield(result['E'], result['H'], result['dkx'], result['dky'], padded_size=shape)
|
||||
numpy.testing.assert_allclose(restored['E'], electric, atol=1e-13)
|
||||
numpy.testing.assert_allclose(restored['H'], magnetic, atol=1e-13)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,4 @@
|
|||
import numpy
|
||||
import pytest
|
||||
|
||||
from ..fdmath import unvec, vec
|
||||
from ..fdfd import functional, operators
|
||||
|
|
@ -121,17 +120,3 @@ def test_poynting_e_cross_h_matches_sparse_operator() -> None:
|
|||
functional_result = functional.poynting_e_cross_h(DXES)(E_FIELD, H_FIELD)
|
||||
|
||||
assert_fields_match(functional_result, matrix_result)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('wide_input', ['electric', 'magnetic', 'widths'])
|
||||
def test_poynting_uses_promoted_precision_before_products(wide_input) -> None:
|
||||
shape = (2, 2, 2)
|
||||
e = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'electric' else numpy.complex64)
|
||||
h = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'magnetic' else numpy.complex64)
|
||||
e[1], h[2] = 1e20, 1e20j
|
||||
widths = [[numpy.ones(n, dtype=numpy.float64 if wide_input == 'widths' else numpy.float32) for n in shape] for _ in range(2)]
|
||||
with numpy.errstate(over='raise', invalid='raise'):
|
||||
result = functional.poynting_e_cross_h(widths)(e, h)
|
||||
expected = unvec(operators.poynting_e_cross(vec(e.astype(complex)), widths) @ vec(h.astype(complex)), shape)
|
||||
assert result.dtype == numpy.complex128
|
||||
assert_fields_close(result, expected, atol=0, rtol=1e-14)
|
||||
|
|
|
|||
|
|
@ -1,10 +1,8 @@
|
|||
from typing import cast
|
||||
|
||||
import numpy
|
||||
import pytest
|
||||
from scipy.sparse.linalg import spsolve
|
||||
|
||||
from ..fdfd import solvers, operators
|
||||
from ..fdfd import solvers
|
||||
from ._solver_cases import solver_plumbing_case
|
||||
from .utils import assert_close
|
||||
|
||||
|
|
@ -128,41 +126,3 @@ def test_generic_without_guess_does_not_inject_x0(monkeypatch) -> None:
|
|||
|
||||
assert 'x0' not in cast(dict[str, object], captured['kwargs'])
|
||||
assert_close(result, case.pr @ numpy.array([1.0, -1.0]))
|
||||
|
||||
|
||||
def test_default_solver_reports_breakdown_in_nonsingular_maxwell_system() -> None:
|
||||
dxes = [[numpy.ones(2), numpy.ones(1), numpy.ones(1)] for _ in range(2)]
|
||||
epsilon = numpy.full(6, 2.0)
|
||||
current = numpy.zeros(6, dtype=complex)
|
||||
current[2] = 1
|
||||
with pytest.raises(RuntimeError, match='QMR failed to converge.*info=-14'):
|
||||
solvers.generic(1.0, dxes, current, epsilon)
|
||||
direct = solvers.generic(1.0, dxes, current, epsilon, matrix_solver=spsolve)
|
||||
assert_close(operators.e_full(1.0, dxes, epsilon) @ direct, -1j * current)
|
||||
|
||||
|
||||
def test_qmr_iteration_limit_is_not_silently_accepted(monkeypatch) -> None:
|
||||
monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', lambda *args, **kwargs: (numpy.zeros(2), 1))
|
||||
with pytest.raises(RuntimeError, match='info=1'):
|
||||
solvers._scipy_qmr(solver_plumbing_case().a0, numpy.ones(2))
|
||||
|
||||
|
||||
@pytest.mark.parametrize('adjoint', [False, True])
|
||||
@pytest.mark.parametrize('direct', [False, True])
|
||||
def test_pec_constraints_remove_sources_and_initial_fields(adjoint, direct) -> None:
|
||||
dxes = [[numpy.ones(2), numpy.ones(1), numpy.ones(1)] for _ in range(2)]
|
||||
epsilon = numpy.full(6, 2.0)
|
||||
current = numpy.array([7, 0, 1, 0, 0, 0], dtype=complex)
|
||||
pec = numpy.array([1, 0, 0, 0, 0, 0])
|
||||
guess = numpy.ones(6, dtype=complex)
|
||||
options = {'rtol': 1e-10}
|
||||
kwargs = {'matrix_solver': spsolve} if direct else {'E_guess': guess, 'matrix_solver_opts': options}
|
||||
field = solvers.generic(0.7, dxes, current, epsilon, pec=pec, adjoint=adjoint, **kwargs)
|
||||
assert field[0] == 0
|
||||
assert numpy.isfinite(field).all()
|
||||
matrix = operators.e_full(0.7, dxes, epsilon, pec=pec)
|
||||
if adjoint:
|
||||
matrix = matrix.conj().T
|
||||
assert_close(matrix @ field, -0.7j * numpy.where(pec, 0, current), atol=1e-10)
|
||||
assert 'x0' not in options
|
||||
assert guess[0] == 1
|
||||
|
|
|
|||
|
|
@ -87,24 +87,6 @@ def test_curls_preserve_small_residual_between_large_derivative_terms(name) -> N
|
|||
assert_fields_close(result, expected, atol=1e-7, rtol=1e-12)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('name', ['curl_forward', 'curl_back'])
|
||||
@pytest.mark.parametrize('complex_field', [False, True])
|
||||
def test_curls_use_common_precision_for_mixed_width_dtypes(name, complex_field) -> None:
|
||||
shape = (1, 2, 2)
|
||||
field = numpy.zeros((3, *shape), dtype=numpy.complex64 if complex_field else numpy.float32)
|
||||
field[2, :, 1, :] = 1e8
|
||||
field[1, :, :, 1] = 1e8
|
||||
widths = [numpy.ones(1), numpy.full(2, 0.7, dtype=numpy.float32), numpy.full(2, float(numpy.float32(0.7)))]
|
||||
wide_widths = [dx.astype(numpy.float64) for dx in widths]
|
||||
expected = unvec(getattr(fd_operators, name)(wide_widths) @ vec(field), shape)
|
||||
assert expected[0, 0, 0, 0] == 0
|
||||
result = getattr(fd_functional, name)(widths)(field)
|
||||
assert_fields_close(result, expected, atol=1e-8, rtol=1e-12)
|
||||
parts = getattr(fd_functional, name + '_parts')(widths)(field)
|
||||
for component, (left, right) in enumerate(parts):
|
||||
assert_close(left + right, expected[component], atol=1e-8, rtol=1e-12)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('name', ['deriv_forward', 'deriv_back'])
|
||||
@pytest.mark.parametrize('axis', [0, 1, 2])
|
||||
@pytest.mark.parametrize('spaced', [False, True])
|
||||
|
|
|
|||
|
|
@ -91,22 +91,6 @@ def test_vec_cross_matches_pointwise_cross_product() -> None:
|
|||
assert_close(matrix_result, expected)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('name', ['deriv_forward', 'deriv_back'])
|
||||
@pytest.mark.parametrize('complex_widths', [False, True])
|
||||
def test_sparse_derivatives_widen_reciprocals_before_overflow(name, complex_widths) -> None:
|
||||
dtype = numpy.complex64 if complex_widths else numpy.float32
|
||||
width = numpy.array(1e-39 * (1 + 1j) if complex_widths else 1e-39, dtype=dtype)
|
||||
widths = [numpy.full(2, width, dtype=dtype), numpy.ones(1), numpy.ones(1)]
|
||||
wide_widths = [dx.astype(numpy.complex128 if complex_widths else numpy.float64) for dx in widths]
|
||||
field = numpy.array([0, 1e-39])
|
||||
with numpy.errstate(over='raise', invalid='raise'):
|
||||
matrix = getattr(operators, name)(widths)[0]
|
||||
result = matrix @ field
|
||||
assert numpy.isfinite(matrix.data).all()
|
||||
expected = getattr(operators, name)(wide_widths)[0] @ field
|
||||
assert_close(result, expected, atol=1e-12, rtol=1e-12)
|
||||
|
||||
|
||||
def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
|
||||
matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
|
||||
expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))
|
||||
|
|
|
|||
|
|
@ -1261,7 +1261,6 @@ def test_richer_complex_basis_reduces_driven_fdfd_monitor_residual() -> None:
|
|||
dxes_2d=dxes_2d,
|
||||
mode_numbers=[0],
|
||||
axis=2,
|
||||
prop_dx=0.15, # Resolve both port modes below the longitudinal band edge.
|
||||
pml_cells=4,
|
||||
pml_buffer_cells=2,
|
||||
source_monitor_cells=4,
|
||||
|
|
|
|||
|
|
@ -1,11 +1,9 @@
|
|||
import numpy
|
||||
import pytest
|
||||
from numpy.linalg import norm
|
||||
from numpy.testing import assert_allclose
|
||||
from scipy import linalg
|
||||
|
||||
from ..fdmath import vec
|
||||
from ..fdfd import waveguide_2d, scpml
|
||||
from ..fdfd import waveguide_2d
|
||||
|
||||
|
||||
OMEGA = 1 / 1500
|
||||
|
|
@ -122,38 +120,6 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
|
|||
assert 0.4 < ratio < 1.8
|
||||
|
||||
|
||||
@pytest.mark.parametrize('absorption', ['none', 'pml', 'material'])
|
||||
def test_sensitivity_matches_central_differences_on_graded_yee_mesh(absorption) -> None:
|
||||
rng = numpy.random.default_rng(4)
|
||||
shape = (4, 3)
|
||||
omega = 0.7
|
||||
epsilon = vec(1 + rng.uniform(0, 3, (3, *shape)))
|
||||
magnetic_widths = [rng.uniform(0.2, 1.2, n) for n in shape]
|
||||
electric_widths = [0.5 * (dx + numpy.roll(dx, -1)) for dx in magnetic_widths]
|
||||
dxes = [electric_widths, magnetic_widths]
|
||||
if absorption == 'pml':
|
||||
dxes = scpml.stretch_with_scpml(dxes, axis=0, polarity=-1, omega=omega, thickness=1, max_imag_ratio=0.3)
|
||||
elif absorption == 'material':
|
||||
epsilon = epsilon * (1 - 0.15j * rng.uniform(0.5, 1.5, epsilon.size))
|
||||
values, vectors = linalg.eig(waveguide_2d.operator_e(omega, dxes, epsilon).toarray())
|
||||
selected = numpy.argmax(values.real)
|
||||
beta = numpy.sqrt(values[selected])
|
||||
e, h = waveguide_2d.normalized_fields_e(vectors[:, selected], beta, omega, dxes, epsilon)
|
||||
sensitivity = waveguide_2d.sensitivity(e, h, beta, omega, dxes, epsilon)
|
||||
delta = 1e-5
|
||||
finite_difference = numpy.empty(epsilon.size, dtype=complex)
|
||||
for index in range(epsilon.size):
|
||||
roots = []
|
||||
for direction in (-1, 1):
|
||||
perturbed = epsilon.copy()
|
||||
perturbed[index] += direction * delta
|
||||
roots_squared = linalg.eigvals(waveguide_2d.operator_e(omega, dxes, perturbed).toarray())
|
||||
candidates = numpy.sqrt(roots_squared)
|
||||
roots.append(candidates[numpy.argmin(numpy.abs(candidates - beta))])
|
||||
finite_difference[index] = (roots[1] - roots[0]) / (2 * delta)
|
||||
assert_allclose(sensitivity, finite_difference, atol=1e-8, rtol=1e-5)
|
||||
|
||||
|
||||
def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu() -> None:
|
||||
epsilon = build_asymmetric_epsilon()
|
||||
mu = build_mu_profile()
|
||||
|
|
|
|||
|
|
@ -36,7 +36,7 @@ SCATTERING_SOURCE_SLICES = (slice(4, 5), slice(None), slice(None))
|
|||
SCATTERING_REFLECT_SLICES = (slice(10, 11), slice(None), slice(None))
|
||||
SCATTERING_TRANSMIT_SLICES = (slice(29, 30), slice(None), slice(None))
|
||||
SCATTERING_STEP_X = 18
|
||||
SCATTERING_WARMUP_PERIODS = 30
|
||||
SCATTERING_WARMUP_PERIODS = 10
|
||||
SCATTERING_ACCUMULATION_PERIODS = 10
|
||||
|
||||
|
||||
|
|
@ -221,21 +221,6 @@ def _build_cpml_params() -> list[list[dict[str, numpy.ndarray | float]]]:
|
|||
]
|
||||
|
||||
|
||||
def _build_cpml_frequency_dxes(base_dxes: list[list[numpy.ndarray]]) -> list[list[numpy.ndarray]]:
|
||||
"""Match the discrete CPML recurrence, including its Yee sample offsets."""
|
||||
widths = [[dx.astype(complex) for dx in family] for family in base_dxes]
|
||||
for axis, faces in enumerate(_build_cpml_params()):
|
||||
for face in faces:
|
||||
# psi[n] = p0 psi[n-1] + p1 p2 D[n]. Its phasor response adds
|
||||
# p1 / (1 - p0 exp(-i omega dt)) to the scaled derivative.
|
||||
# Forward electric derivatives use the H-update CPML parameters.
|
||||
for family, parameter in enumerate(('param_h', 'param_e')):
|
||||
p0, p1, p2 = face[parameter]
|
||||
response = p2 * (1 + p1 / (1 - p0 * numpy.exp(-1j * OMEGA * DT)))
|
||||
widths[family][axis][face['region'][axis]] /= response.ravel()
|
||||
return widths
|
||||
|
||||
|
||||
def _build_complex_pulse_waveform(total_steps: int) -> tuple[numpy.ndarray, complex]:
|
||||
source_phasor, _delay = gaussian_packet(wl=WAVELENGTH, dwl=PULSE_DWL, dt=DT, turn_on=1e-5)
|
||||
aa, cc, ss = source_phasor(numpy.arange(total_steps) + 0.5)
|
||||
|
|
@ -465,10 +450,7 @@ def _run_straight_waveguide_case(variant: str) -> WaveguideCalibrationResult:
|
|||
def _run_width_step_scattering_case() -> WaveguideScatteringResult:
|
||||
epsilon = _build_scattering_epsilon()
|
||||
base_dxes = _build_uniform_dxes(SCATTERING_SHAPE)
|
||||
# A correctly signed backward mode resolves the weak reflected wave;
|
||||
# incident leakage no longer hides PML differences or startup transients.
|
||||
stretched_dxes = _build_cpml_frequency_dxes(base_dxes)
|
||||
discrete_frequency = 2 * numpy.sin(OMEGA * DT / 2) / DT
|
||||
stretched_dxes = _build_stretched_dxes(base_dxes)
|
||||
|
||||
source_mode = waveguide_3d.solve_mode(
|
||||
0,
|
||||
|
|
@ -557,14 +539,14 @@ def _run_width_step_scattering_case() -> WaveguideScatteringResult:
|
|||
e_fdfd = unvec(
|
||||
fdfd.solvers.generic(
|
||||
J=vec(j_ph),
|
||||
omega=discrete_frequency,
|
||||
omega=OMEGA,
|
||||
dxes=stretched_dxes,
|
||||
epsilon=vec(epsilon),
|
||||
matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-7},
|
||||
),
|
||||
SCATTERING_SHAPE[1:],
|
||||
)
|
||||
h_fdfd = functional.e2h(discrete_frequency, stretched_dxes)(e_fdfd)
|
||||
h_fdfd = functional.e2h(OMEGA, stretched_dxes)(e_fdfd)
|
||||
|
||||
reflected_td = complex(vec(e_ph) @ vec(reflected_overlap).conj())
|
||||
reflected_fd = complex(vec(e_fdfd) @ vec(reflected_overlap).conj())
|
||||
|
|
|
|||
|
|
@ -2,49 +2,17 @@ import contextlib
|
|||
import io
|
||||
import numpy
|
||||
from numpy.linalg import norm
|
||||
from scipy import linalg
|
||||
import pytest
|
||||
import warnings
|
||||
|
||||
from ..fdmath import vec, unvec
|
||||
from ..fdmath.operators import deriv_forward, deriv_back
|
||||
from ..fdmath.functional import curl_forward, curl_back
|
||||
from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl
|
||||
|
||||
|
||||
OMEGA = 1 / 1500
|
||||
|
||||
|
||||
@pytest.mark.parametrize('axis', [0, 1, 2])
|
||||
@pytest.mark.parametrize('polarity', [-1, 1])
|
||||
def test_embedded_hybrid_mode_obeys_both_maxwell_curls(axis, polarity) -> None:
|
||||
omega = 0.7
|
||||
# Cyclic rotation preserves the orientation of the vector cross product.
|
||||
permutation = numpy.roll(numpy.arange(3), axis + 1)
|
||||
rng = numpy.random.default_rng(4)
|
||||
canonical = numpy.repeat((1 + rng.uniform(0, 3, (3, 4, 3)))[..., None], 9, axis=3)
|
||||
epsilon = canonical[permutation].transpose(0, *(permutation + 1))
|
||||
shape = epsilon.shape[1:]
|
||||
dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
|
||||
cross_dxes = [[numpy.full(4, 0.4), numpy.full(3, 0.4)] for _ in range(2)]
|
||||
values, vectors = linalg.eig(waveguide_2d.operator_e(omega, cross_dxes, vec(canonical[:, :, :, 0])).toarray())
|
||||
selected = numpy.argmax(values.real)
|
||||
slices = [slice(None)] * 3
|
||||
slices[axis] = slice(4, 5)
|
||||
mode = waveguide_3d.mode_from_2d(
|
||||
vectors[:, selected], numpy.sqrt(values[selected]), omega=omega, dxes=dxes,
|
||||
axis=axis, polarity=polarity, slices=slices, epsilon=epsilon,
|
||||
)
|
||||
e, h = [waveguide_3d.expand_e(mode[field], mode['wavenumber'], dxes, axis, polarity, slices) for field in ('E', 'H')]
|
||||
assert norm(e[axis]) > 0.01 * norm(e)
|
||||
assert norm(h[axis]) > 0.01 * norm(h)
|
||||
region = [slice(None)] * 4
|
||||
region[axis + 1] = slice(1, -1)
|
||||
ce, ch = curl_forward(dxes[0])(e), curl_back(dxes[1])(h)
|
||||
assert relative_residual(ce[tuple(region)], (-1j * omega * h)[tuple(region)]) < 1e-11
|
||||
assert relative_residual(ch[tuple(region)], (1j * omega * epsilon * e)[tuple(region)]) < 1e-11
|
||||
|
||||
|
||||
def build_waveguide_3d_mode(
|
||||
*,
|
||||
slice_start: int,
|
||||
|
|
@ -153,8 +121,8 @@ def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent(
|
|||
line_selection[axis] = slice(None)
|
||||
values = expanded[component][tuple(line_selection)]
|
||||
ratios = values[1:] / values[:-1]
|
||||
intervals = numpy.real(dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1])
|
||||
expected_ratios = numpy.exp(-1j * result['wavenumber'] * intervals)
|
||||
phase_positions = numpy.real(dxes[0][axis]).cumsum()
|
||||
expected_ratios = numpy.exp(-1j * result['wavenumber'] * numpy.diff(phase_positions))
|
||||
|
||||
numpy.testing.assert_allclose(ratios, expected_ratios, rtol=1e-6, atol=1e-9)
|
||||
|
||||
|
|
@ -422,7 +390,7 @@ def test_waveguide_cyl_reconstruction_is_finite_radius_consistent(use_mu: bool)
|
|||
eps_r, eps_y, eps_z = numpy.split(epsilon, 3)
|
||||
mu_r, mu_y, mu_z = numpy.split(mu_full, 3)
|
||||
|
||||
ez_rhs = (
|
||||
ez_rhs = ta * (
|
||||
(Dbx @ (tb * eps_r * er))
|
||||
+ (Dby @ (ta * eps_y * ey))
|
||||
) / eps_z
|
||||
|
|
@ -469,79 +437,8 @@ def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None:
|
|||
rmin=rmin,
|
||||
)
|
||||
|
||||
numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / numpy.average(rb, weights=dxes[0][0]))
|
||||
numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / numpy.average(ra, weights=dxes[1][0]))
|
||||
|
||||
|
||||
@pytest.mark.parametrize('use_mu', [False, True])
|
||||
def test_cylindrical_modes_obey_the_full_ampere_curl(use_mu) -> None:
|
||||
shape, omega, rmin = (4, 3), 0.7, 2.0
|
||||
dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
|
||||
rng = numpy.random.default_rng(4)
|
||||
epsilon = vec(1 + rng.uniform(0, 3, (3, *shape)))
|
||||
mu = vec(1 + rng.uniform(0, 2, (3, *shape))) if use_mu else numpy.ones_like(epsilon)
|
||||
matrix = waveguide_cyl.cylindrical_operator(omega, dxes, epsilon, rmin, mu).toarray()
|
||||
values, vectors = linalg.eig(matrix)
|
||||
selected = numpy.argmax(values.real)
|
||||
beta = numpy.sqrt(values[selected])
|
||||
electric = waveguide_cyl.exy2e(beta * rmin, dxes, rmin, epsilon) @ vectors[:, selected]
|
||||
magnetic = waveguide_cyl.e2h(beta * rmin, omega, dxes, rmin, mu) @ electric
|
||||
ta, tb = (metric.diagonal() for metric in waveguide_cyl.dxes2T(dxes, rmin))
|
||||
dbx, dby = deriv_back(dxes[1])
|
||||
hr, hy, hz = numpy.split(magnetic, 3)
|
||||
# Physical cylindrical Ampere curl, independently of electric reconstruction.
|
||||
curl_h = numpy.concatenate((
|
||||
1j * beta * hy / tb + dby @ hz,
|
||||
-1j * beta * hr / ta - (dbx @ (tb * hz)) / ta,
|
||||
-dby @ hr + dbx @ hy,
|
||||
))
|
||||
numpy.testing.assert_allclose(curl_h, 1j * omega * epsilon * electric, atol=1e-11, rtol=1e-10)
|
||||
|
||||
|
||||
def test_cylindrical_mean_radius_is_invariant_under_constant_density_refinement() -> None:
|
||||
for widths in ([1.0, 3.0], [1.0, 1.0, 1.0, 1.0]):
|
||||
n = len(widths)
|
||||
dxes = [[numpy.array(widths), numpy.ones(1)] for _ in range(2)]
|
||||
electric_xy = numpy.concatenate((numpy.ones(n), numpy.zeros(n)))[None, :].astype(complex)
|
||||
beta = waveguide_cyl.linear_wavenumbers(electric_xy, [12.0], numpy.ones(3 * n), dxes, rmin=10.0)
|
||||
numpy.testing.assert_allclose(beta, 1.0)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('axis', [0, 1, 2])
|
||||
@pytest.mark.parametrize('polarity', [-1, 1])
|
||||
def test_graded_expansion_uses_each_components_actual_yee_intervals(axis, polarity) -> None:
|
||||
shape = [1, 1, 1]
|
||||
shape[axis] = 4
|
||||
dxes = [[numpy.ones(n) for n in shape] for _ in range(2)]
|
||||
dxes[0][axis] = numpy.array([1.0, 2.0, 3.0, 4.0])
|
||||
dxes[1][axis] = numpy.array([0.6, 0.8, 1.0, 1.2])
|
||||
source = [slice(None)] * 3
|
||||
source[axis] = slice(1, 2)
|
||||
field = numpy.zeros((3, *shape), dtype=complex)
|
||||
field[(slice(None), *source)] = 1 + 0.5j
|
||||
expanded = waveguide_3d.expand_e(field, 0.2, dxes, axis, polarity, source)
|
||||
for component in range(3):
|
||||
line = expanded[component].ravel()
|
||||
intervals = dxes[1][axis][1:] if component == axis else dxes[0][axis][:-1]
|
||||
numpy.testing.assert_allclose(line[1:] / line[:-1], numpy.exp(-1j * polarity * 0.2 * intervals))
|
||||
numpy.testing.assert_array_equal(expanded[(slice(None), *source)], field[(slice(None), *source)])
|
||||
|
||||
|
||||
@pytest.mark.parametrize('beta', [2.0, 3.0])
|
||||
def test_grid_band_edge_rejects_power_normalization_and_selects_decay(beta) -> None:
|
||||
shape = (1, 1, 5)
|
||||
dxes = [[numpy.ones(n) for n in shape] for _ in range(2)]
|
||||
kwargs = dict(omega=1.0, dxes=dxes, axis=2, polarity=1,
|
||||
slices=(slice(None), slice(None), slice(0, 1)), epsilon=numpy.full((3, *shape), beta ** 2))
|
||||
with pytest.raises(ValueError, match='grid band edge'):
|
||||
waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, **kwargs)
|
||||
mode = waveguide_3d.mode_from_2d(numpy.array([1, 0], dtype=complex), beta, normalization='lorentz', **kwargs)
|
||||
assert mode['wavenumber'].imag <= 0
|
||||
expanded = waveguide_3d.expand_e(mode['E'], mode['wavenumber'], dxes, 2, 1, kwargs['slices'])
|
||||
ratio = abs(expanded[0, 0, 0, 1] / expanded[0, 0, 0, 0])
|
||||
assert ratio <= 1 + 1e-14
|
||||
if beta > 2:
|
||||
assert ratio < 1
|
||||
numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / rb.mean())
|
||||
numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / ra.mean())
|
||||
|
||||
|
||||
def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None:
|
||||
|
|
|
|||
Loading…
Reference in a new issue