jan/fdfd_tools
2
2
Fork 1
Code Issues Pull requests Releases Wiki Activity

rename to meanas and split fdtd/fdfd

Browse source
This commit is contained in:
Jan Petykiewicz 2019-08-04 13:48:41 -07:00
commit f61bcf3dfa
22 changed files with 519 additions and 439 deletions
Download patch file Download diff file Expand all files Collapse all files

48
meanas/__init__.py Normal file
View file

@ -0,0 +1,48 @@
"""
Electromagnetic simulation tools
This package is intended for building simulation inputs, analyzing
simulation outputs, and running short simulations on unspecialized hardware.
It is designed to provide tooling and a baseline for other, high-performance
purpose- and hardware-specific solvers.
**Contents**
- Finite difference frequency domain (FDFD)
* Library of sparse matrices for representing the electromagnetic wave
equation in 3D, as well as auxiliary matrices for conversion between fields
* Waveguide mode operators
* Waveguide mode eigensolver
* Stretched-coordinate PML boundaries (SCPML)
* Functional versions of most operators
* Anisotropic media (limited to diagonal elements eps_xx, eps_yy, eps_zz, mu_xx, ...)
* Arbitrary distributions of perfect electric and magnetic conductors (PEC / PMC)
- Finite difference time domain (FDTD)
* Basic Maxwell time-steps
* Poynting vector and energy calculation
* Convolutional PMLs
This package does *not* provide a fast matrix solver, though by default
```meanas.fdfd.solvers.generic(...)``` will call
```scipy.sparse.linalg.qmr(...)``` to perform a solve.
For 2D FDFD problems this should be fine; likewise, the waveguide mode
solver uses scipy's eigenvalue solver, with reasonable results.
For solving large (or 3D) FDFD problems, I recommend a GPU-based iterative
solver, such as [opencl_fdfd](https://mpxd.net/code/jan/opencl_fdfd) or
those included in [MAGMA](http://icl.cs.utk.edu/magma/index.html)). Your
solver will need the ability to solve complex symmetric (non-Hermitian)
linear systems, ideally with double precision.
Dependencies:
- numpy
- scipy
"""
from .types import dx_lists_t, field_t, vfield_t, field_updater
from .vectorization import vec, unvec
__author__ = 'Jan Petykiewicz'
version = '0.5'

120
meanas/eigensolvers.py Normal file
View file

@ -0,0 +1,120 @@
"""
Solvers for eigenvalue / eigenvector problems
"""
from typing import Tuple, List
import numpy
from numpy.linalg import norm
from scipy import sparse
import scipy.sparse.linalg as spalg
def power_iteration(operator: sparse.spmatrix,
guess_vector: numpy.ndarray = None,
iterations: int = 20,
) -> Tuple[complex, numpy.ndarray]:
"""
Use power iteration to estimate the dominant eigenvector of a matrix.
:param operator: Matrix to analyze.
:param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
:param iterations: Number of iterations to perform. Default 20.
:return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
"""
if numpy.any(numpy.equal(guess_vector, None)):
v = numpy.random.rand(operator.shape[0])
else:
v = guess_vector
for _ in range(iterations):
v = operator @ v
v /= norm(v)
lm_eigval = v.conj() @ (operator @ v)
return lm_eigval, v
def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperator,
guess_vector: numpy.ndarray,
iterations: int = 40,
tolerance: float = 1e-13,
solver=None,
) -> Tuple[complex, numpy.ndarray]:
"""
Use Rayleigh quotient iteration to refine an eigenvector guess.
:param operator: Matrix to analyze.
:param guess_vector: Eigenvector to refine.
:param iterations: Maximum number of iterations to perform. Default 40.
:param tolerance: Stop iteration if (A - I*eigenvalue) @ v < tolerance.
Default 1e-13.
:param solver: Solver function of the form x = solver(A, b).
By default, use scipy.sparse.spsolve for sparse matrices and
scipy.sparse.bicgstab for general LinearOperator instances.
:return: (eigenvalue, eigenvector)
"""
try:
_test = operator - sparse.eye(operator.shape[0])
shift = lambda eigval: eigval * sparse.eye(operator.shape[0])
if solver is None:
solver = spalg.spsolve
except TypeError:
shift = lambda eigval: spalg.LinearOperator(shape=operator.shape,
dtype=operator.dtype,
matvec=lambda v: eigval * v)
if solver is None:
solver = lambda A, b: spalg.bicgstab(A, b)[0]
v = guess_vector
v /= norm(v)
for _ in range(iterations):
eigval = v.conj() @ (operator @ v)
if norm(operator @ v - eigval * v) < tolerance:
break
shifted_operator = operator - shift(eigval)
v = solver(shifted_operator, v)
v /= norm(v)
return eigval, v
def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator,
how_many: int,
negative: bool = False,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
"""
Find the largest-magnitude positive-only (or negative-only) eigenvalues and
eigenvectors of the provided matrix.
:param operator: Matrix to analyze.
:param how_many: How many eigenvalues to find.
:param negative: Whether to find negative-only eigenvalues.
Default False (positive only).
:return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
eigenvectors[:, k] corresponds to the k-th eigenvalue
"""
# Use power iteration to estimate the dominant eigenvector
lm_eigval, _ = power_iteration(operator)
'''
Shift by the absolute value of the largest eigenvalue, then find a few of the
largest-magnitude (shifted) eigenvalues. A positive shift ensures that we find the
largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the
range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval)
'''
shift = numpy.abs(lm_eigval)
if negative:
shift *= -1
# Try to combine, use general LinearOperator if we fail
try:
shifted_operator = operator + shift * sparse.eye(operator.shape[0])
except TypeError:
shifted_operator = operator + spalg.LinearOperator(shape=operator.shape,
matvec=lambda v: shift * v)
shifted_eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50)
eigenvalues = shifted_eigenvalues - shift
k = eigenvalues.argsort()
return eigenvalues[k], eigenvectors[:, k]

649
meanas/fdfd/bloch.py Normal file
View file

