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

Compare commits

base: jan:ongoing
Branches Tags
jan:master
jan:wip
jan:ongoing
jan:fdtd_extras
jan:own_cg
jan:blochsolver
jan:release
jan:fdtd
..
compare: jan:master
Branches Tags
jan:wip
jan:master
jan:ongoing
jan:fdtd_extras
jan:own_cg
jan:blochsolver
jan:release
jan:fdtd

1 Commits
ongoing ... master

Author SHA1 Message Date
Jan Petykiewicz
03fc9e6d70 deprecate 2020-02-19 19:10:18 -08:00
30 changed files with 770 additions and 1743 deletions
Show Stats Expand all files Collapse all files

3
MANIFEST.in
View File

@ -1,3 +0,0 @@
include README.md
include LICENSE.md
include meanas/VERSION

52
README.md
View File

@ -1,56 +1,52 @@
# meanas
# fdfd_tools
**meanas** is a python package for electromagnetic simulations
** DEPRECATED **
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.
The functionality in this module is now provided by [meanas](https://mpxd.net/code/jan/meanas).
-----------------------
**fdfd_tools** is a python package containing utilities for
creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD)
electromagnetic simulations.
**Contents**
- Finite difference frequency domain (FDFD)
* Library of sparse matrices for representing the electromagnetic wave
* 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
* Waveguide mode solver and waveguide mode operators
* Stretched-coordinate PML boundaries (SCPML)
* Functional versions of most operators
* Anisotropic media (eps_xx, eps_yy, eps_zz, mu_xx, ...)
* Arbitrary distributions of perfect electric and magnetic conductors (PEC / PMC)
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
```fdfd_tools.solvers.generic(...)``` will call
```scipy.sparse.linalg.qmr(...)``` to perform a solve.
For 2D 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
For solving large (or 3D) problems, I recommend a GPU-based iterative
solver, such as [opencl_fdfd](https://mpxd.net/gogs/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.
## Installation
**Requirements:**
* python 3 (tests require 3.7)
* python 3 (written and tested with 3.5)
* numpy
* scipy
Install with pip, via git:
```bash
pip install git+https://mpxd.net/code/jan/meanas.git@release
pip install git+https://mpxd.net/gogs/jan/fdfd_tools.git@release
```
## Use
See `examples/` for some simple examples; you may need additional
packages such as [gridlock](https://mpxd.net/code/jan/gridlock)
See examples/test.py for some simple examples; you may need additional
packages such as [gridlock](https://mpxd.net/gogs/jan/gridlock)
to run the examples.

51
examples/bloch.py
View File

@ -1,44 +1,12 @@
import numpy, scipy, gridlock, meanas
from meanas.fdfd import bloch
import numpy, scipy, gridlock, fdfd_tools
from fdfd_tools import bloch
from numpy.linalg import norm
import logging
from pathlib import Path
logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger(__name__)
WISDOM_FILEPATH = pathlib.Path.home() / '.local' / 'share' / 'pyfftw' / 'wisdom.pickle'
def pyfftw_save_wisdom(path):
path = pathlib.Path(path)
try:
import pyfftw
import pickle
except ImportError as e:
pass
path.parent.mkdir(parents=True, exist_ok=True)
with open(path, 'wb') as f:
pickle.dump(wisdom, f)
def pyfftw_load_wisdom(path):
path = pathlib.Path(path)
try:
import pyfftw
import pickle
except ImportError as e:
pass
try:
with open(path, 'rb') as f:
wisdom = pickle.load(f)
pyfftw.import_wisdom(wisdom)
except FileNotFoundError as e:
pass
logger.info('Drawing grid...')
dx = 40
x_period = 400
y_period = z_period = 2000
@ -62,13 +30,11 @@ g2.shifts = numpy.zeros((6,3))
g2.grids = [numpy.zeros(g.shape) for _ in range(6)]
epsilon = [g.grids[0],] * 3
reciprocal_lattice = numpy.diag(1000/numpy.array([x_period, y_period, z_period])) #cols are vectors
pyfftw_load_wisdom(WISDOM_FILEPATH)
reciprocal_lattice = numpy.diag(1e6/numpy.array([x_period, y_period, z_period])) #cols are vectors
#print('Finding k at 1550nm')
#k, f = bloch.find_k(frequency=1000/1550,
# tolerance=(1000 * (1/1550 - 1/1551)),
#k, f = bloch.find_k(frequency=1/1550,
# tolerance=(1/1550 - 1/1551),
# direction=[1, 0, 0],
# G_matrix=reciprocal_lattice,
# epsilon=epsilon,
@ -76,15 +42,15 @@ pyfftw_load_wisdom(WISDOM_FILEPATH)
#
#print("k={}, f={}, 1/f={}, k/f={}".format(k, f, 1/f, norm(reciprocal_lattice @ k) / f ))
logger.info('Finding f at [0.25, 0, 0]')
print('Finding f at [0.25, 0, 0]')
for k0x in [.25]:
k0 = numpy.array([k0x, 0, 0])
kmag = norm(reciprocal_lattice @ k0)
tolerance = (1000/1550) * 1e-4/1.5 # df = f * dn_eff / n
tolerance = (1e6/1550) * 1e-4/1.5 # df = f * dn_eff / n
logger.info('tolerance {}'.format(tolerance))
n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance**2)
n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance)
v2e = bloch.hmn_2_exyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
v2h = bloch.hmn_2_hxyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
ki = bloch.generate_kmn(k0, reciprocal_lattice, g.shape)
@ -100,4 +66,3 @@ for k0x in [.25]:
n_eff = norm(reciprocal_lattice @ k0) / f
print('kmag/f = n_eff = {} \n wl = {}\n'.format(n_eff, 1/f ))
pyfftw_save_wisdom(WISDOM_FILEPATH)

90
examples/tcyl.py
View File

@ -1,90 +0,0 @@
import importlib
import numpy
from numpy.linalg import norm
from meanas import vec, unvec
from meanas.fdfd import waveguide_mode, functional, scpml
from meanas.fdfd.solvers import generic as generic_solver
import gridlock
from matplotlib import pyplot
__author__ = 'Jan Petykiewicz'
def test1(solver=generic_solver):
dx = 20 # discretization (nm/cell)
pml_thickness = 10 # (number of cells)
wl = 1550 # Excitation wavelength
omega = 2 * numpy.pi / wl
# Device design parameters
w = 800
th = 220
center = [0, 0, 0]
r0 = 8e3
# refractive indices
n_wg = numpy.sqrt(12.6) # ~Si
n_air = 1.0 # air
# Half-dimensions of the simulation grid
y_max = 1200
z_max = 900
xyz_max = numpy.array([800, y_max, z_max]) + (pml_thickness + 2) * dx
# Coordinates of the edges of the cells.
half_edge_coords = [numpy.arange(dx/2, m + dx/2, step=dx) for m in xyz_max]
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
edge_coords[0] = numpy.array([-dx, dx])
# #### Create the grid and draw the device ####
grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3)
grid.draw_cuboid(center=center, dimensions=[8e3, w, th], eps=n_wg**2)
dxes = [grid.dxyz, grid.autoshifted_dxyz()]
for a in (1, 2):
for p in (-1, 1):
dxes = scmpl.stretch_with_scpml(dxes, omega=omega, axis=a, polarity=p,
thickness=pml_thickness)
wg_args = {
'omega': omega,
'dxes': [(d[1], d[2]) for d in dxes],
'epsilon': vec(g.transpose([1, 2, 0]) for g in grid.grids),
'r0': r0,
}
wg_results = waveguide_mode.solve_waveguide_mode_cylindrical(mode_number=0, **wg_args)
E = wg_results['E']
n_eff = wl / (2 * numpy.pi / wg_results['wavenumber'])
print('n =', n_eff)
print('alpha (um^-1) =', -4 * numpy.pi * numpy.imag(n_eff) / (wl * 1e-3))
'''
Plot results
'''
def pcolor(v):
vmax = numpy.max(numpy.abs(v))
pyplot.pcolor(v.T, cmap='seismic', vmin=-vmax, vmax=vmax)
pyplot.axis('equal')
pyplot.colorbar()
pyplot.figure()
pyplot.subplot(2, 2, 1)
pcolor(numpy.real(E[0][:, :]))
pyplot.subplot(2, 2, 2)
pcolor(numpy.real(E[1][:, :]))
pyplot.subplot(2, 2, 3)
pcolor(numpy.real(E[2][:, :]))
pyplot.subplot(2, 2, 4)
pyplot.show()
if __name__ == '__main__':
test1()

