398 lines
13 KiB
Python
398 lines
13 KiB
Python
|
"""
|
||
|
|
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
|
||
|
|
fdfd_tools.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 fdfd_tools.dx_lists_type,
|
||
|
|
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.
|
||
|
|
|
||
|
|
|
||
|
|
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
|
||
|
|
) -> 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 fdfd_tools.operators header
|
||
|
|
:param epsilon: Vectorized dielectric constant
|
||
|
|
:param mu: Vectorized magnetic permeability (default 1 everywhere)..
|
||
|
|
:return: Sparse matrix containing the wave operator
|
||
|
|
"""
|
||
|
|
ce = curl_e(dxes)
|
||
|
|
ch = curl_h(dxes)
|
||
|
|
|
||
|
|
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 = ch @ m_div @ ce - omega**2 * e
|
||
|
|
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 fdfd_tools.operators header
|
||
|
|
: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
|
||
|
|
) -> 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 fdfd_tools.operators header
|
||
|
|
:param epsilon: Vectorized dielectric constant
|
||
|
|
:param mu: Vectorized magnetic permeability (default 1 everywhere)
|
||
|
|
:return: Sparse matrix containing the wave operator
|
||
|
|
"""
|
||
|
|
ec = curl_e(dxes)
|
||
|
|
hc = curl_h(dxes)
|
||
|
|
|
||
|
|
e_div = sparse.diags(1 / epsilon)
|
||
|
|
if numpy.any(numpy.equal(mu, None)):
|
||
|
|
m = sparse.eye(epsilon.size)
|
||
|
|
else:
|
||
|
|
m = sparse.diags(mu)
|
||
|
|
|
||
|
|
A = ec @ e_div @ hc - omega**2 * m
|
||
|
|
return A
|
||
|
|
|
||
|
|
|
||
|
|
def eh_full(omega, dxes, epsilon, mu=None):
|
||
|
|
"""
|
||
|
|
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 fdfd_tools.operators header
|
||
|
|
:param epsilon: Vectorized dielectric constant
|
||
|
|
:param mu: Vectorized magnetic permeability (default 1 everywhere)
|
||
|
|
:return: Sparse matrix containing the wave operator
|
||
|
|
"""
|
||
|
|
A2 = curl_e(dxes)
|
||
|
|
A1 = curl_h(dxes)
|
||
|
|
|
||
|
|
iwe = 1j * omega * sparse.diags(epsilon)
|
||
|
|
iwm = 1j * omega
|
||
|
|
if not numpy.any(numpy.equal(mu, None)):
|
||
|
|
iwm *= sparse.diags(mu)
|
||
|
|
|
||
|
|
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 fdfd_tools.operators header
|
||
|
|
: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 fdfd_tools.operators header
|
||
|
|
: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,
|
||
|
|
) -> 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 fdfd_tools.operators header
|
||
|
|
:param mu: Vectorized magnetic permeability (default 1 everywhere)
|
||
|
|
: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
|
||
|
|
|
||
|
|
return op
|
||
|
|
|
||
|
|
|
||
|
|
def m2j(omega: complex,
|
||
|
|
dxes: dx_lists_t,
|
||
|
|
mu: vfield_t = None):
|
||
|
|
"""
|
||
|
|
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 fdfd_tools.operators header
|
||
|
|
: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]) -> sparse.spmatrix:
|
||
|
|
"""
|
||
|
|
Utility operator for performing a circular shift along a specified axis by 1 element.
|
||
|
|
|
||
|
|
:param axis: Axis to shift along. x=0, y=1, z=2
|
||
|
|
:param shape: Shape of the grid being shifted
|
||
|
|
: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))
|
||
|
|
|
||
|
|
n = numpy.prod(shape)
|
||
|
|
|
||
|
|
shifts = [1 if k == axis else 0 for k in range(3)]
|
||
|
|
shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)]
|
||
|
|
ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
|
||
|
|
|
||
|
|
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.flatten(order='F')))
|
||
|
|
|
||
|
|
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) - sparse.eye(n)
|
||
|
|
|
||
|
|
Ds = [sparse.diags(+1 / dx.flatten(order='F')) @ 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) - sparse.eye(n)
|
||
|
|
|
||
|
|
Ds = [sparse.diags(-1 / dx.flatten(order='F')) @ deriv(a).T
|
||
|
|
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, contining 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 fdfd_tools.operators header
|
||
|
|
: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.flatten(order='F') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
|
||
|
|
dbgx, dbgy, dbgz = [sparse.diags(dx.flatten(order='F'))
|
||
|
|
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 fdfd_tools.operators header
|
||
|
|
: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.flatten(order='F') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
|
||
|
|
dagx, dagy, dagz = [sparse.diags(dx.flatten(order='F'))
|
||
|
|
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
|