@ -0,0 +1,649 @@
'''
Bloch eigenmode solver/operators
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
1/mu * curl(1/eps * curl(H)) = (w/c)^2 H
into
conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
where:
- 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
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 with
k @ h = kx hx + ky hy + kz hz = 0 = hk
h = hk + hm + hn = hm + hn
k = kk + km + kn = kk = |k|
k x h = (ky hz - kz hy,
kz hx - kx hz,
kx hy - ky hx)
= ((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
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) + (kk hm)(nx, ny, nz)
= |k| (hm * (nx, ny, nz) - hn * (mx, my, mz))
where h is shorthand for H_k, (...)_kmn deontes the (k, m, n) basis,
and e.g. hm is the component of h in the m direction.
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)
'''
from typing import Tuple, Callable
import logging
import numpy
from numpy import pi, real, trace
from numpy.fft import fftfreq
import scipy
import scipy.optimize
from scipy.linalg import norm
import scipy.sparse.linalg as spalg
from . import field_t
logger = logging.getLogger(__name__)
try:
import pyfftw.interfaces.numpy_fft
import pyfftw.interfaces
import multiprocessing
logger.info('Using pyfftw')
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(3600)
fftw_args = {
'threads': multiprocessing.cpu_count(),
'overwrite_input': True,
'planner_effort': 'FFTW_EXHAUSTIVE',
}
def fftn(*args, **kwargs):
return pyfftw.interfaces.numpy_fft.fftn(*args, **kwargs, **fftw_args)
def ifftn(*args, **kwargs):
return pyfftw.interfaces.numpy_fft.ifftn(*args, **kwargs, **fftw_args)
except ImportError:
from numpy.fft import fftn, ifftn
logger.info('Using numpy fft')
def generate_kmn(k0: numpy.ndarray,
G_matrix: numpy.ndarray,
shape: numpy.ndarray
) -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]:
"""
Generate a (k, m, n) orthogonal basis for each k-vector in the simulation grid.
:param k0: [k0x, k0y, k0z], Bloch wavevector, in G basis.
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param shape: [nx, ny, nz] shape of the simulation grid.
:return: (|k|, m, n) where |k| has shape tuple(shape) + (1,)
and m, n have shape tuple(shape) + (3,).
All are given in the xyz basis (e.g. |k|[0,0,0] = norm(G_matrix @ k0)).
"""
k0 = numpy.array(k0)
Gi_grids = numpy.meshgrid(*(fftfreq(n, 1/n) for n in shape[:3]), indexing='ij')
Gi = numpy.stack(Gi_grids, axis=3)
k_G = k0[None, None, None, :] - Gi
k_xyz = numpy.rollaxis(G_matrix @ numpy.rollaxis(k_G, 3, 2), 3, 2)
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[:, :, :, 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]
z_non0 = numpy.any(k_xyz != 0, axis=3)
if numpy.any(z_non0):
v = numpy.cross(k_xyz[z_non0], m[z_non0])
n[z_non0, :] = v / norm(v, axis=1)[:, None]
k_mag = norm(k_xyz, axis=3)[:, :, :, None]
return k_mag, m, n
def maxwell_operator(k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
mu: field_t = None
) -> Callable[[numpy.ndarray], numpy.ndarray]:
"""
Generate the Maxwell operator
conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
which is the spatial-frequency-space representation of
1/mu * curl(1/eps * curl(___))
The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
See the module-level docstring for more information.
:param k0: Bloch wavevector, [k0x, k0y, k0z].
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
:param mu: Magnetic permability distribution for the simulation.
Default None (1 everywhere).
:return: Function which applies the maxwell operator to h_mn.
"""
shape = epsilon[0].shape + (1,)
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
epsilon = numpy.stack(epsilon, 3)
if mu is not None:
mu = numpy.stack(mu, 3)
def operator(h: numpy.ndarray):
"""
Maxwell operator for Bloch eigenmode simulation.
h is complex 2-field in (m, n) basis, vectorized
:param h: Raveled h_mn; size (2 * epsilon[0].size).
:return: Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x 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
# cross product and transform into xyz basis
d_xyz = (n * hin_m -
m * hin_n) * k_mag
# divide by epsilon
e_xyz = fftn(ifftn(d_xyz, axes=range(3)) / epsilon, axes=range(3))
# cross product and transform into mn basis
b_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] * -k_mag
b_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] * +k_mag
if mu is None:
h_m, h_n = b_m, b_n
else:
# transform from mn to xyz
b_xyz = (m * b_m[:, :, :, None] +
n * b_n[:, :, :, None])
# divide by mu
h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, 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 numpy.hstack((h_m.ravel(), h_n.ravel()))
return operator
def hmn_2_exyz(k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
) -> Callable[[numpy.ndarray], field_t]:
"""
Generate an operator which converts a vectorized spatial-frequency-space
h_mn into an E-field distribution, i.e.
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)
See the module-level docstring for more information.
:param k0: Bloch wavevector, [k0x, k0y, k0z].
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
:return: Function for converting h_mn into E_xyz
"""
shape = epsilon[0].shape + (1,)
epsilon = numpy.stack(epsilon, 3)
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
def operator(h: numpy.ndarray) -> field_t:
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
d_xyz = (n * hin_m -
m * hin_n) * k_mag
# divide by epsilon
return [ei for ei in numpy.rollaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3)]
return operator
def hmn_2_hxyz(k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t
) -> Callable[[numpy.ndarray], field_t]:
"""
Generate an operator which converts a vectorized spatial-frequency-space
h_mn into an H-field distribution, i.e.
ifft(h_mn)
The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
See the module-level docstring for more information.
:param k0: Bloch wavevector, [k0x, k0y, k0z].
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
Only epsilon[0].shape is used.
:return: Function for converting h_mn into H_xyz
"""
shape = epsilon[0].shape + (1,)
_k_mag, m, n = generate_kmn(k0, G_matrix, shape)
def operator(h: numpy.ndarray):
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
h_xyz = (m * hin_m +
n * hin_n)
return [ifftn(hi) for hi in numpy.rollaxis(h_xyz, 3)]
return operator
def inverse_maxwell_operator_approx(k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
mu: field_t = None
) -> Callable[[numpy.ndarray], numpy.ndarray]:
"""
Generate an approximate inverse of the Maxwell operator,
ik x conv(eps_k, ik x conv(mu_k, ___))
which can be used to improve the speed of ARPACK in shift-invert mode.
See the module-level docstring for more information.
:param k0: Bloch wavevector, [k0x, k0y, k0z].
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
:param mu: Magnetic permability distribution for the simulation.
Default None (1 everywhere).
:return: Function which applies the approximate inverse of the maxwell operator to h_mn.
"""
shape = epsilon[0].shape + (1,)
epsilon = numpy.stack(epsilon, 3)
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
if mu is not None:
mu = numpy.stack(mu, 3)
def operator(h: numpy.ndarray):
"""
Approximate inverse Maxwell operator for Bloch eigenmode simulation.
h is complex 2-field in (m, n) basis, vectorized
:param h: Raveled h_mn; size (2 * epsilon[0].size).
:return: 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[:, :, :, None] +
n * hin_n[:, :, :, None])
# multiply by mu
b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3))
# transform back to mn
b_m = numpy.sum(b_xyz * m, axis=3)
b_n = numpy.sum(b_xyz * n, axis=3)
# cross product and transform into xyz basis
e_xyz = (n * b_m -
m * b_n) / k_mag
# multiply by epsilon
d_xyz = fftn(ifftn(e_xyz, axes=range(3)) * epsilon, axes=range(3))
# cross product and transform into mn basis crossinv_t2c
h_m = numpy.sum(d_xyz * n, axis=3)[:, :, :, None] / +k_mag
h_n = numpy.sum(d_xyz * m, axis=3)[:, :, :, None] / -k_mag
return numpy.hstack((h_m.ravel(), h_n.ravel()))
return operator
def find_k(frequency: float,
tolerance: float,
direction: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
mu: field_t = None,
band: int = 0,
k_min: float = 0,
k_max: float = 0.5,
solve_callback: Callable = None
) -> Tuple[numpy.ndarray, float]:
"""
Search for a bloch vector that has a given frequency.
:param frequency: Target frequency.
:param tolerance: Target frequency tolerance.
:param direction: k-vector direction to search along.
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
:param mu: Magnetic permability distribution for the simulation.
Default None (1 everywhere).
:param band: Which band to search in. Default 0 (lowest frequency).
return: (k, actual_frequency) The found k-vector and its frequency
"""
direction = numpy.array(direction) / norm(direction)
def get_f(k0_mag: float, band: int = 0):
k0 = direction * k0_mag
n, v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
f = numpy.sqrt(numpy.abs(numpy.real(n[band])))
if solve_callback:
solve_callback(k0_mag, n, v, f)
return f
res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency),
(k_min + k_max) / 2,
method='Bounded',
bounds=(k_min, k_max),
options={'xatol': abs(tolerance)})
return res.x * direction, res.fun + frequency
def eigsolve(num_modes: int,
k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
mu: field_t = None,
tolerance: float = 1e-20,
max_iters: int = 10000,
reset_iters: int = 100,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
"""
Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
k0 of the specified structure.
:param k0: Bloch wavevector, [k0x, k0y, k0z].
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
:param epsilon: Dielectric constant distribution for the simulation.
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
:param mu: Magnetic permability distribution for the simulation.
Default None (1 everywhere).
:param tolerance: Solver stops when fractional change in the objective
trace(Z.H @ A @ Z @ inv(Z Z.H)) is smaller than the tolerance
:return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
vector eigenvectors[i, :]
"""
h_size = 2 * epsilon[0].size
kmag = norm(G_matrix @ k0)
'''
Generate the operators
'''
mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
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)
y_shape = (h_size, num_modes)
prev_E = 0
d_scale = 1
prev_traceGtKG = 0
#prev_theta = 0.5
D = numpy.zeros(shape=y_shape, dtype=complex)
y0 = None
if y0 is None:
Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
else:
Z = y0
while True:
Z *= num_modes / norm(Z)
ZtZ = Z.conj().T @ Z
try:
U = numpy.linalg.inv(ZtZ)
except numpy.linalg.LinAlgError:
Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
continue
trace_U = real(trace(U))
if trace_U > 1e8 * num_modes:
Z = Z @ scipy.linalg.sqrtm(U).conj().T
prev_traceGtKG = 0
continue
break
for i in range(max_iters):
ZtZ = Z.conj().T @ Z
U = numpy.linalg.inv(ZtZ)
AZ = scipy_op @ Z
AZU = AZ @ U
ZtAZU = Z.conj().T @ AZU
E_signed = real(trace(ZtAZU))
sgn = numpy.sign(E_signed)
E = numpy.abs(E_signed)
G = (AZU - Z @ U @ ZtAZU) * sgn
if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7):
logger.info('Optimization succeded: {} - 5e-8 < {} * {} / 2'.format(abs(E - prev_E), tolerance, E + prev_E))
break
KG = scipy_iop @ G
traceGtKG = _rtrace_AtB(G, KG)
if prev_traceGtKG == 0 or i % reset_iters == 0:
logger.info('CG reset')
gamma = 0
else:
gamma = traceGtKG / prev_traceGtKG
D = gamma / d_scale * D + KG
d_scale = num_modes / norm(D)
D *= d_scale
ZtAZ = Z.conj().T @ AZ
AD = scipy_op @ D
DtD = D.conj().T @ D
DtAD = D.conj().T @ AD
symZtD = _symmetrize(Z.conj().T @ D)
symZtAD = _symmetrize(Z.conj().T @ AD)
Qi_memo = [None, None]
def Qi_func(theta):
nonlocal Qi_memo
if Qi_memo[0] == theta:
return Qi_memo[1]
c = numpy.cos(theta)
s = numpy.sin(theta)
Q = c*c * ZtZ + s*s * DtD + 2*s*c * symZtD
try:
Qi = numpy.linalg.inv(Q)
except numpy.linalg.LinAlgError:
logger.info('taylor Qi')
# if c or s small, taylor expand
if c < 1e-4 * s and c != 0:
DtDi = numpy.linalg.inv(DtD)
Qi = DtDi / (s*s) - 2*c/(s*s*s) * (DtDi @ (DtDi @ symZtD).conj().T)
elif s < 1e-4 * c and s != 0:
ZtZi = numpy.linalg.inv(ZtZ)
Qi = ZtZi / (c*c) - 2*s/(c*c*c) * (ZtZi @ (ZtZi @ symZtD).conj().T)
else:
raise Exception('Inexplicable singularity in trace_func')
Qi_memo[0] = theta
Qi_memo[1] = Qi
return Qi
def trace_func(theta):
c = numpy.cos(theta)
s = numpy.sin(theta)
Qi = Qi_func(theta)
R = c*c * ZtAZ + s*s * DtAD + 2*s*c * symZtAD
trace = _rtrace_AtB(R, Qi)
return numpy.abs(trace)
'''
def trace_deriv(theta):
Qi = Qi_func(theta)
c2 = numpy.cos(2 * theta)
s2 = numpy.sin(2 * theta)
F = -0.5*s2 * (ZtAZ - DtAD) + c2 * symZtAD
trace_deriv = _rtrace_AtB(Qi, F)
G = Qi @ F.conj().T @ Qi.conj().T
H = -0.5*s2 * (ZtZ - DtD) + c2 * symZtD
trace_deriv -= _rtrace_AtB(G, H)
trace_deriv *= 2
return trace_deriv * sgn
U_sZtD = U @ symZtD
dE = 2.0 * (_rtrace_AtB(U, symZtAD) -
_rtrace_AtB(ZtAZU, U_sZtD))
d2E = 2 * (_rtrace_AtB(U, DtAD) -
_rtrace_AtB(ZtAZU, U @ (DtD - 4 * symZtD @ U_sZtD)) -
4 * _rtrace_AtB(U, symZtAD @ U_sZtD))
# Newton-Raphson to find a root of the first derivative:
theta = -dE/d2E
if d2E < 0 or abs(theta) >= pi:
theta = -abs(prev_theta) * numpy.sign(dE)
# theta, new_E, new_dE = linmin(theta, E, dE, 0.1, min(tolerance, 1e-6), 1e-14, 0, -numpy.sign(dE) * K_PI, trace_func)
theta, n, _, new_E, _, _new_dE = scipy.optimize.line_search(trace_func, trace_deriv, xk=theta, pk=numpy.ones((1,1)), gfk=dE, old_fval=E, c1=min(tolerance, 1e-6), c2=0.1, amax=pi)
'''
result = scipy.optimize.minimize_scalar(trace_func, bounds=(0, pi), tol=tolerance)
new_E = result.fun
theta = result.x
improvement = numpy.abs(E - new_E) * 2 / numpy.abs(E + new_E)
logger.info('linmin improvement {}'.format(improvement))
Z *= numpy.cos(theta)
Z += D * numpy.sin(theta)
prev_traceGtKG = traceGtKG
#prev_theta = theta
prev_E = E
'''
Recover eigenvectors from Z
'''
U = numpy.linalg.inv(Z.conj().T @ Z)
Y = Z @ scipy.linalg.sqrtm(U)
W = Y.conj().T @ (scipy_op @ Y)
eigvals, W_eigvecs = numpy.linalg.eig(W)
eigvecs = Y @ W_eigvecs
for i in range(len(eigvals)):
v = eigvecs[:, i]
n = eigvals[i]
v /= norm(v)
eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v )
f = numpy.sqrt(-numpy.real(n))
df = numpy.sqrt(-numpy.real(n + eigness))
neff_err = kmag * (1/df - 1/f)
logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err))
order = numpy.argsort(numpy.abs(eigvals))
return eigvals[order], eigvecs.T[order]
'''
def linmin(x_guess, f0, df0, x_max, f_tol=0.1, df_tol=min(tolerance, 1e-6), x_tol=1e-14, x_min=0, linmin_func):
if df0 > 0:
x0, f0, df0 = linmin(-x_guess, f0, -df0, -x_max, f_tol, df_tol, x_tol, -x_min, lambda q, dq: -linmin_func(q, dq))
return -x0, f0, -df0
elif df0 == 0:
return 0, f0, df0
else:
x = x_guess
fx = f0
dfx = df0
isave = numpy.zeros((2,), numpy.intc)
dsave = numpy.zeros((13,), float)
x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
x_min, x_max, isave, dsave)
for i in range(int(1e6)):
if task != 'F':
logging.info('search converged in {} iterations'.format(i))
break
fx = f(x, dfx)
x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
x_min, x_max, isave, dsave)
return x, fx, dfx
'''
def _rtrace_AtB(A, B):
return real(numpy.sum(A.conj() * B))
def _symmetrize(A):
return (A + A.conj().T) * 0.5

220
meanas/fdfd/farfield.py Normal file
View file