87
examples/fdfd.py → examples/test.py
View File

@ -2,10 +2,11 @@ import importlib
import numpy
from numpy.linalg import norm
import meanas
from meanas import vec, unvec, fdtd
from meanas.fdfd import waveguide_mode, functional, scpml, operators
from meanas.fdfd.solvers import generic as generic_solver
from fdfd_tools import vec, unvec, waveguide_mode
import fdfd_tools
import fdfd_tools.functional
import fdfd_tools.grid
from fdfd_tools.solvers import generic as generic_solver
import gridlock
@ -56,24 +57,18 @@ def test0(solver=generic_solver):
dxes = [grid.dxyz, grid.autoshifted_dxyz()]
for a in (0, 1, 2):
for p in (-1, 1):
dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega,
dxes = fdfd_tools.grid.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega,
thickness=pml_thickness)
J = [numpy.zeros_like(grid.grids[0], dtype=complex) for _ in range(3)]
J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1
J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1e5
'''
Solve!
'''
sim_args = {
'omega': omega,
'dxes': dxes,
'epsilon': vec(grid.grids),
}
x = solver(J=vec(J), **sim_args)
A = operators.e_full(omega, dxes, vec(grid.grids)).tocsr()
A = fdfd_tools.functional.e_full(omega, dxes, vec(grid.grids)).tocsr()
b = -1j * omega * vec(J)
print('Norm of the residual is ', norm(A @ x - b))
@ -118,7 +113,7 @@ def test1(solver=generic_solver):
dxes = [grid.dxyz, grid.autoshifted_dxyz()]
for a in (0, 1, 2):
for p in (-1, 1):
dxes = scpml.stretch_with_scpml(dxes,omega=omega, axis=a, polarity=p,
dxes = fdfd_tools.grid.stretch_with_scpml(dxes,omega=omega, axis=a, polarity=p,
thickness=pml_thickness)
half_dims = numpy.array([10, 20, 15]) * dx
@ -126,18 +121,17 @@ def test1(solver=generic_solver):
dims[1][0] = dims[0][0]
ind_dims = (grid.pos2ind(dims[0], which_shifts=None).astype(int),
grid.pos2ind(dims[1], which_shifts=None).astype(int))
src_axis = 0
wg_args = {
'omega': omega,
'slices': [slice(i, f+1) for i, f in zip(*ind_dims)],
'dxes': dxes,
'axis': src_axis,
'axis': 0,
'polarity': +1,
}
wg_results = waveguide_mode.solve_waveguide_mode(mode_number=0, omega=omega, epsilon=grid.grids, **wg_args)
J = waveguide_mode.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'],
omega=omega, epsilon=grid.grids, **wg_args)
e_overlap = waveguide_mode.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args)
wg_results = waveguide_mode.solve_waveguide_mode(mode_number=0, **wg_args, epsilon=grid.grids)
J = waveguide_mode.compute_source(**wg_args, **wg_results)
H_overlap = waveguide_mode.compute_overlap_e(**wg_args, **wg_results)
pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3)
# pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
@ -147,19 +141,6 @@ def test1(solver=generic_solver):
# pmcg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
# pmcg.visualize_isosurface()
def pcolor(v):
vmax = numpy.max(numpy.abs(v))
pyplot.pcolor(v, cmap='seismic', vmin=-vmax, vmax=vmax)
pyplot.axis('equal')
pyplot.colorbar()
ss = (1, slice(None), J.shape[2]//2+6, slice(None))
# pyplot.figure()
# pcolor(J3[ss].T.imag)
# pyplot.figure()
# pcolor((numpy.abs(J3).sum(axis=2).sum(axis=0) > 0).astype(float).T)
pyplot.show(block=True)
'''
Solve!
'''
@ -174,7 +155,7 @@ def test1(solver=generic_solver):
x = solver(J=vec(J), **sim_args)
b = -1j * omega * vec(J)
A = operators.e_full(**sim_args).tocsr()
A = fdfd_tools.operators.e_full(**sim_args).tocsr()
print('Norm of the residual is ', norm(A @ x - b))
E = unvec(x, grid.shape)
@ -182,38 +163,40 @@ def test1(solver=generic_solver):
'''
Plot results
'''
def pcolor(v):
vmax = numpy.max(numpy.abs(v))
pyplot.pcolor(v, cmap='seismic', vmin=-vmax, vmax=vmax)
pyplot.axis('equal')
pyplot.colorbar()
center = grid.pos2ind([0, 0, 0], None).astype(int)
pyplot.figure()
pyplot.subplot(2, 2, 1)
pcolor(numpy.real(E[1][center[0], :, :]).T)
pcolor(numpy.real(E[1][center[0], :, :]))
pyplot.subplot(2, 2, 2)
pyplot.plot(numpy.log10(numpy.abs(E[1][:, center[1], center[2]]) + 1e-10))
pyplot.subplot(2, 2, 3)
pcolor(numpy.real(E[1][:, :, center[2]]).T)
pcolor(numpy.real(E[1][:, :, center[2]]))
pyplot.subplot(2, 2, 4)
def poyntings(E):
H = functional.e2h(omega, dxes)(E)
poynting = 0.5 * fdtd.poynting(e=E, h=H.conj()) * dx * dx
cross1 = operators.poynting_e_cross(vec(E), dxes) @ vec(H).conj()
cross2 = operators.poynting_h_cross(vec(H), dxes) @ vec(E).conj() * -1
e = vec(E)
h = fdfd_tools.operators.e2h(omega, dxes) @ e
cross1 = fdfd_tools.operators.poynting_e_cross(e, dxes) @ h.conj()
cross2 = fdfd_tools.operators.poynting_h_cross(h.conj(), dxes) @ e
s1 = unvec(0.5 * numpy.real(cross1), grid.shape)
s2 = unvec(0.5 * numpy.real(cross2), grid.shape)
s0 = poynting.real
# s2 = poynting.imag
return s0, s1, s2
s2 = unvec(0.5 * numpy.real(-cross2), grid.shape)
return s1, s2
s0x, s1x, s2x = poyntings(E)
pyplot.plot(s0x[0].sum(axis=2).sum(axis=1), label='s0')
pyplot.plot(s1x[0].sum(axis=2).sum(axis=1), label='s1')
pyplot.plot(s2x[0].sum(axis=2).sum(axis=1), label='s2')
pyplot.legend()
s1x, s2x = poyntings(E)
pyplot.plot(s1x[0].sum(axis=2).sum(axis=1))
pyplot.plot(s2x[0].sum(axis=2).sum(axis=1))
pyplot.show()
q = []
for i in range(-5, 30):
e_ovl_rolled = numpy.roll(e_overlap, i, axis=1)
q += [numpy.abs(vec(E) @ vec(e_ovl_rolled).conj())]
H_rolled = [numpy.roll(h, i, axis=0) for h in H_overlap]
q += [numpy.abs(vec(E) @ vec(H_rolled))]
pyplot.figure()
pyplot.plot(q)
pyplot.title('Overlap with mode')
@ -226,7 +209,7 @@ def module_available(name):
if __name__ == '__main__':
#test0()
# test0()
if module_available('opencl_fdfd'):
from opencl_fdfd import cg_solver as opencl_solver

2
examples/fdtd.py → examples/test_fdtd.py
View File

@ -10,7 +10,7 @@ import time
import numpy
import h5py
from meanas import fdtd
from fdfd_tools import fdtd
from masque import Pattern, shapes
import gridlock
import pcgen

25
fdfd_tools/__init__.py Normal file
View File

@ -0,0 +1,25 @@
"""
Electromagnetic FDFD simulation tools
Tools for 3D and 2D Electromagnetic Finite Difference Frequency Domain (FDFD)
simulations. These tools handle conversion of fields to/from vector form,
creation of the wave operator matrix, stretched-coordinate PMLs, PECs and PMCs,
field conversion operators, waveguide mode operator, and waveguide mode
solver.
This package only contains a solver for the waveguide mode eigenproblem;
if you want to solve 3D problems you can use your favorite iterative sparse
matrix solver (so long as it can handle complex symmetric [non-Hermitian]
matrices, ideally with double precision).
Dependencies:
- numpy
- scipy
"""
from .vectorization import vec, unvec, field_t, vfield_t
from .grid import dx_lists_t
__author__ = 'Jan Petykiewicz'

416
meanas/fdfd/bloch.py → fdfd_tools/bloch.py
View File

@ -73,46 +73,21 @@ This module contains functions for generating and solving the
'''
from typing import Tuple, Callable
from typing import List, Tuple, Callable, Dict
import logging
import numpy
from numpy import pi, real, trace
from numpy.fft import fftfreq
from numpy.fft import fftn, ifftn, fftfreq
import scipy
import scipy.optimize
from scipy.linalg import norm
import scipy.sparse.linalg as spalg
from .eigensolvers import rayleigh_quotient_iteration
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
@ -280,7 +255,7 @@ def hmn_2_hxyz(k0: numpy.ndarray,
:return: Function for converting h_mn into H_xyz
"""
shape = epsilon[0].shape + (1,)
_k_mag, m, n = generate_kmn(k0, G_matrix, shape)
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)]
@ -354,14 +329,147 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray,
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
h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag
h_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] / -k_mag
return numpy.hstack((h_m.ravel(), h_n.ravel()))
return operator
def eigsolve(num_modes: int,
k0: numpy.ndarray,
G_matrix: numpy.ndarray,
epsilon: field_t,
mu: field_t = None,
tolerance = 1e-8,
) -> 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).
: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)
def rayleigh_quotient(Z: numpy.ndarray, approx_grad: bool = True):
"""
Absolute value of the block Rayleigh quotient, and the associated gradient.
See Johnson and Joannopoulos, Opt. Expr. 8, 3 (2001) for details (full
citation in module docstring).
===
Notes on my understanding of the procedure:
Minimize f(Y) = |trace((Y.H @ A @ Y)|, making use of Y = Z @ inv(Z.H @ Z)^(1/2)
(a polar orthogonalization of Y). This gives f(Z) = |trace(Z.H @ A @ Z @ U)|,
where U = inv(Z.H @ Z). We minimize the absolute value to find the eigenvalues
with smallest magnitude.
The gradient is P @ (A @ Z @ U), where P = (1 - Z @ U @ Z.H) is a projection
onto the space orthonormal to Z. If approx_grad is True, the approximate
inverse of the maxwell operator is used to precondition the gradient.
"""
z = Z.view(dtype=complex).reshape(y_shape)
U = numpy.linalg.inv(z.conj().T @ z)
zU = z @ U
AzU = scipy_op @ zU
zTAzU = z.conj().T @ AzU
f = numpy.real(numpy.trace(zTAzU))
if approx_grad:
df_dy = scipy_iop @ (AzU - zU @ zTAzU)
else:
df_dy = (AzU - zU @ zTAzU)
df_dy_flat = df_dy.view(dtype=float).ravel()
return numpy.abs(f), numpy.sign(f) * df_dy_flat
'''
Use the conjugate gradient method and the approximate gradient calculation to
quickly find approximate eigenvectors.
'''
result = scipy.optimize.minimize(rayleigh_quotient,
numpy.random.rand(*y_shape, 2),
jac=True,
method='L-BFGS-B',
tol=1e-20,
options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30})
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, True),
result.x,
jac=True,
method='L-BFGS-B',
tol=1e-20,
options={'maxiter': 2000, 'gtol':0, 'disp':True})
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
result.x,
jac=True,
method='L-BFGS-B',
tol=1e-20,
options={'maxiter': 2000, 'gtol':0, 'disp':True})
for i in range(20):
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
result.x,
jac=True,
method='L-BFGS-B',
tol=1e-20,
options={'maxiter': 70, 'gtol':0, 'disp':True})
if result.nit == 0:
# We took 0 steps, so re-running won't help
break
z = result.x.view(dtype=complex).reshape(y_shape)
'''
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 find_k(frequency: float,
tolerance: float,
direction: numpy.ndarray,
@ -371,7 +479,6 @@ def find_k(frequency: float,
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.
@ -392,10 +499,8 @@ def find_k(frequency: float,
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)
n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon)
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),
@ -406,244 +511,3 @@ def find_k(frequency: float,
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

0
meanas/eigensolvers.py → fdfd_tools/eigensolvers.py
View File

10
meanas/fdfd/farfield.py → fdfd_tools/farfield.py
View File

@ -6,11 +6,9 @@ import numpy
from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
from numpy import pi
from .. import field_t
def near_to_farfield(E_near: field_t,
H_near: field_t,
def near_to_farfield(E_near: List[numpy.ndarray],
H_near: List[numpy.ndarray],
dx: float,
dy: float,
padded_size: List[int] = None
@ -117,8 +115,8 @@ def near_to_farfield(E_near: field_t,
def far_to_nearfield(E_far: field_t,
H_far: field_t,
def far_to_nearfield(E_far: List[numpy.ndarray],
H_far: List[numpy.ndarray],
dkx: float,
dky: float,
padded_size: List[int] = None

239
fdfd_tools/fdtd.py Normal file
View File

@ -0,0 +1,239 @@
from typing import List, Callable, Tuple, Dict
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 = None) -> functional_matrix:
"""
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:
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 = None) -> functional_matrix:
"""
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:
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 maxwell_e(dt: float, dxes: dx_lists_t = None) -> functional_matrix:
curl_h_fun = curl_h(dxes)
def me_fun(e: field_t, h: field_t, epsilon: field_t):
ch = curl_h_fun(h)
for ei, ci, epsi in zip(e, ch, epsilon):
ei += dt * ci / epsi
return e
return me_fun
def maxwell_h(dt: float, dxes: dx_lists_t = None) -> functional_matrix:
curl_e_fun = curl_e(dxes)
def mh_fun(e: field_t, h: field_t):
ce = curl_e_fun(e)
for hi, ci in zip(h, ce):
hi -= dt * ci
return h
return mh_fun
def conducting_boundary(direction: int,
polarity: int
) -> Tuple[functional_matrix, functional_matrix]:
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))
def cpml(direction:int,
polarity: int,
dt: float,
epsilon: field_t,
thickness: int = 8,
epsilon_eff: float = 1,
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')
m = (3.5, 1)
sigma_max = 0.8 * (m[0] + 1) / numpy.sqrt(epsilon_eff)
alpha_max = 0 # TODO: Decide what to do about non-zero 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):
sigma = ((x / thickness) ** m[0]) * sigma_max
alpha = ((1 - x / thickness) ** m[1]) * alpha_max
p0 = numpy.exp(-(sigma + alpha) * dt)
p1 = sigma / (sigma + alpha) * (p0 - 1)
return p0[expand_slice], p1[expand_slice]
p0e, p1e = par(xe)
p0h, p1h = 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!')
if direction == 1:
se = 1
else:
se = -1
# TODO check if epsilon is uniform?
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],
}
def pml_e(e: field_t, h: field_t, epsilon: field_t) -> Tuple[field_t, field_t]:
psi_e[0] *= p0e
psi_e[0] += p1e * (h[v][region] - numpy.roll(h[v], 1, axis=direction)[region])
psi_e[1] *= p0e
psi_e[1] += p1e * (h[u][region] - numpy.roll(h[u], 1, axis=direction)[region])
e[u][region] += se * dt * psi_e[0] / epsilon[u][region]
e[v][region] -= se * dt * psi_e[1] / epsilon[v][region]
return e, h
def pml_h(e: field_t, h: field_t) -> Tuple[field_t, field_t]:
psi_h[0] *= p0h
psi_h[0] += p1h * (numpy.roll(e[v], -1, axis=direction)[region] - e[v][region])
psi_h[1] *= p0h
psi_h[1] += p1h * (numpy.roll(e[u], -1, axis=direction)[region] - e[u][region])
h[u][region] -= se * dt * psi_h[0]
h[v][region] += se * dt * psi_h[1]
return e, h
return pml_e, pml_h, fields

103
meanas/fdfd/functional.py → fdfd_tools/functional.py
View File

@ -2,13 +2,13 @@
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)
The functions generated here expect inputs in the form E = [E_x, E_y, E_z], where each
component E_* is an ndarray of equal shape.
"""
from typing import List, Callable
import numpy
from .. import dx_lists_t, field_t
from . import dx_lists_t, field_t
__author__ = 'Jan Petykiewicz'
@ -20,7 +20,7 @@ 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
: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
"""
dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij')
@ -29,13 +29,9 @@ def curl_h(dxes: dx_lists_t) -> functional_matrix:
return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax]
def ch_fun(h: field_t) -> field_t:
e = numpy.empty_like(h)
e[0] = dh(h[2], 1)
e[0] -= dh(h[1], 2)
e[1] = dh(h[0], 2)
e[1] -= dh(h[2], 0)
e[2] = dh(h[1], 0)
e[2] -= dh(h[0], 1)
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
@ -45,7 +41,7 @@ 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
: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
"""
dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij')
@ -54,13 +50,9 @@ def curl_e(dxes: dx_lists_t) -> functional_matrix:
return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax]
def ce_fun(e: field_t) -> field_t:
h = numpy.empty_like(e)
h[0] = de(e[2], 1)
h[0] -= de(e[1], 2)
h[1] = de(e[0], 2)
h[1] -= de(e[2], 0)
h[2] = de(e[1], 0)
h[2] -= de(e[0], 1)
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
@ -77,7 +69,7 @@ def e_full(omega: complex,
(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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function implementing the wave operator A(E) -> E
@ -87,11 +79,11 @@ def e_full(omega: complex,
def op_1(e):
curls = ch(ce(e))
return curls - omega ** 2 * epsilon * e
return [c - omega ** 2 * e * x for c, e, x in zip(curls, epsilon, e)]
def op_mu(e):
curls = ch(mu * ce(e))
return curls - omega ** 2 * epsilon * 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
@ -108,7 +100,7 @@ def eh_full(omega: complex,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function implementing the wave operator A(E, H) -> (E, H)
@ -117,12 +109,12 @@ def eh_full(omega: complex,
ce = curl_e(dxes)
def op_1(e, h):
return (ch(h) - 1j * omega * epsilon * e,
ce(e) + 1j * omega * 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 (ch(h) - 1j * omega * epsilon * e,
ce(e) + 1j * omega * mu * 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
@ -139,68 +131,19 @@ def e2h(omega: complex,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:param mu: Magnetic permeability (default 1 everywhere)
:return: Function for converting E to H
"""
A2 = curl_e(dxes)
def e2h_1_1(e):
return A2(e) / (-1j * omega)
return [y / (-1j * omega) for y in A2(e)]
def e2h_mu(e):
return A2(e) / (-1j * omega * mu)
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):
J = ch(m / mu) / (-1j * omega)
return J
def m2j_1(m):
J = ch(m) / (-1j * omega)
return J
if numpy.any(numpy.equal(mu, None)):
return m2j_1
else:
return m2j_mu
def e_tfsf_source(TF_region: field_t,
omega: complex,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None,
) -> functional_matrix:
"""
Operator that turuns an E-field distribution into a total-field/scattered-field
(TFSF) source.
"""
# TODO documentation
A = e_full(omega, dxes, epsilon, mu)
def op(e):
neg_iwj = A(TF_region * e) - TF_region * A(e)
return neg_iwj / (-1j * omega)