@ -0,0 +1,220 @@
"""
Functions for performing near-to-farfield transformation (and the reverse).
"""
from typing import Dict, List
import numpy
from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
from numpy import pi
def near_to_farfield(E_near: List[numpy.ndarray],
H_near: List[numpy.ndarray],
dx: float,
dy: float,
padded_size: List[int] = None
) -> Dict[str]:
"""
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.
:param E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).
:param H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).
:param dx: Cell size along x-dimension, in units of wavelength.
:param dy: Cell size along y-dimension, in units of wavelength.
:param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
Powers of 2 are most efficient for FFT computation.
Default is the smallest power of 2 larger than the input, for each axis.
:returns: Dict with keys
'E_far': Normalized E-field farfield; multiply by
(i k exp(-i k r) / (4 pi r)) to get the actual field value.
'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).
'dkx', 'dky': step size for kx and ky, normalized to wavelength.
'theta': arctan2(ky, kx) corresponding to each (kx, ky).
This is the angle in the x-y plane, counterclockwise from above, starting from +x.
'phi': arccos(kz / k) corresponding to each (kx, ky).
This is the angle away from +z.
"""
if not len(E_near) == 2:
raise Exception('E_near must be a length-2 list of ndarrays')
if not len(H_near) == 2:
raise Exception('H_near must be a length-2 list of ndarrays')
s = E_near[0].shape
if not all(s == f.shape for f in E_near + H_near):
raise Exception('All fields must be the same shape!')
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)
En_fft = [fftshift(fft2(fftshift(Eni), s=padded_size)) for Eni in E_near]
Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_size)) for Hni in H_near]
# Propagation vectors kx, ky
k = 2 * pi
kxx = 2 * pi * fftshift(fftfreq(padded_size[0], dx))
kyy = 2 * pi * fftshift(fftfreq(padded_size[1], dy))
kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
kxy2 = kx * kx + ky * ky
kxy = numpy.sqrt(kxy2)
kz = numpy.sqrt(k * k - kxy2)
sin_th = ky / kxy
cos_th = kx / kxy
cos_phi = kz / k
sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
# Normalized vector potentials N, L
N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th,
Hn_fft[1] * sin_th + Hn_fft[0] * cos_th]
L = [ En_fft[1] * cos_phi * cos_th - En_fft[0] * cos_phi * sin_th,
-En_fft[1] * sin_th - En_fft[0] * cos_th]
E_far = [-L[1] - N[0],
L[0] - N[1]]
H_far = [-E_far[1],
E_far[0]]
theta = numpy.arctan2(ky, kx)
phi = numpy.arccos(cos_phi)
theta[numpy.logical_and(kx == 0, ky == 0)] = 0
phi[numpy.logical_and(kx == 0, ky == 0)] = 0
# Zero fields beyond valid (phi, theta)
invalid_ind = kxy2 >= k * k
theta[invalid_ind] = 0
phi[invalid_ind] = 0
for i in range(2):
E_far[i][invalid_ind] = 0
H_far[i][invalid_ind] = 0
outputs = {
'E': E_far,
'H': H_far,
'dkx': kx[1]-kx[0],
'dky': ky[1]-ky[0],
'kx': kx,
'ky': ky,
'theta': theta,
'phi': phi,
}
return outputs
def far_to_nearfield(E_far: List[numpy.ndarray],
H_far: List[numpy.ndarray],
dkx: float,
dky: float,
padded_size: List[int] = None
) -> Dict[str]:
"""
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.
:param 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))
:param 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))
:param dkx: kx discretization, in units of wavelength.
:param dky: ky discretization, in units of wavelength.
:param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
Powers of 2 are most efficient for FFT computation.
Default is the smallest power of 2 larger than the input, for each axis.
:returns: Dict with keys
'E': E-field nearfield
'H': H-field nearfield
'dx', 'dy': spatial discretization, normalized to wavelength (dimensionless)
"""
if not len(E_far) == 2:
raise Exception('E_far must be a length-2 list of ndarrays')
if not len(H_far) == 2:
raise Exception('H_far must be a length-2 list of ndarrays')
s = E_far[0].shape
if not all(s == f.shape for f in E_far + H_far):
raise Exception('All fields must be the same shape!')
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)
k = 2 * pi
kxs = fftshift(fftfreq(s[0], 1/(s[0] * dkx)))
kys = fftshift(fftfreq(s[0], 1/(s[1] * dky)))
kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
kxy2 = kx * kx + ky * ky
kxy = numpy.sqrt(kxy2)
kz = numpy.sqrt(k * k - kxy2)
sin_th = ky / kxy
cos_th = kx / kxy
cos_phi = kz / k
sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
# Zero fields beyond valid (phi, theta)
invalid_ind = kxy2 >= k * k
theta[invalid_ind] = 0
phi[invalid_ind] = 0
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],
-0.5 * E_far[0]]
N = [L[1],
-L[0]]
En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th)/cos_phi,
-(-L[0] * cos_th + L[1] * cos_phi * sin_th)/cos_phi]
Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th)/cos_phi,
(-N[0] * cos_th + N[1] * cos_phi * sin_th)/cos_phi]
for i in range(2):
En_fft[i][cos_phi == 0] = 0
Hn_fft[i][cos_phi == 0] = 0
E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_size)) for Ei in En_fft]
H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_size)) for Hi in Hn_fft]
dx = 2 * pi / (s[0] * dkx)
dy = 2 * pi / (s[0] * dky)
outputs = {
'E': E_near,
'H': H_near,
'dx': dx,
'dy': dy,
}
return outputs

182
meanas/fdfd/functional.py Normal file
View file

@ -0,0 +1,182 @@
"""
Functional versions of many FDFD operators. These can be useful for performing
FDFD calculations without needing to construct large matrices in memory.
The functions generated here expect field inputs with shape (3, X, Y, Z),
e.g. E = [E_x, E_y, E_z] where each component has shape (X, Y, Z)
"""
from typing import List, Callable
import numpy
from . import dx_lists_t, field_t
__author__ = 'Jan Petykiewicz'
functional_matrix = Callable[[field_t], field_t]
def curl_h(dxes: dx_lists_t) -> functional_matrix:
"""
Curl operator for use with the H field.
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Function for taking the discretized curl of the H-field, F(H) -> curlH
"""
dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij')
def dh(f, ax):
return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax]
def ch_fun(h: field_t) -> field_t:
e = [dh(h[2], 1) - dh(h[1], 2),
dh(h[0], 2) - dh(h[2], 0),
dh(h[1], 0) - dh(h[0], 1)]
return e
return ch_fun
def curl_e(dxes: dx_lists_t) -> functional_matrix:
"""
Curl operator for use with the E field.
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Function for taking the discretized curl of the E-field, F(E) -> curlE
"""
dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij')
def de(f, ax):
return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax]
def ce_fun(e: field_t) -> field_t:
h = [de(e[2], 1) - de(e[1], 2),
de(e[0], 2) - de(e[2], 0),
de(e[1], 0) - de(e[0], 1)]
return h
return ce_fun
def e_full(omega: complex,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None
) -> functional_matrix:
"""
Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field,
with wave equation
(del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function implementing the wave operator A(E) -> E
"""
ch = curl_h(dxes)
ce = curl_e(dxes)
def op_1(e):
curls = ch(ce(e))
return [c - omega ** 2 * e * x for c, e, x in zip(curls, epsilon, e)]
def op_mu(e):
curls = ch([m * y for m, y in zip(mu, ce(e))])
return [c - omega ** 2 * p * x for c, p, x in zip(curls, epsilon, e)]
if numpy.any(numpy.equal(mu, None)):
return op_1
else:
return op_mu
def eh_full(omega: complex,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None
) -> functional_matrix:
"""
Wave operator for full (both E and H) field representation.
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function implementing the wave operator A(E, H) -> (E, H)
"""
ch = curl_h(dxes)
ce = curl_e(dxes)
def op_1(e, h):
return ([c - 1j * omega * p * x for c, p, x in zip(ch(h), epsilon, e)],
[c + 1j * omega * y for c, y in zip(ce(e), h)])
def op_mu(e, h):
return ([c - 1j * omega * p * x for c, p, x in zip(ch(h), epsilon, e)],
[c + 1j * omega * m * y for c, m, y in zip(ce(e), mu, h)])
if numpy.any(numpy.equal(mu, None)):
return op_1
else:
return op_mu
def e2h(omega: complex,
dxes: dx_lists_t,
mu: field_t = None,
) -> functional_matrix:
"""
Utility operator for converting the E field into the H field.
For use with e_full -- assumes that there is no magnetic current M.
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function for converting E to H
"""
A2 = curl_e(dxes)
def e2h_1_1(e):
return [y / (-1j * omega) for y in A2(e)]
def e2h_mu(e):
return [y / (-1j * omega * m) for y, m in zip(A2(e), mu)]
if numpy.any(numpy.equal(mu, None)):
return e2h_1_1
else:
return e2h_mu
def m2j(omega: complex,
dxes: dx_lists_t,
mu: field_t = None,
) -> functional_matrix:
"""
Utility operator for converting magnetic current (M) distribution
into equivalent electric current distribution (J).
For use with e.g. e_full().
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function for converting M to J
"""
ch = curl_h(dxes)
def m2j_mu(m):
m_mu = [m[k] / mu[k] for k in range[3]]
J = [Ji / (-1j * omega) for Ji in ch(m_mu)]
return J
def m2j_1(m):
J = [Ji / (-1j * omega) for Ji in ch(m)]
return J
if numpy.any(numpy.equal(mu, None)):
return m2j_1
else:
return m2j_mu

503
meanas/fdfd/operators.py Normal file
View file

@ -0,0 +1,503 @@
"""
Sparse matrix operators for use with electromagnetic wave equations.
These functions return sparse-matrix (scipy.sparse.spmatrix) representations of
a variety of operators, intended for use with E and H fields vectorized using the
meanas.vec() and .unvec() functions (column-major/Fortran ordering).
E- and H-field values are defined on a Yee cell; epsilon values should be calculated for
cells centered at each E component (mu at each H component).
Many of these functions require a 'dxes' parameter, of type meanas.dx_lists_type; see
the meanas.types submodule for details.
The following operators are included:
- E-only wave operator
- H-only wave operator
- EH wave operator
- Curl for use with E, H fields
- E to H conversion
- M to J conversion
- Poynting cross products
Also available:
- Circular shifts
- Discrete derivatives
- Averaging operators
- Cross product matrices
"""
from typing import List, Tuple
import numpy
import scipy.sparse as sparse
from . import vec, dx_lists_t, vfield_t
__author__ = 'Jan Petykiewicz'
def e_full(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
pec: vfield_t = None,
pmc: vfield_t = None,
) -> sparse.spmatrix:
"""
Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field,
with wave equation
(del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
To make this matrix symmetric, use the preconditions from e_full_preconditioners().
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Vectorized dielectric constant
:param mu: Vectorized magnetic permeability (default 1 everywhere).
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
as containing a perfect electrical conductor (PEC).
The PEC is applied per-field-component (ie, pec.size == epsilon.size)
:param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
as containing a perfect magnetic conductor (PMC).
The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
:return: Sparse matrix containing the wave operator
"""
ce = curl_e(dxes)
ch = curl_h(dxes)
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
else:
pe = sparse.diags(numpy.where(pec, 0, 1)) # Set pe to (not PEC)
if numpy.any(numpy.equal(pmc, None)):
pm = sparse.eye(epsilon.size)
else:
pm = sparse.diags(numpy.where(pmc, 0, 1)) # set pm to (not PMC)
e = sparse.diags(epsilon)
if numpy.any(numpy.equal(mu, None)):
m_div = sparse.eye(epsilon.size)
else:
m_div = sparse.diags(1 / mu)
op = pe @ (ch @ pm @ m_div @ ce - omega**2 * e) @ pe
return op
def e_full_preconditioners(dxes: dx_lists_t
) -> Tuple[sparse.spmatrix, sparse.spmatrix]:
"""
Left and right preconditioners (Pl, Pr) for symmetrizing the e_full wave operator.
The preconditioned matrix A_symm = (Pl @ A @ Pr) is complex-symmetric
(non-Hermitian unless there is no loss or PMLs).
The preconditioner matrices are diagonal and complex, with Pr = 1 / Pl
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Preconditioner matrices (Pl, Pr)
"""
p_squared = [dxes[0][0][:, None, None] * dxes[1][1][None, :, None] * dxes[1][2][None, None, :],
dxes[1][0][:, None, None] * dxes[0][1][None, :, None] * dxes[1][2][None, None, :],
dxes[1][0][:, None, None] * dxes[1][1][None, :, None] * dxes[0][2][None, None, :]]
p_vector = numpy.sqrt(vec(p_squared))
P_left = sparse.diags(p_vector)
P_right = sparse.diags(1 / p_vector)
return P_left, P_right
def h_full(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
pec: vfield_t = None,
pmc: vfield_t = None,
) -> sparse.spmatrix:
"""
Wave operator del x (1/epsilon * del x) - omega**2 * mu, for use with H-field,
with wave equation
(del x (1/epsilon * del x) - omega**2 * mu) H = i * omega * M
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Vectorized dielectric constant
:param mu: Vectorized magnetic permeability (default 1 everywhere)
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
as containing a perfect electrical conductor (PEC).
The PEC is applied per-field-component (ie, pec.size == epsilon.size)
:param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
as containing a perfect magnetic conductor (PMC).
The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
:return: Sparse matrix containing the wave operator
"""
ec = curl_e(dxes)
hc = curl_h(dxes)
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
else:
pe = sparse.diags(numpy.where(pec, 0, 1)) # set pe to (not PEC)
if numpy.any(numpy.equal(pmc, None)):
pm = sparse.eye(epsilon.size)
else:
pm = sparse.diags(numpy.where(pmc, 0, 1)) # Set pe to (not PMC)
e_div = sparse.diags(1 / epsilon)
if mu is None:
m = sparse.eye(epsilon.size)
else:
m = sparse.diags(mu)
A = pm @ (ec @ pe @ e_div @ hc - omega**2 * m) @ pm
return A
def eh_full(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
pec: vfield_t = None,
pmc: vfield_t = None
) -> sparse.spmatrix:
"""
Wave operator for [E, H] field representation. This operator implements Maxwell's
equations without cancelling out either E or H. The operator is
[[-i * omega * epsilon, del x],
[del x, i * omega * mu]]
for use with a field vector of the form hstack(vec(E), vec(H)).
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Vectorized dielectric constant
:param mu: Vectorized magnetic permeability (default 1 everywhere)
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
as containing a perfect electrical conductor (PEC).
The PEC is applied per-field-component (i.e., pec.size == epsilon.size)
:param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
as containing a perfect magnetic conductor (PMC).
The PMC is applied per-field-component (i.e., pmc.size == epsilon.size)
:return: Sparse matrix containing the wave operator
"""
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
else:
pe = sparse.diags(numpy.where(pec, 0, 1)) # set pe to (not PEC)
if numpy.any(numpy.equal(pmc, None)):
pm = sparse.eye(epsilon.size)
else:
pm = sparse.diags(numpy.where(pmc, 0, 1)) # set pm to (not PMC)
iwe = pe @ (1j * omega * sparse.diags(epsilon)) @ pe
iwm = 1j * omega
if not numpy.any(numpy.equal(mu, None)):
iwm *= sparse.diags(mu)
iwm = pm @ iwm @ pm
A1 = pe @ curl_h(dxes) @ pm
A2 = pm @ curl_e(dxes) @ pe
A = sparse.bmat([[-iwe, A1],
[A2, iwm]])
return A
def curl_h(dxes: dx_lists_t) -> sparse.spmatrix:
"""
Curl operator for use with the H field.
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Sparse matrix for taking the discretized curl of the H-field
"""
return cross(deriv_back(dxes[1]))
def curl_e(dxes: dx_lists_t) -> sparse.spmatrix:
"""
Curl operator for use with the E field.
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Sparse matrix for taking the discretized curl of the E-field
"""
return cross(deriv_forward(dxes[0]))
def e2h(omega: complex,
dxes: dx_lists_t,
mu: vfield_t = None,
pmc: vfield_t = None,
) -> sparse.spmatrix:
"""
Utility operator for converting the E field into the H field.
For use with e_full -- assumes that there is no magnetic current M.
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param mu: Vectorized magnetic permeability (default 1 everywhere)
:param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
as containing a perfect magnetic conductor (PMC).
The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
:return: Sparse matrix for converting E to H
"""
op = curl_e(dxes) / (-1j * omega)
if not numpy.any(numpy.equal(mu, None)):
op = sparse.diags(1 / mu) @ op
if not numpy.any(numpy.equal(pmc, None)):
op = sparse.diags(numpy.where(pmc, 0, 1)) @ op
return op
def m2j(omega: complex,
dxes: dx_lists_t,
mu: vfield_t = None
) -> sparse.spmatrix:
"""
Utility operator for converting M field into J.
Converts a magnetic current M into an electric current J.
For use with eg. e_full.
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param mu: Vectorized magnetic permeability (default 1 everywhere)
:return: Sparse matrix for converting E to H
"""
op = curl_h(dxes) / (1j * omega)
if not numpy.any(numpy.equal(mu, None)):
op = op @ sparse.diags(1 / mu)
return op
def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
"""
Utility operator for performing a circular shift along a specified axis by a
specified number of elements.
:param axis: Axis to shift along. x=0, y=1, z=2
:param shape: Shape of the grid being shifted
:param shift_distance: Number of cells to shift by. May be negative. Default 1.
:return: Sparse matrix for performing the circular shift
"""
if len(shape) not in (2, 3):
raise Exception('Invalid shape: {}'.format(shape))
if axis not in range(len(shape)):
raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
shifts = [abs(shift_distance) if a == axis else 0 for a in range(3)]
shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)]
ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
n = numpy.prod(shape)
i_ind = numpy.arange(n)
j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
d = sparse.csr_matrix(vij, shape=(n, n))
if shift_distance < 0:
d = d.T
return d
def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
"""
Utility operator for performing an n-element shift along a specified axis, with mirror
boundary conditions applied to the cells beyond the receding edge.
:param axis: Axis to shift along. x=0, y=1, z=2
:param shape: Shape of the grid being shifted
:param shift_distance: Number of cells to shift by. May be negative. Default 1.
:return: Sparse matrix for performing the circular shift
"""
if len(shape) not in (2, 3):
raise Exception('Invalid shape: {}'.format(shape))
if axis not in range(len(shape)):
raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
if shift_distance >= shape[axis]:
raise Exception('Shift ({}) is too large for axis {} of size {}'.format(
shift_distance, axis, shape[axis]))
def mirrored_range(n, s):
v = numpy.arange(n) + s
v = numpy.where(v >= n, 2 * n - v - 1, v)
v = numpy.where(v < 0, - 1 - v, v)
return v
shifts = [shift_distance if a == axis else 0 for a in range(3)]
shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts)]
ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
n = numpy.prod(shape)
i_ind = numpy.arange(n)
j_ind = ijk[0] + ijk[1] * shape[0]
if len(shape) == 3:
j_ind += ijk[2] * shape[0] * shape[1]
vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
d = sparse.csr_matrix(vij, shape=(n, n))
return d
def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]:
"""
Utility operators for taking discretized derivatives (forward variant).
:param dx_e: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...].
:return: List of operators for taking forward derivatives along each axis.
"""
shape = [s.size for s in dx_e]
n = numpy.prod(shape)
dx_e_expanded = numpy.meshgrid(*dx_e, indexing='ij')
def deriv(axis):
return rotation(axis, shape, 1) - sparse.eye(n)
Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a)
for a, dx in enumerate(dx_e_expanded)]
return Ds
def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]:
"""
Utility operators for taking discretized derivatives (backward variant).
:param dx_h: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...].
:return: List of operators for taking forward derivatives along each axis.
"""
shape = [s.size for s in dx_h]
n = numpy.prod(shape)
dx_h_expanded = numpy.meshgrid(*dx_h, indexing='ij')
def deriv(axis):
return rotation(axis, shape, -1) - sparse.eye(n)
Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a)
for a, dx in enumerate(dx_h_expanded)]
return Ds
def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix:
"""
Cross product operator
:param B: List [Bx, By, Bz] of sparse matrices corresponding to the x, y, z
portions of the operator on the left side of the cross product.
:return: Sparse matrix corresponding to (B x), where x is the cross product
"""
n = B[0].shape[0]
zero = sparse.csr_matrix((n, n))
return sparse.bmat([[zero, -B[2], B[1]],
[B[2], zero, -B[0]],
[-B[1], B[0], zero]])
def vec_cross(b: vfield_t) -> sparse.spmatrix:
"""
Vector cross product operator
:param b: Vector on the left side of the cross product
:return: Sparse matrix corresponding to (b x), where x is the cross product
"""
B = [sparse.diags(c) for c in numpy.split(b, 3)]
return cross(B)
def avgf(axis: int, shape: List[int]) -> sparse.spmatrix:
"""
Forward average operator (x4 = (x4 + x5) / 2)
:param axis: Axis to average along (x=0, y=1, z=2)
:param shape: Shape of the grid to average
:return: Sparse matrix for forward average operation
"""
if len(shape) not in (2, 3):
raise Exception('Invalid shape: {}'.format(shape))
n = numpy.prod(shape)
return 0.5 * (sparse.eye(n) + rotation(axis, shape))
def avgb(axis: int, shape: List[int]) -> sparse.spmatrix:
"""
Backward average operator (x4 = (x4 + x3) / 2)
:param axis: Axis to average along (x=0, y=1, z=2)
:param shape: Shape of the grid to average
:return: Sparse matrix for backward average operation
"""
return avgf(axis, shape).T
def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
"""
Operator for computing the Poynting vector, containing the (E x) portion of the Poynting vector.
:param e: Vectorized E-field for the ExH cross product
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Sparse matrix containing (E x) portion of Poynting cross product
"""
shape = [len(dx) for dx in dxes[0]]
fx, fy, fz = [avgf(i, shape) for i in range(3)]
bx, by, bz = [avgb(i, shape) for i in range(3)]
dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C'))
for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
Ex, Ey, Ez = [sparse.diags(ei * da) for ei, da in zip(numpy.split(e, 3), dxag)]
n = numpy.prod(shape)
zero = sparse.csr_matrix((n, n))
P = sparse.bmat(
[[ zero, -fx @ Ez @ bz @ dbgy, fx @ Ey @ by @ dbgz],
[ fy @ Ez @ bz @ dbgx, zero, -fy @ Ex @ bx @ dbgz],
[-fz @ Ey @ by @ dbgx, fz @ Ex @ bx @ dbgy, zero]])
return P
def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
"""
Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector.
:param h: Vectorized H-field for the HxE cross product
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:return: Sparse matrix containing (H x) portion of Poynting cross product
"""
shape = [len(dx) for dx in dxes[0]]
fx, fy, fz = [avgf(i, shape) for i in range(3)]
bx, by, bz = [avgb(i, shape) for i in range(3)]
dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
dagx, dagy, dagz = [sparse.diags(dx.ravel(order='C'))
for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
Hx, Hy, Hz = [sparse.diags(hi * db) for hi, db in zip(numpy.split(h, 3), dxbg)]
n = numpy.prod(shape)
zero = sparse.csr_matrix((n, n))
P = sparse.bmat(
[[ zero, -by @ Hz @ fx @ dagy, bz @ Hy @ fx @ dagz],
[ bx @ Hz @ fy @ dagx, zero, -bz @ Hx @ fy @ dagz],
[-bx @ Hy @ fz @ dagx, by @ Hx @ fz @ dagy, zero]])
return P

169
meanas/fdfd/scpml.py Normal file
View file

@ -0,0 +1,169 @@
"""
Functions for creating stretched coordinate PMLs.
"""
from typing import List, Callable
import numpy
__author__ = 'Jan Petykiewicz'
s_function_type = Callable[[float], float]
def prepare_s_function(ln_R: float = -16,
m: float = 4
) -> s_function_type:
"""
Create an s_function to pass to the SCPML functions. This is used when you would like to
customize the PML parameters.
:param ln_R: Natural logarithm of the desired reflectance
:param m: Polynomial order for the PML (imaginary part increases as distance ** m)
:return: An s_function, which takes an ndarray (distances) and returns an ndarray (complex part
of the cell width; needs to be divided by sqrt(epilon_effective) * real(omega))
before use.
"""
def s_factor(distance: numpy.ndarray) -> numpy.ndarray:
s_max = (m + 1) * ln_R / 2 # / 2 because we assume periodic boundaries
return s_max * (distance ** m)
return s_factor
def uniform_grid_scpml(shape: numpy.ndarray or List[int],
thicknesses: numpy.ndarray or List[int],
omega: float,
epsilon_effective: float = 1.0,
s_function: s_function_type = None,
) -> dx_lists_t:
"""
Create dx arrays for a uniform grid with a cell width of 1 and a pml.
If you want something more fine-grained, check out stretch_with_scpml(...).
:param shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick)
:param thicknesses: [th_x, th_y, th_z] Thickness of the PML in each direction.
Both polarities are added.
Each th_ of pml is applied twice, once on each edge of the grid along the given axis.
th_* may be zero, in which case no pml is added.
:param omega: Angular frequency for the simulation
:param epsilon_effective: Effective epsilon of the PML. Match this to the material
at the edge of your grid.
Default 1.
:param s_function: s_function created by prepare_s_function(...), allowing
customization of pml parameters.
Default uses prepare_s_function() with no parameters.
:return: Complex cell widths (dx_lists)
"""
if s_function is None:
s_function = prepare_s_function()
# Normalized distance to nearest boundary
def l(u, n, t):
return ((t - u).clip(0) + (u - (n - t)).clip(0)) / t
dx_a = [numpy.array(numpy.inf)] * 3
dx_b = [numpy.array(numpy.inf)] * 3
# divide by this to adjust for epsilon_effective and omega
s_correction = numpy.sqrt(epsilon_effective) * numpy.real(omega)
for k, th in enumerate(thicknesses):
s = shape[k]
if th > 0:
sr = numpy.arange(s)
dx_a[k] = 1 + 1j * s_function(l(sr, s, th)) / s_correction
dx_b[k] = 1 + 1j * s_function(l(sr+0.5, s, th)) / s_correction
else:
dx_a[k] = numpy.ones((s,))
dx_b[k] = numpy.ones((s,))
return [dx_a, dx_b]
def stretch_with_scpml(dxes: dx_lists_t,
axis: int,
polarity: int,
omega: float,
epsilon_effective: float = 1.0,
thickness: int = 10,
s_function: s_function_type = None,
) -> dx_lists_t:
"""
Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis.
:param dxes: dx_tuple with coordinates to stretch
:param axis: axis to stretch (0=x, 1=y, 2=z)
:param polarity: direction to stretch (-1 for -ve, +1 for +ve)
:param omega: Angular frequency for the simulation
:param epsilon_effective: Effective epsilon of the PML. Match this to the material at the
edge of your grid. Default 1.
:param thickness: number of cells to use for pml (default 10)
:param s_function: s_function created by prepare_s_function(...), allowing customization
of pml parameters. Default uses prepare_s_function() with no parameters.
:return: Complex cell widths
"""
if s_function is None:
s_function = prepare_s_function()
dx_ai = dxes[0][axis].astype(complex)
dx_bi = dxes[1][axis].astype(complex)
pos = numpy.hstack((0, dx_ai.cumsum()))
pos_a = (pos[:-1] + pos[1:]) / 2
pos_b = pos[:-1]
# divide by this to adjust for epsilon_effective and omega
s_correction = numpy.sqrt(epsilon_effective) * numpy.real(omega)
if polarity > 0:
# front pml
bound = pos[thickness]
d = bound - pos[0]
def l_d(x):
return (bound - x) / (bound - pos[0])
slc = slice(thickness)
else:
# back pml
bound = pos[-thickness - 1]
d = pos[-1] - bound
def l_d(x):
return (x - bound) / (pos[-1] - bound)
if thickness == 0:
slc = slice(None)
else:
slc = slice(-thickness, None)
dx_ai[slc] *= 1 + 1j * s_function(l_d(pos_a[slc])) / d / s_correction
dx_bi[slc] *= 1 + 1j * s_function(l_d(pos_b[slc])) / d / s_correction
dxes[0][axis] = dx_ai
dxes[1][axis] = dx_bi
return dxes
def generate_periodic_dx(pos: List[numpy.ndarray]) -> dx_lists_t:
"""
Given a list of 3 ndarrays cell centers, creates the cell width parameters for a periodic grid.
:param pos: List of 3 ndarrays of cell centers
:return: (dx_a, dx_b) cell widths (no pml)
"""
if len(pos) != 3:
raise Exception('Must have len(pos) == 3')
dx_a = [numpy.array(numpy.inf)] * 3
dx_b = [numpy.array(numpy.inf)] * 3
for i, p_orig in enumerate(pos):
p = numpy.array(p_orig, dtype=float)
if p.size != 1:
p_shifted = numpy.hstack((p[1:], p[-1] + (p[1] - p[0])))
dx_a[i] = numpy.diff(p)
dx_b[i] = numpy.diff((p + p_shifted) / 2)
return dx_a, dx_b

118
meanas/fdfd/solvers.py Normal file
View file

@ -0,0 +1,118 @@
"""
Solvers for FDFD problems.
"""
from typing import List, Callable, Dict, Any
import logging
import numpy
from numpy.linalg import norm
import scipy.sparse.linalg
from . import operators
logger = logging.getLogger(__name__)
def _scipy_qmr(A: scipy.sparse.csr_matrix,
b: numpy.ndarray,
**kwargs
) -> numpy.ndarray:
"""
Wrapper for scipy.sparse.linalg.qmr
:param A: Sparse matrix
:param b: Right-hand-side vector
:param kwargs: Passed as **kwargs to the wrapped function
:return: Guess for solution (returned even if didn't converge)
"""
'''
Report on our progress
'''
iter = 0
def log_residual(xk):
nonlocal iter
iter += 1
if iter % 100 == 0:
logger.info('Solver residual at iteration {} : {}'.format(iter, norm(A @ xk - b)))
if 'callback' in kwargs:
def augmented_callback(xk):
log_residual(xk)
kwargs['callback'](xk)
kwargs['callback'] = augmented_callback
else:
kwargs['callback'] = log_residual
'''
Run the actual solve
'''
x, _ = scipy.sparse.linalg.qmr(A, b, **kwargs)
return x
def generic(omega: complex,
dxes: List[List[numpy.ndarray]],
J: numpy.ndarray,
epsilon: numpy.ndarray,
mu: numpy.ndarray = None,
pec: numpy.ndarray = None,
pmc: numpy.ndarray = None,
adjoint: bool = False,
matrix_solver: Callable[..., numpy.ndarray] = _scipy_qmr,
matrix_solver_opts: Dict[str, Any] = None,
) -> numpy.ndarray:
"""
Conjugate gradient FDFD solver using CSR sparse matrices.
All ndarray arguments should be 1D array, as returned by meanas.vec().
:param omega: Complex frequency to solve at.
:param dxes: [[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]] (complex cell sizes)
:param J: Electric current distribution (at E-field locations)
:param epsilon: Dielectric constant distribution (at E-field locations)
:param mu: Magnetic permeability distribution (at H-field locations)
:param pec: Perfect electric conductor distribution
(at E-field locations; non-zero value indicates PEC is present)
:param pmc: Perfect magnetic conductor distribution
(at H-field locations; non-zero value indicates PMC is present)
:param adjoint: If true, solves the adjoint problem.
:param matrix_solver: Called as matrix_solver(A, b, **matrix_solver_opts) -> x
Where A: scipy.sparse.csr_matrix
b: numpy.ndarray
x: numpy.ndarray
Default is a wrapped version of scipy.sparse.linalg.qmr()
which doesn't return convergence info and logs the residual
every 100 iterations.
:param matrix_solver_opts: Passed as kwargs to matrix_solver(...)
:return: E-field which solves the system.
"""
if matrix_solver_opts is None:
matrix_solver_opts = dict()
b0 = -1j * omega * J
A0 = operators.e_full(omega, dxes, epsilon=epsilon, mu=mu, pec=pec, pmc=pmc)
Pl, Pr = operators.e_full_preconditioners(dxes)
if adjoint:
A = (Pl @ A0 @ Pr).H
b = Pr.H @ b0
else:
A = Pl @ A0 @ Pr
b = Pl @ b0
x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
if adjoint:
x0 = Pl.H @ x
else:
x0 = Pr @ x
return x0

377
meanas/fdfd/waveguide.py Normal file
View file

@ -0,0 +1,377 @@
"""
Various operators and helper functions for solving for waveguide modes.
Assuming a z-dependence of the from exp(-i * wavenumber * z), we can simplify Maxwell's
equations in the absence of sources to the form
A @ [H_x, H_y] = wavenumber**2 * [H_x, H_y]
with A =
omega**2 * epsilon * mu +
epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
[[Dx], [Dy]] / mu * [Dx, Dy] * mu
which is the form used in this file.
As the z-dependence is known, all the functions in this file assume a 2D grid
(ie. dxes = [[[dx_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dx_h_0, ...], [dy_h_0, ...]]])
with propagation along the z axis.
"""
from typing import List, Tuple
import numpy
from numpy.linalg import norm
import scipy.sparse as sparse
from . import vec, unvec, dx_lists_t, field_t, vfield_t
from . import operators
__author__ = 'Jan Petykiewicz'
def operator(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
) -> sparse.spmatrix:
"""
Waveguide operator of the form
omega**2 * epsilon * mu +
epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
[[Dx], [Dy]] / mu * [Dx, Dy] * mu
for use with a field vector of the form [H_x, H_y].
This operator can be used to form an eigenvalue problem of the form
A @ [H_x, H_y] = wavenumber**2 * [H_x, H_y]
which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * z)
z-dependence is assumed for the fields).
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Sparse matrix representation of the operator
"""
if numpy.any(numpy.equal(mu, None)):
mu = numpy.ones_like(epsilon)
Dfx, Dfy = operators.deriv_forward(dxes[0])
Dbx, Dby = operators.deriv_back(dxes[1])
eps_parts = numpy.split(epsilon, 3)
eps_yx = sparse.diags(numpy.hstack((eps_parts[1], eps_parts[0])))
eps_z_inv = sparse.diags(1 / eps_parts[2])
mu_parts = numpy.split(mu, 3)
mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
mu_z_inv = sparse.diags(1 / mu_parts[2])
op = omega ** 2 * eps_yx @ mu_xy + \
eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \
sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy
return op
def normalized_fields(v: numpy.ndarray,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> Tuple[vfield_t, vfield_t]:
"""
Given a vector v containing the vectorized H_x and H_y fields,
returns normalized, vectorized E and H fields for the system.
:param v: Vector containing H_x and H_y fields
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Normalized, vectorized (e, h) containing all vector components.
"""
e = v2e(v, wavenumber, omega, dxes, epsilon, mu=mu)
h = v2h(v, wavenumber, dxes, mu=mu)
shape = [s.size for s in dxes[0]]
dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)]
E = unvec(e, shape)
H = unvec(h, shape)
S1 = E[0] * numpy.roll(numpy.conj(H[1]), 1, axis=0) * dxes_real[0][1] * dxes_real[1][0]
S2 = E[1] * numpy.roll(numpy.conj(H[0]), 1, axis=1) * dxes_real[0][0] * dxes_real[1][1]
S = 0.25 * ((S1 + numpy.roll(S1, 1, axis=0)) -
(S2 + numpy.roll(S2, 1, axis=1)))
P = 0.5 * numpy.real(S.sum())
assert P > 0, 'Found a mode propagating in the wrong direction! P={}'.format(P)
energy = epsilon * e.conj() * e
norm_amplitude = 1 / numpy.sqrt(P)
norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric
# Try to break symmetry to assign a consistent sign [experimental]
E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum())
logger.debug('norm_angle = {}'.format(norm_angle))
logger.debug('norm_sign = {}'.format(sign)
norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
e *= norm_factor
h *= norm_factor
return e, h
def v2h(v: numpy.ndarray,
wavenumber: complex,
dxes: dx_lists_t,
mu: vfield_t = None
) -> vfield_t:
"""
Given a vector v containing the vectorized H_x and H_y fields,
returns a vectorized H including all three H components.
:param v: Vector containing H_x and H_y fields
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Vectorized H field with all vector components
"""
Dfx, Dfy = operators.deriv_forward(dxes[0])
op = sparse.hstack((Dfx, Dfy))
if not numpy.any(numpy.equal(mu, None)):
mu_parts = numpy.split(mu, 3)
mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
mu_z_inv = sparse.diags(1 / mu_parts[2])
op = mu_z_inv @ op @ mu_xy
w = op @ v / (1j * wavenumber)
return numpy.hstack((v, w)).flatten()
def v2e(v: numpy.ndarray,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> vfield_t:
"""
Given a vector v containing the vectorized H_x and H_y fields,
returns a vectorized E containing all three E components
:param v: Vector containing H_x and H_y fields
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Vectorized E field with all vector components.
"""
h2eop = h2e(wavenumber, omega, dxes, epsilon)
return h2eop @ v2h(v, wavenumber, dxes, mu)
def e2h(wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
mu: vfield_t = None
) -> sparse.spmatrix:
"""
Returns an operator which, when applied to a vectorized E eigenfield, produces
the vectorized H eigenfield.
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Sparse matrix representation of the operator
"""
op = curl_e(wavenumber, dxes) / (-1j * omega)
if not numpy.any(numpy.equal(mu, None)):
op = sparse.diags(1 / mu) @ op
return op
def h2e(wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t
) -> sparse.spmatrix:
"""
Returns an operator which, when applied to a vectorized H eigenfield, produces
the vectorized E eigenfield.
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:return: Sparse matrix representation of the operator
"""
op = sparse.diags(1 / (1j * omega * epsilon)) @ curl_h(wavenumber, dxes)
return op
def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
"""
Discretized curl operator for use with the waveguide E field.
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:return: Sparse matrix representation of the operator
"""
n = 1
for d in dxes[0]:
n *= len(d)
Bz = -1j * wavenumber * sparse.eye(n)
Dfx, Dfy = operators.deriv_forward(dxes[0])
return operators.cross([Dfx, Dfy, Bz])
def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
"""
Discretized curl operator for use with the waveguide H field.
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:return: Sparse matrix representation of the operator
"""
n = 1
for d in dxes[1]:
n *= len(d)
Bz = -1j * wavenumber * sparse.eye(n)
Dbx, Dby = operators.deriv_back(dxes[1])
return operators.cross([Dbx, Dby, Bz])
def h_err(h: vfield_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> float:
"""
Calculates the relative error in the H field
:param h: Vectorized H field
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Relative error norm(OP @ h) / norm(h)
"""
ce = curl_e(wavenumber, dxes)
ch = curl_h(wavenumber, dxes)
eps_inv = sparse.diags(1 / epsilon)
if numpy.any(numpy.equal(mu, None)):
op = ce @ eps_inv @ ch @ h - omega ** 2 * h
else:
op = ce @ eps_inv @ ch @ h - omega ** 2 * (mu * h)
return norm(op) / norm(h)
def e_err(e: vfield_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> float:
"""
Calculates the relative error in the E field
:param e: Vectorized E field
:param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Relative error norm(OP @ e) / norm(e)
"""
ce = curl_e(wavenumber, dxes)
ch = curl_h(wavenumber, dxes)
if numpy.any(numpy.equal(mu, None)):
op = ch @ ce @ e - omega ** 2 * (epsilon * e)
else:
mu_inv = sparse.diags(1 / mu)
op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e)
return norm(op) / norm(e)
def cylindrical_operator(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
r0: float,
) -> sparse.spmatrix:
"""
Cylindrical coordinate waveguide operator of the form
TODO
for use with a field vector of the form [E_r, E_y].
This operator can be used to form an eigenvalue problem of the form
A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y]
which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * theta)
theta-dependence is assumed for the fields).
:param omega: The angular frequency of the system
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
:param epsilon: Vectorized dielectric constant grid
:param r0: Radius of curvature for the simulation. This should be the minimum value of
r within the simulation domain.
:return: Sparse matrix representation of the operator
"""
Dfx, Dfy = operators.deriv_forward(dxes[0])
Dbx, Dby = operators.deriv_back(dxes[1])
rx = r0 + numpy.cumsum(dxes[0][0])
ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0])
tx = rx/r0
ty = ry/r0
Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1)))
Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1)))
eps_parts = numpy.split(epsilon, 3)
eps_x = sparse.diags(eps_parts[0])
eps_y = sparse.diags(eps_parts[1])
eps_z_inv = sparse.diags(1 / eps_parts[2])
pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby))
pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx))
a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy
a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx
b0 = Dbx @ Ty @ Dfy
b1 = Dby @ Ty @ Dfx
diag = sparse.block_diag
op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
- (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1))
return op