3
meanas/fdfd/scpml.py → fdfd_tools/grid.py
View File

@ -5,11 +5,10 @@ Functions for creating stretched coordinate PMLs.
from typing import List, Callable
import numpy
from .. import dx_lists_t
__author__ = 'Jan Petykiewicz'
dx_lists_t = List[List[numpy.ndarray]]
s_function_type = Callable[[float], float]

134
meanas/fdfd/operators.py → fdfd_tools/operators.py
View File

@ -3,13 +3,17 @@ 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).
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 meanas.dx_lists_type; see
the meanas.types submodule for details.
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:
@ -32,7 +36,7 @@ from typing import List, Tuple
import numpy
import scipy.sparse as sparse
from .. import vec, dx_lists_t, vfield_t
from . import vec, dx_lists_t, vfield_t
__author__ = 'Jan Petykiewicz'
@ -53,7 +57,7 @@ def e_full(omega: complex,
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 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).
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
@ -97,7 +101,7 @@ def e_full_preconditioners(dxes: dx_lists_t
The preconditioner matrices are diagonal and complex, with Pr = 1 / Pl
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
: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, :],
@ -123,7 +127,7 @@ def h_full(omega: complex,
(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 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)
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
@ -173,7 +177,7 @@ def eh_full(omega: complex,
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 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)
:param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
@ -212,7 +216,7 @@ 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
: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]))
@ -222,7 +226,7 @@ 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
: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]))
@ -238,7 +242,7 @@ def e2h(omega: complex,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
: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).
@ -266,7 +270,7 @@ def m2j(omega: complex,
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 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
"""
@ -341,7 +345,9 @@ def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> spa
n = numpy.prod(shape)
i_ind = numpy.arange(n)
j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
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')))
@ -445,26 +451,30 @@ def avgb(axis: int, shape: List[int]) -> sparse.spmatrix:
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.
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 meanas.types
: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 = [rotation(i, shape, 1) for i in range(3)]
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')]
dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
Ex, Ey, Ez = [ei * da for ei, da in zip(numpy.split(e, 3), dxag)]
dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C'))
for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
block_diags = [[ None, fx @ -Ez, fx @ Ey],
[ fy @ Ez, None, fy @ -Ex],
[ fz @ -Ey, fz @ Ex, None]]
block_matrix = sparse.bmat([[sparse.diags(x) if x is not None else None for x in row]
for row in block_diags])
P = block_matrix @ sparse.diags(numpy.concatenate(dxag))
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
@ -473,75 +483,25 @@ 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
: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 = [rotation(i, shape, 1) for i in range(3)]
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')]
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)]
P = (sparse.bmat(
[[ None, -Hz @ fx, Hy @ fx],
[ Hz @ fy, None, -Hx @ fy],
[-Hy @ fz, Hx @ fz, None]])
@ sparse.diags(numpy.concatenate(dxag)))
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
def e_tfsf_source(TF_region: vfield_t,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
) -> sparse.spmatrix:
"""
Operator that turns an E-field distribution into a total-field/scattered-field
(TFSF) source.
"""
# TODO documentation
A = e_full(omega, dxes, epsilon, mu)
Q = sparse.diags(TF_region)
return (A @ Q - Q @ A) / (-1j * omega)
def e_boundary_source(mask: vfield_t,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
periodic_mask_edges: bool = False,
) -> sparse.spmatrix:
"""
Operator that turns an E-field distrubtion into a current (J) distribution
along the edges (external and internal) of the provided mask. This is just an
e_tfsf_source with an additional masking step.
"""
full = e_tfsf_source(TF_region=mask, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu)
shape = [len(dxe) for dxe in dxes[0]]
jmask = numpy.zeros_like(mask, dtype=bool)
if periodic_mask_edges:
shift = lambda axis, polarity: rotation(axis=axis, shape=shape, shift_distance=polarity)
else:
shift = lambda axis, polarity: shift_with_mirror(axis=axis, shape=shape, shift_distance=polarity)
for axis in (0, 1, 2):
for polarity in (-1, +1):
r = shift(axis, polarity) - sparse.eye(numpy.prod(shape)) # shifted minus original
r3 = sparse.block_diag((r, r, r))
jmask = numpy.logical_or(jmask, numpy.abs(r3 @ mask))
# jmask = ((numpy.roll(mask, -1, axis=0) != mask) |
# (numpy.roll(mask, +1, axis=0) != mask) |
# (numpy.roll(mask, -1, axis=1) != mask) |
# (numpy.roll(mask, +1, axis=1) != mask) |
# (numpy.roll(mask, -1, axis=2) != mask) |
# (numpy.roll(mask, +1, axis=2) != mask))
return sparse.diags(jmask.astype(int)) @ full

2
meanas/fdfd/solvers.py → fdfd_tools/solvers.py
View File

@ -70,7 +70,7 @@ def generic(omega: complex,
"""
Conjugate gradient FDFD solver using CSR sparse matrices.
All ndarray arguments should be 1D array, as returned by meanas.vec().
All ndarray arguments should be 1D array, as returned by fdfd_tools.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)