475
meanas/fdfd/waveguide_mode.py Normal file
View file

@ -0,0 +1,475 @@
from typing import Dict, List
import numpy
import scipy.sparse as sparse
from . import vec, unvec, dx_lists_t, vfield_t, field_t
from . import operators, waveguide, functional
from .eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
def solve_waveguide_mode_2d(mode_number: int,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
wavenumber_correction: bool = True,
) -> Dict[str, complex or field_t]:
"""
Given a 2d region, attempts to solve for the eigenmode with the specified mode number.
:param mode_number: Number of the mode, 0-indexed.
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:param wavenumber_correction: Whether to correct the wavenumber to
account for numerical dispersion (default True)
:return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
"""
'''
Solve for the largest-magnitude eigenvalue of the real operator
'''
dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
A_r = waveguide.operator(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
eigvals, eigvecs = signed_eigensolve(A_r, mode_number+3)
v = eigvecs[:, -(mode_number + 1)]
'''
Now solve for the eigenvector of the full operator, using the real operator's
eigenvector as an initial guess for Rayleigh quotient iteration.
'''
A = waveguide.operator(omega, dxes, epsilon, mu)
eigval, v = rayleigh_quotient_iteration(A, v)
# Calculate the wave-vector (force the real part to be positive)
wavenumber = numpy.sqrt(eigval)
wavenumber *= numpy.sign(numpy.real(wavenumber))
e, h = waveguide.normalized_fields(v, wavenumber, omega, dxes, epsilon, mu)
'''
Perform correction on wavenumber to account for numerical dispersion.
See Numerical Dispersion in Taflove's FDTD book.
This correction term reduces the error in emitted power, but additional
error is introduced into the E_err and H_err terms. This effect becomes
more pronounced as the wavenumber increases.
'''
if wavenumber_correction:
dx_mean = (numpy.hstack(dxes[0]) + numpy.hstack(dxes[1])).mean() / 2 #TODO figure out what dx to use here
wavenumber -= 2 * numpy.sin(numpy.real(wavenumber * dx_mean / 2)) / dx_mean - numpy.real(wavenumber)
shape = [d.size for d in dxes[0]]
fields = {
'wavenumber': wavenumber,
'E': unvec(e, shape),
'H': unvec(h, shape),
}
return fields
def solve_waveguide_mode(mode_number: int,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
epsilon: field_t,
mu: field_t = None,
wavenumber_correction: bool = True
) -> Dict[str, complex or numpy.ndarray]:
"""
Given a 3D grid, selects a slice from the grid and attempts to
solve for an eigenmode propagating through that slice.
:param mode_number: Number of the mode, 0-indexed
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param axis: Propagation axis (0=x, 1=y, 2=z)
:param polarity: Propagation direction (+1 for +ve, -1 for -ve)
:param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
as the waveguide cross-section. slices[axis] should select only one
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:param wavenumber_correction: Whether to correct the wavenumber to
account for numerical dispersion (default True)
:return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
"""
if mu is None:
mu = numpy.ones_like(epsilon)
slices = tuple(slices)
'''
Solve the 2D problem in the specified plane
'''
# Define rotation to set z as propagation direction
order = numpy.roll(range(3), 2 - axis)
reverse_order = numpy.roll(range(3), axis - 2)
# Reduce to 2D and solve the 2D problem
args_2d = {
'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes],
'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]),
'mu': vec([mu[i][slices].transpose(order) for i in order]),
'wavenumber_correction': wavenumber_correction,
}
fields_2d = solve_waveguide_mode_2d(mode_number, omega=omega, **args_2d)
'''
Apply corrections and expand to 3D
'''
# Scale based on dx in propagation direction
dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes])
# Adjust for propagation direction
fields_2d['E'][2] *= polarity
fields_2d['H'][2] *= polarity
# Apply phase shift to H-field
d_prop = 0.5 * sum(dxab_forward)
fields_2d['H'] *= numpy.exp(-polarity * 1j * 0.5 * fields_2d['wavenumber'] * d_prop)
# Expand E, H to full epsilon space we were given
E = numpy.zeros_like(epsilon, dtype=complex)
H = numpy.zeros_like(epsilon, dtype=complex)
for a, o in enumerate(reverse_order):
E[(a, *slices)] = fields_2d['E'][o][:, :, None].transpose(reverse_order)
H[(a, *slices)] = fields_2d['H'][o][:, :, None].transpose(reverse_order)
results = {
'wavenumber': fields_2d['wavenumber'],
'H': H,
'E': E,
}
return results
def compute_source(E: field_t,
H: field_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
mu: field_t = None,
) -> field_t:
"""
Given an eigenmode obtained by solve_waveguide_mode, returns the current source distribution
necessary to position a unidirectional source at the slice location.
:param E: E-field of the mode
:param H: H-field of the mode (advanced by half of a Yee cell from E)
:param wavenumber: Wavenumber of the mode
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param axis: Propagation axis (0=x, 1=y, 2=z)
:param polarity: Propagation direction (+1 for +ve, -1 for -ve)
:param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
as the waveguide cross-section. slices[axis] should select only one
:param mu: Magnetic permeability (default 1 everywhere)
:return: J distribution for the unidirectional source
"""
if mu is None:
mu = numpy.ones(3)
J = numpy.zeros_like(E, dtype=complex)
M = numpy.zeros_like(E, dtype=complex)
src_order = numpy.roll(range(3), -axis)
exp_iphi = numpy.exp(1j * polarity * wavenumber * dxes[1][axis][slices[axis]])
J[src_order[1]] = +exp_iphi * H[src_order[2]] * polarity
J[src_order[2]] = -exp_iphi * H[src_order[1]] * polarity
rollby = -1 if polarity > 0 else 0
M[src_order[1]] = +numpy.roll(E[src_order[2]], rollby, axis=axis)
M[src_order[2]] = -numpy.roll(E[src_order[1]], rollby, axis=axis)
m2j = functional.m2j(omega, dxes, mu)
Jm = m2j(M)
Jtot = J + Jm
return Jtot
def compute_overlap_e(E: field_t,
H: field_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
mu: field_t = None,
) -> field_t:
"""
Given an eigenmode obtained by solve_waveguide_mode, calculates overlap_e for the
mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
[assumes reflection symmetry].
overlap_e makes use of the e2h operator to collapse the above expression into
(vec(E) @ vec(overlap_e)), allowing for simple calculation of the mode overlap.
:param E: E-field of the mode
:param H: H-field of the mode (advanced by half of a Yee cell from E)
:param wavenumber: Wavenumber of the mode
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
:param axis: Propagation axis (0=x, 1=y, 2=z)
:param polarity: Propagation direction (+1 for +ve, -1 for -ve)
:param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
as the waveguide cross-section. slices[axis] should select only one
:param mu: Magnetic permeability (default 1 everywhere)
:return: overlap_e for calculating the mode overlap
"""
slices = tuple(slices)
cross_plane = [slice(None)] * 4
cross_plane[axis + 1] = slices[axis]
cross_plane = tuple(cross_plane)
# Determine phase factors for parallel slices
a_shape = numpy.roll([-1, 1, 1], axis)
a_E = numpy.real(dxes[0][axis]).cumsum()
a_H = numpy.real(dxes[1][axis]).cumsum()
iphi = -polarity * 1j * wavenumber
phase_E = numpy.exp(iphi * (a_E - a_E[slices[axis]])).reshape(a_shape)
phase_H = numpy.exp(iphi * (a_H - a_H[slices[axis]])).reshape(a_shape)
# Expand our slice to the entire grid using the calculated phase factors
Ee = phase_E * E[cross_plane]
He = phase_H * H[cross_plane]
# Write out the operator product for the mode orthogonality integral
domain = numpy.zeros_like(E[0], dtype=int)
domain[slices] = 1
npts = E[0].size
dn = numpy.zeros(npts * 3, dtype=int)
dn[0:npts] = 1
dn = numpy.roll(dn, npts * axis)
e2h = operators.e2h(omega, dxes, mu)
ds = sparse.diags(vec([domain]*3))
h_cross_ = operators.poynting_h_cross(vec(He), dxes)
e_cross_ = operators.poynting_e_cross(vec(Ee), dxes)
overlap_e = dn @ ds @ (-h_cross_ + e_cross_ @ e2h)
# Normalize
dx_forward = dxes[0][axis][slices[axis]]
norm_factor = numpy.abs(overlap_e @ vec(Ee))
overlap_e /= norm_factor * dx_forward
return unvec(overlap_e, E[0].shape)
def solve_waveguide_mode_cylindrical(mode_number: int,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
r0: float,
wavenumber_correction: bool = True,
) -> Dict[str, complex or field_t]:
"""
Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode
of the bent waveguide with the specified mode number.
:param mode_number: Number of the mode, 0-indexed
:param omega: Angular frequency of the simulation
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types.
The first coordinate is assumed to be r, the second is y.
:param epsilon: Dielectric constant
:param r0: Radius of curvature for the simulation. This should be the minimum value of
r within the simulation domain.
:param wavenumber_correction: Whether to correct the wavenumber to
account for numerical dispersion (default True)
:return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
"""
'''
Solve for the largest-magnitude eigenvalue of the real operator
'''
dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0)
eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3)
v = eigvecs[:, -(mode_number+1)]
'''
Now solve for the eigenvector of the full operator, using the real operator's
eigenvector as an initial guess for Rayleigh quotient iteration.
'''
A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0)
eigval, v = rayleigh_quotient_iteration(A, v)
# Calculate the wave-vector (force the real part to be positive)
wavenumber = numpy.sqrt(eigval)
wavenumber *= numpy.sign(numpy.real(wavenumber))
'''
Perform correction on wavenumber to account for numerical dispersion.
See Numerical Dispersion in Taflove's FDTD book.
This correction term reduces the error in emitted power, but additional
error is introduced into the E_err and H_err terms. This effect becomes
more pronounced as the wavenumber increases.
'''
if wavenumber_correction:
wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - numpy.real(wavenumber)
shape = [d.size for d in dxes[0]]
v = numpy.hstack((v, numpy.zeros(shape[0] * shape[1])))
fields = {
'wavenumber': wavenumber,
'E': unvec(v, shape),
# 'E': unvec(e, shape),
# 'H': unvec(h, shape),
}
return fields
def compute_source_q(E: field_t,
H: field_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
mu: field_t = None,
) -> field_t:
A1f = functional.curl_h(dxes)
A2f = functional.curl_e(dxes)
J = A1f(H)
M = A2f(-E)
m2j = functional.m2j(omega, dxes, mu)
Jm = m2j(M)
Jtot = J + Jm
return Jtot, J, M
def compute_source_e(QE: field_t,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
epsilon: field_t,
mu: field_t = None,
) -> field_t:
"""
Want (AQ-QA) E = -iwJ, where Q is a mask
If E is an eigenmode, AE = 0 so just AQE = -iwJ
Really only need E in 4 cells along axis (0, 0, Emode1, Emode2), find AE (1 fdtd step), then use center 2 cells as src
"""
slices = tuple(slices)
# Trim a cell from each end of the propagation axis
slices_reduced = list(slices)
slices_reduced[axis] = slice(slices[axis].start + 1, slices[axis].stop - 1)
slices_reduced = tuple(slice(None), *slices_reduced)
# Don't actually need to mask out E here since it needs to be pre-masked (QE)
A = functional.e_full(omega, dxes, epsilon, mu)
J4 = A(QE) / (-1j * omega) #J4 is 4-cell result of -iwJ = A QE
J = numpy.zeros_like(J4)
J[slices_reduced] = J4[slices_reduced]
return J
def compute_source_wg(E: field_t,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
epsilon: field_t,
mu: field_t = None,
) -> field_t:
slices = tuple(slices)
Etgt, _slices2 = compute_overlap_ce(E=E, wavenumber=wavenumber,
dxes=dxes, axis=axis, polarity=polarity,
slices=slices)
slices4 = list(slices)
slices4[axis] = slice(slices[axis].start - 4 * polarity, slices[axis].start)
slices4 = tuple(slices4)
J = compute_source_e(QE=Etgt,
omega=omega, dxes=dxes, axis=axis,
polarity=polarity, slices=slices4,
epsilon=epsilon, mu=mu)
return J
def compute_overlap_ce(E: field_t,
wavenumber: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
) -> field_t:
slices = tuple(slices)
Ee = expand_wgmode_e(E=E, wavenumber=wavenumber,
dxes=dxes, axis=axis, polarity=polarity,
slices=slices)
start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
slices2 = list(slices)
slices2[axis] = slice(start, stop)
slices2 = tuple(slice(None), slices2)
Etgt = numpy.zeros_like(Ee)
Etgt[slices2] = Ee[slices2]
Etgt /= (Etgt.conj() * Etgt).sum()
return Etgt, slices2
def expand_wgmode_e(E: field_t,
wavenumber: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
) -> field_t:
slices = tuple(slices)
# Determine phase factors for parallel slices
a_shape = numpy.roll([1, -1, 1, 1], axis)
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)
# Expand our slice to the entire grid using the phase factors
Ee = numpy.zeros_like(E)
slices_exp = list(slices)
slices_exp[axis] = slice(E.shape[axis + 1])
slices_exp = (slice(None), *slices_exp)
slices_in = tuple(slice(None), *slices)
Ee[slices_exp] = phase_E * numpy.array(E)[slices_in]
return Ee

9
meanas/fdtd/__init__.py Normal file
View file

@ -0,0 +1,9 @@
"""
Basic FDTD functionality
"""
from .base import maxwell_e, maxwell_h
from .pml import cpml
from .energy import (poynting, poynting_divergence, energy_hstep, energy_estep,
delta_energy_h2e, delta_energy_h2e, delta_energy_j)
from .boundaries import conducting_boundary

87
meanas/fdtd/base.py Normal file
View file

@ -0,0 +1,87 @@
"""
Basic FDTD field updates
"""
from typing import List, Callable, Tuple, Dict
import numpy
from .. import dx_lists_t, field_t, field_updater
__author__ = 'Jan Petykiewicz'
def curl_h(dxes: dx_lists_t = None) -> field_updater:
"""
Curl operator for use with the H field.
:param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:return: Function for taking the discretized curl of the H-field, F(H) -> curlH
"""
if dxes:
dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij')
def dh(f, ax):
return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax]
else:
def dh(f, ax):
return f - numpy.roll(f, 1, axis=ax)
def ch_fun(h: field_t) -> field_t:
output = numpy.empty_like(h)
output[0] = dh(h[2], 1)
output[1] = dh(h[0], 2)
output[2] = dh(h[1], 0)
output[0] -= dh(h[1], 2)
output[1] -= dh(h[2], 0)
output[2] -= dh(h[0], 1)
return output
return ch_fun
def curl_e(dxes: dx_lists_t = None) -> field_updater:
"""
Curl operator for use with the E field.
:param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:return: Function for taking the discretized curl of the E-field, F(E) -> curlE
"""
if dxes is not None:
dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij')
def de(f, ax):
return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax]
else:
def de(f, ax):
return numpy.roll(f, -1, axis=ax) - f
def ce_fun(e: field_t) -> field_t:
output = numpy.empty_like(e)
output[0] = de(e[2], 1)
output[1] = de(e[0], 2)
output[2] = de(e[1], 0)
output[0] -= de(e[1], 2)
output[1] -= de(e[2], 0)
output[2] -= de(e[0], 1)
return output
return ce_fun
def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater:
curl_h_fun = curl_h(dxes)
def me_fun(e: field_t, h: field_t, epsilon: field_t):
e += dt * curl_h_fun(h) / epsilon
return e
return me_fun
def maxwell_h(dt: float, dxes: dx_lists_t = None) -> field_updater:
curl_e_fun = curl_e(dxes)
def mh_fun(e: field_t, h: field_t):
h -= dt * curl_e_fun(e)
return h
return mh_fun