12
meanas/vectorization.py → fdfd_tools/vectorization.py
View File

@ -4,14 +4,16 @@ 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'
# Types
field_t = List[numpy.ndarray] # vector field (eg. [E_x, E_y, E_z]
vfield_t = numpy.ndarray # linearized vector field
def vec(f: field_t) -> vfield_t:
"""
@ -25,7 +27,7 @@ def vec(f: field_t) -> vfield_t:
"""
if numpy.any(numpy.equal(f, None)):
return None
return numpy.ravel(f, order='C')
return numpy.hstack(tuple((fi.ravel(order='C') for fi in f)))
def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t:
@ -43,5 +45,5 @@ def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t:
"""
if numpy.any(numpy.equal(v, None)):
return None
return v.reshape((3, *shape), order='C')
return [vi.reshape(shape, order='C') for vi in numpy.split(v, 3)]

232
meanas/fdfd/waveguide.py → fdfd_tools/waveguide.py
View File

@ -17,46 +17,20 @@ 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.
"""
# TODO update module docs
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 vec, unvec, dx_lists_t, field_t, vfield_t
from . import operators
__author__ = 'Jan Petykiewicz'
def operator_e(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
) -> sparse.spmatrix:
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_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1])))
eps_z_inv = sparse.diags(1 / eps_parts[2])
mu_parts = numpy.split(mu, 3)
mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0])))
mu_z_inv = sparse.diags(1 / mu_parts[2])
op = omega * omega * mu_yx @ eps_xy + \
mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \
sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy
return op
def operator_h(omega: complex,
def operator(omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
@ -77,7 +51,7 @@ def operator_h(omega: complex,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Sparse matrix representation of the operator
@ -96,101 +70,51 @@ def operator_h(omega: complex,
mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
mu_z_inv = sparse.diags(1 / mu_parts[2])
op = omega * omega * eps_yx @ mu_xy + \
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_e(e_xy: numpy.ndarray,
def normalized_fields(v: numpy.ndarray,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
prop_phase: float = 0,
mu: vfield_t = None
) -> Tuple[vfield_t, vfield_t]:
"""
Given a vector e_xy containing the vectorized E_x and E_y fields,
Given a vector v containing the vectorized H_x and H_y fields,
returns normalized, vectorized E and H fields for the system.
:param e_xy: Vector containing E_x and E_y fields
:param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:param prop_phase: Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction.
Default 0 (continuous propagation direction, i.e. dz->0).
:return: Normalized, vectorized (e, h) containing all vector components.
"""
e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy
h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy
e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
mu=mu, prop_phase=prop_phase)
return e_norm, h_norm
e = v2e(v, wavenumber, omega, dxes, epsilon, mu=mu)
h = v2h(v, wavenumber, dxes, mu=mu)
def normalized_fields_h(h_xy: numpy.ndarray,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
prop_phase: float = 0,
) -> Tuple[vfield_t, vfield_t]:
"""
Given a vector e_xy containing the vectorized E_x and E_y fields,
returns normalized, vectorized E and H fields for the system.
:param e_xy: Vector containing E_x and E_y fields
:param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
: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)
:param dxes_prop: Grid cell width in the propagation direction. Default 0 (continuous).
:return: Normalized, vectorized (e, h) containing all vector components.
"""
e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy
h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy
e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
mu=mu, prop_phase=prop_phase)
return e_norm, h_norm
def _normalized_fields(e: numpy.ndarray,
h: numpy.ndarray,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
prop_phase: float = 0,
) -> Tuple[vfield_t, vfield_t]:
# TODO documentation
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)
phase = numpy.exp(-1j * prop_phase / 2)
S1 = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0]
S2 = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1]
P = numpy.real(S1.sum() - S2.sum())
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 TODO]
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())
norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
norm_angle = -numpy.angle(e[e.size//2])
norm_factor = norm_amplitude * numpy.exp(1j * norm_angle)
e *= norm_factor
h *= norm_factor
@ -198,104 +122,56 @@ def _normalized_fields(e: numpy.ndarray,
return e, h
def exy2h(wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> sparse.spmatrix:
"""
Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields,
into a vectorized H containing all three H components
:param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
: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 representing the operator
"""
e2hop = e2h(wavenumber=wavenumber, omega=omega, dxes=dxes, mu=mu)
return e2hop @ exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon)
def hxy2e(wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None
) -> sparse.spmatrix:
"""
Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields,
into a vectorized E containing all three E components
:param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy`
: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 representing the operator
"""
h2eop = h2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon)
return h2eop @ hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu)
def hxy2h(wavenumber: complex,
def v2h(v: numpy.ndarray,
wavenumber: complex,
dxes: dx_lists_t,
mu: vfield_t = None
) -> sparse.spmatrix:
) -> vfield_t:
"""
Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields,
into a vectorized H containing all three H components
Given a vector v containing the vectorized H_x and H_y fields,
returns a vectorized H including all three H components.
:param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy`
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
: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 fdfd_tools.operators header (2D)
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Sparse matrix representing the operator
:return: Vectorized H field with all vector components
"""
Dfx, Dfy = operators.deriv_forward(dxes[0])
hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber)
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])
hxy2hz = mu_z_inv @ hxy2hz @ mu_xy
op = mu_z_inv @ op @ mu_xy
n_pts = dxes[1][0].size * dxes[1][1].size
op = sparse.vstack((sparse.eye(2 * n_pts),
hxy2hz))
return op
w = op @ v / (1j * wavenumber)
return numpy.hstack((v, w)).flatten()
def exy2e(wavenumber: complex,
def v2e(v: numpy.ndarray,
wavenumber: complex,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
) -> sparse.spmatrix:
mu: vfield_t = None
) -> vfield_t:
"""
Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields,
into a vectorized E containing all three E components
Given a vector v containing the vectorized H_x and H_y fields,
returns a vectorized E containing all three E components
:param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
:param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
: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 fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:return: Sparse matrix representing the operator
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Vectorized E field with all vector components.
"""
Dbx, Dby = operators.deriv_back(dxes[1])
exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber)
if not numpy.any(numpy.equal(epsilon, None)):
epsilon_parts = numpy.split(epsilon, 3)
epsilon_xy = sparse.diags(numpy.hstack((epsilon_parts[0], epsilon_parts[1])))
epsilon_z_inv = sparse.diags(1 / epsilon_parts[2])
exy2ez = epsilon_z_inv @ exy2ez @ epsilon_xy
n_pts = dxes[0][0].size * dxes[0][1].size
op = sparse.vstack((sparse.eye(2 * n_pts),
exy2ez))
return op
h2eop = h2e(wavenumber, omega, dxes, epsilon)
return h2eop @ v2h(v, wavenumber, dxes, mu)
def e2h(wavenumber: complex,
@ -309,7 +185,7 @@ def e2h(wavenumber: complex,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Sparse matrix representation of the operator
"""
@ -330,7 +206,7 @@ def h2e(wavenumber: complex,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:return: Sparse matrix representation of the operator
"""
@ -343,7 +219,7 @@ 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)
:param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:return: Sparse matrix representation of the operator
"""
n = 1
@ -360,7 +236,7 @@ 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)
:param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:return: Sparse matrix representation of the operator
"""
n = 1
@ -385,7 +261,7 @@ def h_err(h: vfield_t,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Relative error norm(OP @ h) / norm(h)
@ -416,7 +292,7 @@ def e_err(e: vfield_t,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
:param epsilon: Vectorized dielectric constant grid
:param mu: Vectorized magnetic permeability grid (default 1 everywhere)
:return: Relative error norm(OP @ e) / norm(e)
@ -452,7 +328,7 @@ def cylindrical_operator(omega: complex,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (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.

249
meanas/fdfd/waveguide_mode.py → fdfd_tools/waveguide_mode.py
View File

@ -1,55 +1,73 @@
from typing import Dict, List, Tuple
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 vec, unvec, dx_lists_t, vfield_t, field_t
from . import operators, waveguide, functional
from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
from .eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
def vsolve_waveguide_mode_2d(mode_number: int,
def solve_waveguide_mode_2d(mode_number: int,
omega: complex,
dxes: dx_lists_t,
epsilon: vfield_t,
mu: vfield_t = None,
mode_margin: int = 2,
) -> Tuple[vfield_t, complex]:
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
:param epsilon: Dielectric constant
:param mu: Magnetic permeability (default 1 everywhere)
:param mode_margin: The eigensolver will actually solve for (mode_number + mode_margin)
modes, but only return the target mode. Increasing this value can improve the solver's
ability to find the correct mode. Default 2.
:return: (e_xy, wavenumber)
: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_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
A_r = waveguide.operator(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
eigvals, eigvecs = signed_eigensolve(A_r, mode_number + mode_margin)
e_xy = eigvecs[:, -(mode_number + 1)]
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_e(omega, dxes, epsilon, mu)
eigval, e_xy = rayleigh_quotient_iteration(A, e_xy)
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))
return e_xy, 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 beta increases.
'''
if wavenumber_correction:
wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - 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,
@ -60,6 +78,7 @@ def solve_waveguide_mode(mode_number: 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
@ -67,19 +86,19 @@ def solve_waveguide_mode(mode_number: int,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
: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)
mu = [numpy.ones_like(epsilon[0])] * 3
'''
Solve the 2D problem in the specified plane
@ -88,62 +107,57 @@ def solve_waveguide_mode(mode_number: int,
order = numpy.roll(range(3), 2 - axis)
reverse_order = numpy.roll(range(3), axis - 2)
# Find dx in propagation direction
dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes])
dx_prop = 0.5 * sum(dxab_forward)[0]
# Reduce to 2D and solve the 2D problem
args_2d = {
'omega': omega,
'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,
}
e_xy, wavenumber_2d = vsolve_waveguide_mode_2d(mode_number, **args_2d)
fields_2d = solve_waveguide_mode_2d(mode_number, omega=omega, **args_2d)
'''
Apply corrections and expand to 3D
'''
# Correct wavenumber to account for numerical dispersion.
wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2)
print(wavenumber_2d / wavenumber)
shape = [d.size for d in args_2d['dxes'][0]]
ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber)
e = unvec(ve, shape)
h = unvec(vh, shape)
# 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
h *= polarity
fields_2d['E'][2] *= polarity
fields_2d['H'][2] *= polarity
# Apply phase shift to H-field
h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
e[2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
d_prop = 0.5 * sum(dxab_forward)
for a in range(3):
fields_2d['H'][a] *= 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)
E = [None]*3
H = [None]*3
for a, o in enumerate(reverse_order):
E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order)
H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order)
E[a] = numpy.zeros_like(epsilon[0], dtype=complex)
H[a] = numpy.zeros_like(epsilon[0], dtype=complex)
E[a][slices] = fields_2d['E'][o][:, :, None].transpose(reverse_order)
H[a][slices] = fields_2d['H'][o][:, :, None].transpose(reverse_order)
results = {
'wavenumber': wavenumber,
'wavenumber_2d': wavenumber_2d,
'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],
epsilon: field_t,
mu: field_t = None,
) -> field_t:
"""
@ -151,9 +165,10 @@ def compute_source(E: field_t,
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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
: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
@ -161,34 +176,45 @@ def compute_source(E: field_t,
:param mu: Magnetic permeability (default 1 everywhere)
:return: J distribution for the unidirectional source
"""
E_expanded = expand_wgmode_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
polarity=polarity, slices=slices)
if mu is None:
mu = [1] * 3
smask = [slice(None)] * 4
if polarity > 0:
smask[axis + 1] = slice(slices[axis].start, None)
else:
smask[axis + 1] = slice(None, slices[axis].stop)
J = [None]*3
M = [None]*3
mask = numpy.zeros_like(E_expanded, dtype=int)
mask[tuple(smask)] = 1
src_order = numpy.roll(range(3), axis)
exp_iphi = numpy.exp(1j * polarity * wavenumber * dxes[1][axis][slices[axis]])
J[src_order[0]] = numpy.zeros_like(E[0])
J[src_order[1]] = +exp_iphi * H[src_order[2]] * polarity
J[src_order[2]] = -exp_iphi * H[src_order[1]] * polarity
M[src_order[0]] = numpy.zeros_like(E[0])
M[src_order[1]] = +numpy.roll(E[src_order[2]], -1, axis=axis)
M[src_order[2]] = -numpy.roll(E[src_order[1]], -1, axis=axis)
A1f = functional.curl_h(dxes)
Jm_iw = A1f([M[k] / mu[k] for k in range(3)])
for k in range(3):
J[k] += Jm_iw[k] / (-1j * omega)
masked_e2j = operators.e_boundary_source(mask=vec(mask), omega=omega, dxes=dxes, epsilon=vec(epsilon), mu=vec(mu))
J = unvec(masked_e2j @ vec(E_expanded), E.shape[1:])
return J
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],
) -> field_t: # TODO DOCS
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].i
[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.
@ -197,7 +223,7 @@ def compute_overlap_e(E: field_t,
: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 dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
: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
@ -205,22 +231,47 @@ def compute_overlap_e(E: field_t,
:param mu: Magnetic permeability (default 1 everywhere)
:return: overlap_e for calculating the mode overlap
"""
slices = tuple(slices)
cross_plane = [slice(None)] * 3
cross_plane[axis] = slices[axis]
Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes,
axis=axis, polarity=polarity, slices=slices)
# 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)
start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
# Expand our slice to the entire grid using the calculated phase factors
Ee = [None]*3
He = [None]*3
for k in range(3):
Ee[k] = phase_E * E[k][tuple(cross_plane)]
He[k] = phase_H * H[k][tuple(cross_plane)]
slices2 = list(slices)
slices2[axis] = slice(start, stop)
slices2 = (slice(None), *slices2)
Etgt = numpy.zeros_like(Ee)
Etgt[slices2] = Ee[slices2]
# Write out the operator product for the mode orthogonality integral
domain = numpy.zeros_like(E[0], dtype=int)
domain[slices] = 1
Etgt /= (Etgt.conj() * Etgt).sum()
return Etgt
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,
@ -228,19 +279,21 @@ def solve_waveguide_mode_cylindrical(mode_number: int,
dxes: dx_lists_t,
epsilon: vfield_t,
r0: float,
wavenumber_correction: bool = True,
) -> Dict[str, complex or field_t]:
"""
TODO: fixup
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.
:param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header.
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}
"""
@ -251,59 +304,37 @@ def solve_waveguide_mode_cylindrical(mode_number: int,
A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0)
eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3)
e_xy = eigvecs[:, -(mode_number+1)]
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, e_xy = rayleigh_quotient_iteration(A, e_xy)
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))
# TODO: Perform correction on wavenumber to account for numerical dispersion.
'''
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]]
e_xy = numpy.hstack((e_xy, numpy.zeros(shape[0] * shape[1])))
v = numpy.hstack((v, numpy.zeros(shape[0] * shape[1])))
fields = {
'wavenumber': wavenumber,
'E': unvec(e_xy, shape),
'E': unvec(v, shape),
# 'E': unvec(e, shape),
# 'H': unvec(h, shape),
}
return fields
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
E_expanded = numpy.zeros_like(E)
slices_exp = list(slices)
slices_exp[axis] = slice(E.shape[axis + 1])
slices_exp = (slice(None), *slices_exp)
slices_in = (slice(None), *slices)
E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
return E_expanded

1
meanas/VERSION
View File

@ -1 +0,0 @@
0.5

52
meanas/__init__.py
View File

@ -1,52 +0,0 @@
"""
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
"""
import pathlib
from .types import dx_lists_t, field_t, vfield_t, field_updater
from .vectorization import vec, unvec
__author__ = 'Jan Petykiewicz'
with open(pathlib.Path(__file__).parent / 'VERSION', 'r') as f:
__version__ = f.read().strip()

0
meanas/fdfd/__init__.py
View File

9
meanas/fdtd/__init__.py
View File

@ -1,9 +0,0 @@
"""
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
View File

@ -1,87 +0,0 @@
"""
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 meanas.types
: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 meanas.types
: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
View File

@ -1,68 +0,0 @@
"""
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
View File

@ -1,84 +0,0 @@
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
View File

@ -1,122 +0,0 @@
"""
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

0
meanas/test/__init__.py
View File

293
meanas/test/test_fdtd.py
View File

@ -1,293 +0,0 @@
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
View File

@ -1,22 +0,0 @@
"""
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]

33
setup.py
View File

@ -1,41 +1,18 @@
#!/usr/bin/env python3
#!/usr/bin/env python
from setuptools import setup, find_packages
with open('README.md', 'r') as f:
long_description = f.read()
with open('meanas/VERSION', 'r') as f:
version = f.read().strip()
setup(name='meanas',
version=version,
description='Electromagnetic simulation tools',
long_description=long_description,
long_description_content_type='text/markdown',
setup(name='fdfd_tools',
version='0.4',
description='FDFD Electromagnetic simulation tools',
author='Jan Petykiewicz',
author_email='anewusername@gmail.com',
url='https://mpxd.net/code/jan/fdfd_tools',
url='https://mpxd.net/gogs/jan/fdfd_tools',
packages=find_packages(),
package_data={
'meanas': ['VERSION']
},
install_requires=[
'numpy',
'scipy',
],
extras_require={
'test': [
'pytest',
'dataclasses',
],
},
classifiers=[
'Programming Language :: Python :: 3',
'Development Status :: 4 - Beta',
'Intended Audience :: Developers',
'Intended Audience :: Science/Research',
'License :: OSI Approved :: GNU Affero General Public License v3',
'Topic :: Scientific/Engineering :: Physics',
],
)