68
meanas/fdtd/boundaries.py Normal file
View file

@ -0,0 +1,68 @@
"""
Boundary conditions
"""
from typing import List, Callable, Tuple, Dict
import numpy
from .. import dx_lists_t, field_t, field_updater
def conducting_boundary(direction: int,
polarity: int
) -> Tuple[field_updater, field_updater]:
dirs = [0, 1, 2]
if direction not in dirs:
raise Exception('Invalid direction: {}'.format(direction))
dirs.remove(direction)
u, v = dirs
if polarity < 0:
boundary_slice = [slice(None)] * 3
shifted1_slice = [slice(None)] * 3
boundary_slice[direction] = 0
shifted1_slice[direction] = 1
def en(e: field_t):
e[direction][boundary_slice] = 0
e[u][boundary_slice] = e[u][shifted1_slice]
e[v][boundary_slice] = e[v][shifted1_slice]
return e
def hn(h: field_t):
h[direction][boundary_slice] = h[direction][shifted1_slice]
h[u][boundary_slice] = 0
h[v][boundary_slice] = 0
return h
return en, hn
elif polarity > 0:
boundary_slice = [slice(None)] * 3
shifted1_slice = [slice(None)] * 3
shifted2_slice = [slice(None)] * 3
boundary_slice[direction] = -1
shifted1_slice[direction] = -2
shifted2_slice[direction] = -3
def ep(e: field_t):
e[direction][boundary_slice] = -e[direction][shifted2_slice]
e[direction][shifted1_slice] = 0
e[u][boundary_slice] = e[u][shifted1_slice]
e[v][boundary_slice] = e[v][shifted1_slice]
return e
def hp(h: field_t):
h[direction][boundary_slice] = h[direction][shifted1_slice]
h[u][boundary_slice] = -h[u][shifted2_slice]
h[u][shifted1_slice] = 0
h[v][boundary_slice] = -h[v][shifted2_slice]
h[v][shifted1_slice] = 0
return h
return ep, hp
else:
raise Exception('Bad polarity: {}'.format(polarity))

84
meanas/fdtd/energy.py Normal file
View file

@ -0,0 +1,84 @@
from typing import List, Callable, Tuple, Dict
import numpy
from .. import dx_lists_t, field_t, field_updater
def poynting(e, h):
s = (numpy.roll(e[1], -1, axis=0) * h[2] - numpy.roll(e[2], -1, axis=0) * h[1],
numpy.roll(e[2], -1, axis=1) * h[0] - numpy.roll(e[0], -1, axis=1) * h[2],
numpy.roll(e[0], -1, axis=2) * h[1] - numpy.roll(e[1], -1, axis=2) * h[0])
return numpy.array(s)
def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes
if dxes is None:
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
if s is None:
s = poynting(e, h)
ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) / numpy.sqrt(dxes[0][0] * dxes[1][0])[:, None, None] +
(s[1] - numpy.roll(s[1], 1, axis=1)) / numpy.sqrt(dxes[0][1] * dxes[1][1])[None, :, None] +
(s[2] - numpy.roll(s[2], 1, axis=2)) / numpy.sqrt(dxes[0][2] * dxes[1][2])[None, None, :] )
return ds
def energy_hstep(e0, h1, e2, epsilon=None, mu=None, dxes=None):
u = dxmul(e0 * e2, h1 * h1, epsilon, mu, dxes)
return u
def energy_estep(h0, e1, h2, epsilon=None, mu=None, dxes=None):
u = dxmul(e1 * e1, h0 * h2, epsilon, mu, dxes)
return u
def delta_energy_h2e(dt, e0, h1, e2, h3, epsilon=None, mu=None, dxes=None):
"""
This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2)
"""
de = e2 * (e2 - e0) / dt
dh = h1 * (h3 - h1) / dt
du = dxmul(de, dh, epsilon, mu, dxes)
return du
def delta_energy_e2h(dt, h0, e1, h2, e3, epsilon=None, mu=None, dxes=None):
"""
This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2)
"""
de = e1 * (e3 - e1) / dt
dh = h2 * (h2 - h0) / dt
du = dxmul(de, dh, epsilon, mu, dxes)
return du
def delta_energy_j(j0, e1, dxes=None):
if dxes is None:
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
du = ((j0 * e1).sum(axis=0) *
dxes[0][0][:, None, None] *
dxes[0][1][None, :, None] *
dxes[0][2][None, None, :])
return du
def dxmul(ee, hh, epsilon=None, mu=None, dxes=None):
if epsilon is None:
epsilon = 1
if mu is None:
mu = 1
if dxes is None:
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
result = ((ee * epsilon).sum(axis=0) *
dxes[0][0][:, None, None] *
dxes[0][1][None, :, None] *
dxes[0][2][None, None, :] +
(hh * mu).sum(axis=0) *
dxes[1][0][:, None, None] *
dxes[1][1][None, :, None] *
dxes[1][2][None, None, :])
return result

122
meanas/fdtd/pml.py Normal file
View file

@ -0,0 +1,122 @@
"""
PML implementations
"""
# TODO retest pmls!
from typing import List, Callable, Tuple, Dict
import numpy
from .. import dx_lists_t, field_t, field_updater
__author__ = 'Jan Petykiewicz'
def cpml(direction:int,
polarity: int,
dt: float,
epsilon: field_t,
thickness: int = 8,
ln_R_per_layer: float = -1.6,
epsilon_eff: float = 1,
mu_eff: float = 1,
m: float = 3.5,
ma: float = 1,
cfs_alpha: float = 0,
dtype: numpy.dtype = numpy.float32,
) -> Tuple[Callable, Callable, Dict[str, field_t]]:
if direction not in range(3):
raise Exception('Invalid direction: {}'.format(direction))
if polarity not in (-1, 1):
raise Exception('Invalid polarity: {}'.format(polarity))
if thickness <= 2:
raise Exception('It would be wise to have a pml with 4+ cells of thickness')
if epsilon_eff <= 0:
raise Exception('epsilon_eff must be positive')
sigma_max = -ln_R_per_layer / 2 * (m + 1)
kappa_max = numpy.sqrt(epsilon_eff * mu_eff)
alpha_max = cfs_alpha
transverse = numpy.delete(range(3), direction)
u, v = transverse
xe = numpy.arange(1, thickness+1, dtype=float)
xh = numpy.arange(1, thickness+1, dtype=float)
if polarity > 0:
xe -= 0.5
elif polarity < 0:
xh -= 0.5
xe = xe[::-1]
xh = xh[::-1]
else:
raise Exception('Bad polarity!')
expand_slice = [None] * 3
expand_slice[direction] = slice(None)
def par(x):
scaling = (x / thickness) ** m
sigma = scaling * sigma_max
kappa = 1 + scaling * (kappa_max - 1)
alpha = ((1 - x / thickness) ** ma) * alpha_max
p0 = numpy.exp(-(sigma / kappa + alpha) * dt)
p1 = sigma / (sigma + kappa * alpha) * (p0 - 1)
p2 = 1 / kappa
return p0[expand_slice], p1[expand_slice], p2[expand_slice]
p0e, p1e, p2e = par(xe)
p0h, p1h, p2h = par(xh)
region = [slice(None)] * 3
if polarity < 0:
region[direction] = slice(None, thickness)
elif polarity > 0:
region[direction] = slice(-thickness, None)
else:
raise Exception('Bad polarity!')
se = 1 if direction == 1 else -1
# TODO check if epsilon is uniform in pml region?
shape = list(epsilon[0].shape)
shape[direction] = thickness
psi_e = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)]
psi_h = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)]
fields = {
'psi_e_u': psi_e[0],
'psi_e_v': psi_e[1],
'psi_h_u': psi_h[0],
'psi_h_v': psi_h[1],
}
# Note that this is kinda slow -- would be faster to reuse dHv*p2h for the original
# H update, but then you have multiple arrays and a monolithic (field + pml) update operation
def pml_e(e: field_t, h: field_t, epsilon: field_t) -> Tuple[field_t, field_t]:
dHv = h[v][region] - numpy.roll(h[v], 1, axis=direction)[region]
dHu = h[u][region] - numpy.roll(h[u], 1, axis=direction)[region]
psi_e[0] *= p0e
psi_e[0] += p1e * dHv * p2e
psi_e[1] *= p0e
psi_e[1] += p1e * dHu * p2e
e[u][region] += se * dt / epsilon[u][region] * (psi_e[0] + (p2e - 1) * dHv)
e[v][region] -= se * dt / epsilon[v][region] * (psi_e[1] + (p2e - 1) * dHu)
return e, h
def pml_h(e: field_t, h: field_t) -> Tuple[field_t, field_t]:
dEv = (numpy.roll(e[v], -1, axis=direction)[region] - e[v][region])
dEu = (numpy.roll(e[u], -1, axis=direction)[region] - e[u][region])
psi_h[0] *= p0h
psi_h[0] += p1h * dEv * p2h
psi_h[1] *= p0h
psi_h[1] += p1h * dEu * p2h
h[u][region] -= se * dt * (psi_h[0] + (p2h - 1) * dEv)
h[v][region] += se * dt * (psi_h[1] + (p2h - 1) * dEu)
return e, h
return pml_e, pml_h, fields

292
meanas/test/test_fdtd.py Normal file
View file

@ -0,0 +1,292 @@
import numpy
import pytest
import dataclasses
from typing import List, Tuple
from numpy.testing import assert_allclose, assert_array_equal
from meanas import fdtd
prng = numpy.random.RandomState(12345)
def assert_fields_close(a, b, *args, **kwargs):
numpy.testing.assert_allclose(a, b, verbose=False, err_msg='Fields did not match:\n{}\n{}'.format(numpy.rollaxis(a, -1),
numpy.rollaxis(b, -1)), *args, **kwargs)
def assert_close(a, b, *args, **kwargs):
numpy.testing.assert_allclose(a, b, *args, **kwargs)
def test_initial_fields(sim):
# Make sure initial fields didn't change
e0 = sim.es[0]
h0 = sim.hs[0]
j0 = sim.js[0]
mask = (j0 != 0)
assert_fields_close(e0[mask], j0[mask] / sim.epsilon[mask])
assert not e0[~mask].any()
assert not h0.any()
def test_initial_energy(sim):
"""
Assumes fields start at 0 before J0 is added
"""
j0 = sim.js[0]
e0 = sim.es[0]
h0 = sim.hs[0]
h1 = sim.hs[1]
mask = (j0 != 0)
dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0)
u0 = (j0 * j0.conj() / sim.epsilon * dV).sum(axis=0)
args = {'dxes': sim.dxes,
'epsilon': sim.epsilon}
# Make sure initial energy and E dot J are correct
energy0 = fdtd.energy_estep(h0=h0, e1=e0, h2=h1, **args)
e0_dot_j0 = fdtd.delta_energy_j(j0=j0, e1=e0, dxes=sim.dxes)
assert_fields_close(energy0, u0)
assert_fields_close(e0_dot_j0, u0)
def test_energy_conservation(sim):
"""
Assumes fields start at 0 before J0 is added
"""
e0 = sim.es[0]
j0 = sim.js[0]
u = fdtd.delta_energy_j(j0=j0, e1=e0, dxes=sim.dxes).sum()
args = {'dxes': sim.dxes,
'epsilon': sim.epsilon}
for ii in range(1, 8):
u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args)
u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args)
delta_j_A = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii-1], dxes=sim.dxes)
delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii], dxes=sim.dxes)
u += delta_j_A.sum()
assert_close(u_hstep.sum(), u)
u += delta_j_B.sum()
assert_close(u_estep.sum(), u)
def test_poynting_divergence(sim):
args = {'dxes': sim.dxes,
'epsilon': sim.epsilon}
dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0)
u_eprev = None
for ii in range(1, 8):
u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args)
u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args)
delta_j_B = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii], dxes=sim.dxes)
du_half_h2e = u_estep - u_hstep - delta_j_B
div_s_h2e = sim.dt * fdtd.poynting_divergence(e=sim.es[ii], h=sim.hs[ii], dxes=sim.dxes) * dV
assert_fields_close(du_half_h2e, -div_s_h2e)
if u_eprev is None:
u_eprev = u_estep
continue
# previous half-step
delta_j_A = fdtd.delta_energy_j(j0=sim.js[ii], e1=sim.es[ii-1], dxes=sim.dxes)
du_half_e2h = u_hstep - u_eprev - delta_j_A
div_s_e2h = sim.dt * fdtd.poynting_divergence(e=sim.es[ii-1], h=sim.hs[ii], dxes=sim.dxes) * dV
assert_fields_close(du_half_e2h, -div_s_e2h)
u_eprev = u_estep
def test_poynting_planes(sim):
mask = (sim.js[0] != 0)
if mask.sum() > 1:
pytest.skip('test_poynting_planes can only test single point sources')
args = {'dxes': sim.dxes,
'epsilon': sim.epsilon}
dV = numpy.prod(numpy.meshgrid(*sim.dxes[0], indexing='ij'), axis=0)
mx = numpy.roll(mask, (-1, -1), axis=(0, 1))
my = numpy.roll(mask, -1, axis=2)
mz = numpy.roll(mask, (+1, -1), axis=(0, 3))
px = numpy.roll(mask, -1, axis=0)
py = mask.copy()
pz = numpy.roll(mask, +1, axis=0)
u_eprev = None
for ii in range(1, 8):
u_hstep = fdtd.energy_hstep(e0=sim.es[ii-1], h1=sim.hs[ii], e2=sim.es[ii], **args)
u_estep = fdtd.energy_estep(h0=sim.hs[ii], e1=sim.es[ii], h2=sim.hs[ii + 1], **args)
s_h2e = -fdtd.poynting(e=sim.es[ii], h=sim.hs[ii]) * sim.dt
s_h2e[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :]
s_h2e[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :]
s_h2e[2] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][1][None, :, None]
planes = [s_h2e[px].sum(), -s_h2e[mx].sum(),
s_h2e[py].sum(), -s_h2e[my].sum(),
s_h2e[pz].sum(), -s_h2e[mz].sum()]
assert_close(sum(planes), (u_estep - u_hstep).sum())
if u_eprev is None:
u_eprev = u_estep
continue
s_e2h = -fdtd.poynting(e=sim.es[ii - 1], h=sim.hs[ii]) * sim.dt
s_e2h[0] *= sim.dxes[0][1][None, :, None] * sim.dxes[0][2][None, None, :]
s_e2h[1] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][2][None, None, :]
s_e2h[2] *= sim.dxes[0][0][:, None, None] * sim.dxes[0][1][None, :, None]
planes = [s_e2h[px].sum(), -s_e2h[mx].sum(),
s_e2h[py].sum(), -s_e2h[my].sum(),
s_e2h[pz].sum(), -s_e2h[mz].sum()]
assert_close(sum(planes), (u_hstep - u_eprev).sum())
# previous half-step
u_eprev = u_estep
#####################################
# Test fixtures
#####################################
@pytest.fixture(scope='module',
params=[(5, 5, 1),
(5, 1, 5),
(5, 5, 5),
# (7, 7, 7),
])
def shape(request):
yield (3, *request.param)
@pytest.fixture(scope='module', params=[0.3])
def dt(request):
yield request.param
@pytest.fixture(scope='module', params=[1.0, 1.5])
def epsilon_bg(request):
yield request.param
@pytest.fixture(scope='module', params=[1.0, 2.5])
def epsilon_fg(request):
yield request.param
@pytest.fixture(scope='module', params=['center', '000', 'random'])
def epsilon(request, shape, epsilon_bg, epsilon_fg):
is3d = (numpy.array(shape) == 1).sum() == 0
if is3d:
if request.param == '000':
pytest.skip('Skipping 000 epsilon because test is 3D (for speed)')
if epsilon_bg != 1:
pytest.skip('Skipping epsilon_bg != 1 because test is 3D (for speed)')
if epsilon_fg not in (1.0, 2.0):
pytest.skip('Skipping epsilon_fg not in (1, 2) because test is 3D (for speed)')
epsilon = numpy.full(shape, epsilon_bg, dtype=float)
if request.param == 'center':
epsilon[:, shape[1]//2, shape[2]//2, shape[3]//2] = epsilon_fg
elif request.param == '000':
epsilon[:, 0, 0, 0] = epsilon_fg
elif request.param == 'random':
epsilon[:] = prng.uniform(low=min(epsilon_bg, epsilon_fg),
high=max(epsilon_bg, epsilon_fg),
size=shape)
yield epsilon
@pytest.fixture(scope='module', params=[1.0])#, 1.5])
def j_mag(request):
yield request.param
@pytest.fixture(scope='module', params=['center', 'random'])
def j_distribution(request, shape, j_mag):
j = numpy.zeros(shape)
if request.param == 'center':
j[:, shape[1]//2, shape[2]//2, shape[3]//2] = j_mag
elif request.param == '000':
j[:, 0, 0, 0] = j_mag
elif request.param == 'random':
j[:] = prng.uniform(low=-j_mag, high=j_mag, size=shape)
yield j
@pytest.fixture(scope='module', params=[1.0, 1.5])
def dx(request):
yield request.param
@pytest.fixture(scope='module', params=['uniform'])
def dxes(request, shape, dx):
if request.param == 'uniform':
dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)]
yield dxes
@pytest.fixture(scope='module',
params=[(0,),
(0, 4, 8),
]
)
def j_steps(request):
yield request.param
@dataclasses.dataclass()
class SimResult:
shape: Tuple[int]
dt: float
dxes: List[List[numpy.ndarray]]
epsilon: numpy.ndarray
j_distribution: numpy.ndarray
j_steps: Tuple[int]
es: List[numpy.ndarray] = dataclasses.field(default_factory=list)
hs: List[numpy.ndarray] = dataclasses.field(default_factory=list)
js: List[numpy.ndarray] = dataclasses.field(default_factory=list)
@pytest.fixture(scope='module')
def sim(request, shape, epsilon, dxes, dt, j_distribution, j_steps):
is3d = (numpy.array(shape) == 1).sum() == 0
if is3d:
if dt != 0.3:
pytest.skip('Skipping dt != 0.3 because test is 3D (for speed)')
sim = SimResult(
shape=shape,
dt=dt,
dxes=dxes,
epsilon=epsilon,
j_distribution=j_distribution,
j_steps=j_steps,
)
e = numpy.zeros_like(epsilon)
h = numpy.zeros_like(epsilon)
assert 0 in j_steps
j_zeros = numpy.zeros_like(j_distribution)
eh2h = fdtd.maxwell_h(dt=dt, dxes=dxes)
eh2e = fdtd.maxwell_e(dt=dt, dxes=dxes)
for tt in range(10):
e = e.copy()
h = h.copy()
eh2h(e, h)
eh2e(e, h, epsilon)
if tt in j_steps:
e += j_distribution / epsilon
sim.js.append(j_distribution)
else:
sim.js.append(j_zeros)
sim.es.append(e)
sim.hs.append(h)
return sim

22
meanas/types.py Normal file
View file

@ -0,0 +1,22 @@
"""
Types shared across multiple submodules
"""
import numpy
from typing import List, Callable
# Field types
field_t = numpy.ndarray # vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])
vfield_t = numpy.ndarray # linearized vector field (vector of length 3*X*Y*Z)
'''
'dxes' datastructure which contains grid cell width information in the following format:
[[[dx_e_0, dx_e_1, ...], [dy_e_0, ...], [dz_e_0, ...]],
[[dx_h_0, dx_h_1, ...], [dy_h_0, ...], [dz_h_0, ...]]]
where dx_e_0 is the x-width of the x=0 cells, as used when calculating dE/dx,
and dy_h_0 is the y-width of the y=0 cells, as used when calculating dH/dy, etc.
'''
dx_lists_t = List[List[numpy.ndarray]]
field_updater = Callable[[field_t], field_t]

47
meanas/vectorization.py Normal file
View file

@ -0,0 +1,47 @@
"""
Functions for moving between a vector field (list of 3 ndarrays, [f_x, f_y, f_z])
and a 1D array representation of that field [f_x0, f_x1, f_x2,... f_y0,... f_z0,...].
Vectorized versions of the field use row-major (ie., C-style) ordering.
"""
from typing import List
import numpy
from .types import field_t, vfield_t
__author__ = 'Jan Petykiewicz'
def vec(f: field_t) -> vfield_t:
"""
Create a 1D ndarray from a 3D vector field which spans a 1-3D region.
Returns None if called with f=None.
:param f: A vector field, [f_x, f_y, f_z] where each f_ component is a 1 to
3D ndarray (f_* should all be the same size). Doesn't fail with f=None.
:return: A 1D ndarray containing the linearized field (or None)
"""
if numpy.any(numpy.equal(f, None)):
return None
return numpy.ravel(f, order='C')
def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t:
"""
Perform the inverse of vec(): take a 1D ndarray and output a 3D field
of form [f_x, f_y, f_z] where each of f_* is a len(shape)-dimensional
ndarray.
Returns None if called with v=None.
:param v: 1D ndarray representing a 3D vector field of shape shape (or None)
:param shape: shape of the vector field
:return: [f_x, f_y, f_z] where each f_ is a len(shape) dimensional ndarray
(or None)
"""
if numpy.any(numpy.equal(v, None)):
return None
return v.reshape((3, *shape), order='C')