Compare commits
92 Commits
Author | SHA1 | Date |
---|---|---|
Jan Petykiewicz | 6850fe532f | 5 years ago |
Jan Petykiewicz | 6333dbd110 | 5 years ago |
Jan Petykiewicz | 0ad289e5ac | 5 years ago |
Jan Petykiewicz | 487bdd61ec | 5 years ago |
jan | a1a7aa5511 | 5 years ago |
jan | 9cd2dd131b | 5 years ago |
jan | c92656bed8 | 5 years ago |
jan | 10e275611d | 5 years ago |
jan | 2289c6d116 | 5 years ago |
jan | bedec489aa | 5 years ago |
jan | f04c0daf28 | 5 years ago |
jan | b5ad284966 | 5 years ago |
Jan Petykiewicz | 8eac9df76e | 5 years ago |
Jan Petykiewicz | a4c2239ad9 | 5 years ago |
Jan Petykiewicz | e99019b37f | 5 years ago |
Jan Petykiewicz | f4bac9598d | 5 years ago |
Jan Petykiewicz | d2d4220313 | 5 years ago |
Jan Petykiewicz | 337cee8018 | 5 years ago |
Jan Petykiewicz | 5f96474497 | 5 years ago |
Jan Petykiewicz | 7b56aa363f | 5 years ago |
Jan Petykiewicz | 3887a8804f | 5 years ago |
Jan Petykiewicz | d6a34b280e | 5 years ago |
Jan Petykiewicz | 1860d754cd | 5 years ago |
Jan Petykiewicz | 7006b5e6e4 | 5 years ago |
Jan Petykiewicz | c306bb1f46 | 5 years ago |
Jan Petykiewicz | af8efd00eb | 5 years ago |
Jan Petykiewicz | 41bec05d4e | 5 years ago |
Jan Petykiewicz | 2787908640 | 5 years ago |
Jan Petykiewicz | 054ac994d5 | 5 years ago |
Jan Petykiewicz | 0503e9d6ef | 5 years ago |
Jan Petykiewicz | b466ed02ea | 5 years ago |
Jan Petykiewicz | ccdb423ba2 | 5 years ago |
Jan Petykiewicz | aade81c21e | 5 years ago |
Jan Petykiewicz | 07c94617fe | 5 years ago |
Jan Petykiewicz | 1a04bab361 | 5 years ago |
Jan Petykiewicz | 2c91ea249f | 5 years ago |
Jan Petykiewicz | 3429120993 | 5 years ago |
Jan Petykiewicz | 938c4c9a35 | 5 years ago |
Jan Petykiewicz | 5951f2bdb1 | 5 years ago |
Jan Petykiewicz | 94ff3f7853 | 5 years ago |
Jan Petykiewicz | f61bcf3dfa | 5 years ago |
Jan Petykiewicz | 25cb83089d | 5 years ago |
Jan Petykiewicz | 3d07969fd2 | 5 years ago |
Jan Petykiewicz | 557a3b0d9c | 5 years ago |
Jan Petykiewicz | 32055ec8d3 | 5 years ago |
Jan Petykiewicz | 06a491a960 | 5 years ago |
Jan Petykiewicz | 1489308837 | 5 years ago |
Jan Petykiewicz | 1793e8cc37 | 5 years ago |
Jan Petykiewicz | 1d9c9644ee | 5 years ago |
Jan Petykiewicz | 56a1349959 | 5 years ago |
Jan Petykiewicz | f2d061c921 | 5 years ago |
Jan Petykiewicz | b4bbfdb730 | 5 years ago |
Jan Petykiewicz | 39c05d2cab | 5 years ago |
Jan Petykiewicz | 89976647f2 | 5 years ago |
Jan Petykiewicz | 7092c13088 | 5 years ago |
Jan Petykiewicz | f1fc308d25 | 5 years ago |
Jan Petykiewicz | 7f8a326114 | 5 years ago |
Jan Petykiewicz | 30ddeb7b73 | 5 years ago |
Jan Petykiewicz | 223b202d03 | 5 years ago |
Jan Petykiewicz | fb3c88a78d | 5 years ago |
Jan Petykiewicz | 950e70213a | 5 years ago |
Jan Petykiewicz | f858cb8bbb | 5 years ago |
Jan Petykiewicz | 2cec4fabaf | 5 years ago |
Jan Petykiewicz | a528effd89 | 5 years ago |
Jan Petykiewicz | 935b2c9a80 | 5 years ago |
Jan Petykiewicz | 79e14af4db | 5 years ago |
Jan Petykiewicz | dd4e6f294f | 5 years ago |
Jan Petykiewicz | ecaf9fa3d0 | 5 years ago |
Jan Petykiewicz | 099966f291 | 5 years ago |
Jan Petykiewicz | a8a5a69e1a | 5 years ago |
Jan Petykiewicz | 557e748356 | 5 years ago |
Jan Petykiewicz | 9d1d8fe869 | 5 years ago |
Jan Petykiewicz | 8e634e35df | 5 years ago |
Jan Petykiewicz | 4c2035c882 | 5 years ago |
Jan Petykiewicz | d462ae9412 | 5 years ago |
Jan Petykiewicz | 2acbda4764 | 5 years ago |
Jan Petykiewicz | 3a5d75cde4 | 5 years ago |
Jan Petykiewicz | 2b3a74b737 | 5 years ago |
Jan Petykiewicz | 5dd26915fc | 5 years ago |
Jan Petykiewicz | c3f248a73c | 5 years ago |
Jan Petykiewicz | 001c32a2e0 | 5 years ago |
Jan Petykiewicz | 41cd94fe48 | 5 years ago |
Jan Petykiewicz | c7d4c4a8e6 | 5 years ago |
jan | 1f9a9949c0 | 6 years ago |
jan | 323bcf88ad | 6 years ago |
jan | ee9abb77d9 | 6 years ago |
jan | c1f65f61c1 | 6 years ago |
jan | e8f836c908 | 6 years ago |
jan | 0e47fdd5fb | 6 years ago |
jan | e02040c709 | 6 years ago |
jan | c4cbdff751 | 6 years ago |
jan | 4067766478 | 6 years ago |
@ -0,0 +1,3 @@
|
|||||||
|
include README.md
|
||||||
|
include LICENSE.md
|
||||||
|
include meanas/VERSION
|
@ -1,46 +1,56 @@
|
|||||||
# fdfd_tools
|
# meanas
|
||||||
|
|
||||||
**fdfd_tools** is a python package containing utilities for
|
**meanas** is a python package for electromagnetic simulations
|
||||||
creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD)
|
|
||||||
electromagnetic simulations.
|
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**
|
**Contents**
|
||||||
* Library of sparse matrices for representing the electromagnetic wave
|
- Finite difference frequency domain (FDFD)
|
||||||
equation in 3D, as well as auxiliary matrices for conversion between fields
|
* Library of sparse matrices for representing the electromagnetic wave
|
||||||
* Waveguide mode solver and waveguide mode operators
|
equation in 3D, as well as auxiliary matrices for conversion between fields
|
||||||
* Stretched-coordinate PML boundaries (SCPML)
|
* Waveguide mode operators
|
||||||
* Functional versions of most operators
|
* Waveguide mode eigensolver
|
||||||
* Anisotropic media (eps_xx, eps_yy, eps_zz, mu_xx, ...)
|
* Stretched-coordinate PML boundaries (SCPML)
|
||||||
* Arbitrary distributions of perfect electric and magnetic conductors (PEC / PMC)
|
* 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
|
This package does *not* provide a fast matrix solver, though by default
|
||||||
```fdfd_tools.solvers.generic(...)``` will call
|
`meanas.fdfd.solvers.generic(...)` will call
|
||||||
```scipy.sparse.linalg.qmr(...)``` to perform a solve.
|
`scipy.sparse.linalg.qmr(...)` to perform a solve.
|
||||||
For 2D problems this should be fine; likewise, the waveguide mode
|
For 2D FDFD problems this should be fine; likewise, the waveguide mode
|
||||||
solver uses scipy's eigenvalue solver, with reasonable results.
|
solver uses scipy's eigenvalue solver, with reasonable results.
|
||||||
|
|
||||||
For solving large (or 3D) problems, I recommend a GPU-based iterative
|
For solving large (or 3D) FDFD problems, I recommend a GPU-based iterative
|
||||||
solver, such as [opencl_fdfd](https://mpxd.net/gogs/jan/opencl_fdfd) or
|
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
|
those included in [MAGMA](http://icl.cs.utk.edu/magma/index.html)). Your
|
||||||
solver will need the ability to solve complex symmetric (non-Hermitian)
|
solver will need the ability to solve complex symmetric (non-Hermitian)
|
||||||
linear systems, ideally with double precision.
|
linear systems, ideally with double precision.
|
||||||
|
|
||||||
|
|
||||||
## Installation
|
## Installation
|
||||||
|
|
||||||
**Requirements:**
|
**Requirements:**
|
||||||
* python 3 (written and tested with 3.5)
|
* python 3 (tests require 3.7)
|
||||||
* numpy
|
* numpy
|
||||||
* scipy
|
* scipy
|
||||||
|
|
||||||
|
|
||||||
Install with pip, via git:
|
Install with pip, via git:
|
||||||
```bash
|
```bash
|
||||||
pip install git+https://mpxd.net/gogs/jan/fdfd_tools.git@release
|
pip install git+https://mpxd.net/code/jan/meanas.git@release
|
||||||
```
|
```
|
||||||
|
|
||||||
## Use
|
## Use
|
||||||
|
|
||||||
See examples/test.py for some simple examples; you may need additional
|
See `examples/` for some simple examples; you may need additional
|
||||||
packages such as [gridlock](https://mpxd.net/gogs/jan/gridlock)
|
packages such as [gridlock](https://mpxd.net/code/jan/gridlock)
|
||||||
to run the examples.
|
to run the examples.
|
||||||
|
@ -0,0 +1,90 @@
|
|||||||
|
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()
|
@ -1,25 +0,0 @@
|
|||||||
"""
|
|
||||||
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'
|
|
@ -1,239 +0,0 @@
|
|||||||
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
|
|
@ -0,0 +1 @@
|
|||||||
|
0.5
|
@ -0,0 +1,52 @@
|
|||||||
|
"""
|
||||||
|
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,0 +1,9 @@
|
|||||||
|
"""
|
||||||
|
Basic FDTD functionality
|
||||||
|
"""
|
||||||
|
|
||||||
|
from .base import maxwell_e, maxwell_h
|
||||||
|
from .pml import cpml
|
||||||
|
from .energy import (poynting, poynting_divergence, energy_hstep, energy_estep,
|
||||||
|
delta_energy_h2e, delta_energy_h2e, delta_energy_j)
|
||||||
|
from .boundaries import conducting_boundary
|
@ -0,0 +1,87 @@
|
|||||||
|
"""
|
||||||
|
Basic FDTD field updates
|
||||||
|
"""
|
||||||
|
from typing import List, Callable, Tuple, Dict
|
||||||
|
import numpy
|
||||||
|
|
||||||
|
from .. import dx_lists_t, field_t, field_updater
|
||||||
|
|
||||||
|
__author__ = 'Jan Petykiewicz'
|
||||||
|
|
||||||
|
|
||||||
|
def curl_h(dxes: dx_lists_t = None) -> field_updater:
|
||||||
|
"""
|
||||||
|
Curl operator for use with the H field.
|
||||||
|
|
||||||
|
:param dxes: Grid parameters [dx_e, dx_h] as described in 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
|
@ -0,0 +1,68 @@
|
|||||||
|
"""
|
||||||
|
Boundary conditions
|
||||||
|
"""
|
||||||
|
|
||||||
|
from typing import List, Callable, Tuple, Dict
|
||||||
|
import numpy
|
||||||
|
|
||||||
|
from .. import dx_lists_t, field_t, field_updater
|
||||||
|
|
||||||
|
|
||||||
|
def conducting_boundary(direction: int,
|
||||||
|
polarity: int
|
||||||
|
) -> Tuple[field_updater, field_updater]:
|
||||||
|
dirs = [0, 1, 2]
|
||||||
|
if direction not in dirs:
|
||||||
|
raise Exception('Invalid direction: {}'.format(direction))
|
||||||
|
dirs.remove(direction)
|
||||||
|
u, v = dirs
|
||||||
|
|
||||||
|
if polarity < 0:
|
||||||
|
boundary_slice = [slice(None)] * 3
|
||||||
|
shifted1_slice = [slice(None)] * 3
|
||||||
|
boundary_slice[direction] = 0
|
||||||
|
shifted1_slice[direction] = 1
|
||||||
|
|
||||||
|
def en(e: field_t):
|
||||||
|
e[direction][boundary_slice] = 0
|
||||||
|
e[u][boundary_slice] = e[u][shifted1_slice]
|
||||||
|
e[v][boundary_slice] = e[v][shifted1_slice]
|
||||||
|
return e
|
||||||
|
|
||||||
|
def hn(h: field_t):
|
||||||
|
h[direction][boundary_slice] = h[direction][shifted1_slice]
|
||||||
|
h[u][boundary_slice] = 0
|
||||||
|
h[v][boundary_slice] = 0
|
||||||
|
return h
|
||||||
|
|
||||||
|
return en, hn
|
||||||
|
|
||||||
|
elif polarity > 0:
|
||||||
|
boundary_slice = [slice(None)] * 3
|
||||||
|
shifted1_slice = [slice(None)] * 3
|
||||||
|
shifted2_slice = [slice(None)] * 3
|
||||||
|
boundary_slice[direction] = -1
|
||||||
|
shifted1_slice[direction] = -2
|
||||||
|
shifted2_slice[direction] = -3
|
||||||
|
|
||||||
|
def ep(e: field_t):
|
||||||
|
e[direction][boundary_slice] = -e[direction][shifted2_slice]
|
||||||
|
e[direction][shifted1_slice] = 0
|
||||||
|
e[u][boundary_slice] = e[u][shifted1_slice]
|
||||||
|
e[v][boundary_slice] = e[v][shifted1_slice]
|
||||||
|
return e
|
||||||
|
|
||||||
|
def hp(h: field_t):
|
||||||
|
h[direction][boundary_slice] = h[direction][shifted1_slice]
|
||||||
|
h[u][boundary_slice] = -h[u][shifted2_slice]
|
||||||
|
h[u][shifted1_slice] = 0
|
||||||
|
h[v][boundary_slice] = -h[v][shifted2_slice]
|
||||||
|
h[v][shifted1_slice] = 0
|
||||||
|
return h
|
||||||
|
|
||||||
|
return ep, hp
|
||||||
|
|
||||||
|
else:
|
||||||
|
raise Exception('Bad polarity: {}'.format(polarity))
|
||||||
|
|
||||||
|
|
@ -0,0 +1,84 @@
|
|||||||
|
from typing import List, Callable, Tuple, Dict
|
||||||
|
import numpy
|
||||||
|
|
||||||
|
from .. import dx_lists_t, field_t, field_updater
|
||||||
|
|
||||||
|
|
||||||
|
def poynting(e, h):
|
||||||
|
s = (numpy.roll(e[1], -1, axis=0) * h[2] - numpy.roll(e[2], -1, axis=0) * h[1],
|
||||||
|
numpy.roll(e[2], -1, axis=1) * h[0] - numpy.roll(e[0], -1, axis=1) * h[2],
|
||||||
|
numpy.roll(e[0], -1, axis=2) * h[1] - numpy.roll(e[1], -1, axis=2) * h[0])
|
||||||
|
return numpy.array(s)
|
||||||
|
|
||||||
|
|
||||||
|
def poynting_divergence(s=None, *, e=None, h=None, dxes=None): # TODO dxes
|
||||||
|
if dxes is None:
|
||||||
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
||||||
|
|
||||||
|
if s is None:
|
||||||
|
s = poynting(e, h)
|
||||||
|
|
||||||
|
ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) / numpy.sqrt(dxes[0][0] * dxes[1][0])[:, None, None] +
|
||||||
|
(s[1] - numpy.roll(s[1], 1, axis=1)) / numpy.sqrt(dxes[0][1] * dxes[1][1])[None, :, None] +
|
||||||
|
(s[2] - numpy.roll(s[2], 1, axis=2)) / numpy.sqrt(dxes[0][2] * dxes[1][2])[None, None, :] )
|
||||||
|
return ds
|
||||||
|
|
||||||
|
|
||||||
|
def energy_hstep(e0, h1, e2, epsilon=None, mu=None, dxes=None):
|
||||||
|
u = dxmul(e0 * e2, h1 * h1, epsilon, mu, dxes)
|
||||||
|
return u
|
||||||
|
|
||||||
|
|
||||||
|
def energy_estep(h0, e1, h2, epsilon=None, mu=None, dxes=None):
|
||||||
|
u = dxmul(e1 * e1, h0 * h2, epsilon, mu, dxes)
|
||||||
|
return u
|
||||||
|
|
||||||
|
|
||||||
|
def delta_energy_h2e(dt, e0, h1, e2, h3, epsilon=None, mu=None, dxes=None):
|
||||||
|
"""
|
||||||
|
This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2)
|
||||||
|
"""
|
||||||
|
de = e2 * (e2 - e0) / dt
|
||||||
|
dh = h1 * (h3 - h1) / dt
|
||||||
|
du = dxmul(de, dh, epsilon, mu, dxes)
|
||||||
|
return du
|
||||||
|
|
||||||
|
|
||||||
|
def delta_energy_e2h(dt, h0, e1, h2, e3, epsilon=None, mu=None, dxes=None):
|
||||||
|
"""
|
||||||
|
This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2)
|
||||||
|
"""
|
||||||
|
de = e1 * (e3 - e1) / dt
|
||||||
|
dh = h2 * (h2 - h0) / dt
|
||||||
|
du = dxmul(de, dh, epsilon, mu, dxes)
|
||||||
|
return du
|
||||||
|
|
||||||
|
|
||||||
|
def delta_energy_j(j0, e1, dxes=None):
|
||||||
|
if dxes is None:
|
||||||
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
||||||
|
|
||||||
|
du = ((j0 * e1).sum(axis=0) *
|
||||||
|
dxes[0][0][:, None, None] *
|
||||||
|
dxes[0][1][None, :, None] *
|
||||||
|
dxes[0][2][None, None, :])
|
||||||
|
return du
|
||||||
|
|
||||||
|
|
||||||
|
def dxmul(ee, hh, epsilon=None, mu=None, dxes=None):
|
||||||
|
if epsilon is None:
|
||||||
|
epsilon = 1
|
||||||
|
if mu is None:
|
||||||
|
mu = 1
|
||||||
|
if dxes is None:
|
||||||
|
dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
|
||||||
|
|
||||||
|
result = ((ee * epsilon).sum(axis=0) *
|
||||||
|
dxes[0][0][:, None, None] *
|
||||||
|
dxes[0][1][None, :, None] *
|
||||||
|
dxes[0][2][None, None, :] +
|
||||||
|
(hh * mu).sum(axis=0) *
|
||||||
|
dxes[1][0][:, None, None] *
|
||||||
|
dxes[1][1][None, :, None] *
|
||||||
|
dxes[1][2][None, None, :])
|
||||||
|
return result
|
@ -0,0 +1,122 @@
|
|||||||
|
"""
|
||||||
|
PML implementations
|
||||||
|
|
||||||
|
"""
|
||||||
|
# TODO retest pmls!
|
||||||
|
|
||||||
|
from typing import List, Callable, Tuple, Dict
|
||||||
|
import numpy
|
||||||
|
|
||||||
|
from .. import dx_lists_t, field_t, field_updater
|
||||||
|
|
||||||
|
|
||||||
|
__author__ = 'Jan Petykiewicz'
|
||||||
|
|
||||||
|
|
||||||
|
def cpml(direction:int,
|
||||||
|
polarity: int,
|
||||||
|
dt: float,
|
||||||
|
epsilon: field_t,
|
||||||
|
thickness: int = 8,
|
||||||
|
ln_R_per_layer: float = -1.6,
|
||||||
|
epsilon_eff: float = 1,
|
||||||
|
mu_eff: float = 1,
|
||||||
|
m: float = 3.5,
|
||||||
|
ma: float = 1,
|
||||||
|
cfs_alpha: float = 0,
|
||||||
|
dtype: numpy.dtype = numpy.float32,
|
||||||
|
) -> Tuple[Callable, Callable, Dict[str, field_t]]:
|
||||||
|
|
||||||
|
if direction not in range(3):
|
||||||
|
raise Exception('Invalid direction: {}'.format(direction))
|
||||||
|
|
||||||
|
if polarity not in (-1, 1):
|
||||||
|
raise Exception('Invalid polarity: {}'.format(polarity))
|
||||||
|
|
||||||
|
if thickness <= 2:
|
||||||
|
raise Exception('It would be wise to have a pml with 4+ cells of thickness')
|
||||||
|
|
||||||
|
if epsilon_eff <= 0:
|
||||||
|
raise Exception('epsilon_eff must be positive')
|
||||||
|
|
||||||
|
sigma_max = -ln_R_per_layer / 2 * (m + 1)
|
||||||
|
kappa_max = numpy.sqrt(epsilon_eff * mu_eff)
|
||||||
|
alpha_max = cfs_alpha
|
||||||
|
transverse = numpy.delete(range(3), direction)
|
||||||
|
u, v = transverse
|
||||||
|
|
||||||
|
xe = numpy.arange(1, thickness+1, dtype=float)
|
||||||
|
xh = numpy.arange(1, thickness+1, dtype=float)
|
||||||
|
if polarity > 0:
|
||||||
|
xe -= 0.5
|
||||||
|
elif polarity < 0:
|
||||||
|
xh -= 0.5
|
||||||
|
xe = xe[::-1]
|
||||||
|
xh = xh[::-1]
|
||||||
|
else:
|
||||||
|
raise Exception('Bad polarity!')
|
||||||
|
|
||||||
|
expand_slice = [None] * 3
|
||||||
|
expand_slice[direction] = slice(None)
|
||||||
|
|
||||||
|
def par(x):
|
||||||
|
scaling = (x / thickness) ** m
|
||||||
|
sigma = scaling * sigma_max
|
||||||
|
kappa = 1 + scaling * (kappa_max - 1)
|
||||||
|
alpha = ((1 - x / thickness) ** ma) * alpha_max
|
||||||
|
p0 = numpy.exp(-(sigma / kappa + alpha) * dt)
|
||||||
|
p1 = sigma / (sigma + kappa * alpha) * (p0 - 1)
|
||||||
|
p2 = 1 / kappa
|
||||||
|
return p0[expand_slice], p1[expand_slice], p2[expand_slice]
|
||||||
|
|
||||||
|
p0e, p1e, p2e = par(xe)
|
||||||
|
p0h, p1h, p2h = par(xh)
|
||||||
|
|
||||||
|
region = [slice(None)] * 3
|
||||||
|
if polarity < 0:
|
||||||
|
region[direction] = slice(None, thickness)
|
||||||
|
elif polarity > 0:
|
||||||
|
region[direction] = slice(-thickness, None)
|
||||||
|
else:
|
||||||
|
raise Exception('Bad polarity!')
|
||||||
|
|
||||||
|
se = 1 if direction == 1 else -1
|
||||||
|
|
||||||
|
# TODO check if epsilon is uniform in pml region?
|
||||||
|
shape = list(epsilon[0].shape)
|
||||||
|
shape[direction] = thickness
|
||||||
|
psi_e = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)]
|
||||||
|
psi_h = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)]
|
||||||
|
|
||||||
|
fields = {
|
||||||
|
'psi_e_u': psi_e[0],
|
||||||
|
'psi_e_v': psi_e[1],
|
||||||
|
'psi_h_u': psi_h[0],
|
||||||
|
'psi_h_v': psi_h[1],
|
||||||
|
}
|
||||||
|
|
||||||
|
# Note that this is kinda slow -- would be faster to reuse dHv*p2h for the original
|
||||||
|
# H update, but then you have multiple arrays and a monolithic (field + pml) update operation
|
||||||
|
def pml_e(e: field_t, h: field_t, epsilon: field_t) -> Tuple[field_t, field_t]:
|
||||||
|
dHv = h[v][region] - numpy.roll(h[v], 1, axis=direction)[region]
|
||||||
|
dHu = h[u][region] - numpy.roll(h[u], 1, axis=direction)[region]
|
||||||
|
psi_e[0] *= p0e
|
||||||
|
psi_e[0] += p1e * dHv * p2e
|
||||||
|
psi_e[1] *= p0e
|
||||||
|
psi_e[1] += p1e * dHu * p2e
|
||||||
|
e[u][region] += se * dt / epsilon[u][region] * (psi_e[0] + (p2e - 1) * dHv)
|
||||||
|
e[v][region] -= se * dt / epsilon[v][region] * (psi_e[1] + (p2e - 1) * dHu)
|
||||||
|
return e, h
|
||||||
|
|
||||||
|
def pml_h(e: field_t, h: field_t) -> Tuple[field_t, field_t]:
|
||||||
|
dEv = (numpy.roll(e[v], -1, axis=direction)[region] - e[v][region])
|
||||||
|
dEu = (numpy.roll(e[u], -1, axis=direction)[region] - e[u][region])
|
||||||
|
psi_h[0] *= p0h
|
||||||
|
psi_h[0] += p1h * dEv * p2h
|
||||||
|
psi_h[1] *= p0h
|
||||||
|
psi_h[1] += p1h * dEu * p2h
|
||||||
|
h[u][region] -= se * dt * (psi_h[0] + (p2h - 1) * dEv)
|
||||||
|
h[v][region] += se * dt * (psi_h[1] + (p2h - 1) * dEu)
|
||||||
|
return e, h
|
||||||
|
|
||||||
|
return pml_e, pml_h, fields
|
@ -0,0 +1,293 @@
|
|||||||
|
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
|
||||||
|
|
||||||
|
|
@ -0,0 +1,22 @@
|
|||||||
|
"""
|
||||||
|
Types shared across multiple submodules
|
||||||
|
"""
|
||||||
|
import numpy
|
||||||
|
from typing import List, Callable
|
||||||
|
|
||||||
|
|
||||||
|
# Field types
|
||||||
|
field_t = numpy.ndarray # vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])
|
||||||
|
vfield_t = numpy.ndarray # linearized vector field (vector of length 3*X*Y*Z)
|
||||||
|
|
||||||
|
'''
|
||||||
|
'dxes' datastructure which contains grid cell width information in the following format:
|
||||||
|
[[[dx_e_0, dx_e_1, ...], [dy_e_0, ...], [dz_e_0, ...]],
|
||||||
|
[[dx_h_0, dx_h_1, ...], [dy_h_0, ...], [dz_h_0, ...]]]
|
||||||
|
where dx_e_0 is the x-width of the x=0 cells, as used when calculating dE/dx,
|
||||||
|
and dy_h_0 is the y-width of the y=0 cells, as used when calculating dH/dy, etc.
|
||||||
|
'''
|
||||||
|
dx_lists_t = List[List[numpy.ndarray]]
|
||||||
|
|
||||||
|
|
||||||
|
field_updater = Callable[[field_t], field_t]
|
@ -1,18 +1,41 @@
|
|||||||
#!/usr/bin/env python
|
#!/usr/bin/env python3
|
||||||
|
|
||||||
from setuptools import setup, find_packages
|
from setuptools import setup, find_packages
|
||||||
|
|
||||||
setup(name='fdfd_tools',
|
with open('README.md', 'r') as f:
|
||||||
version='0.4',
|
long_description = f.read()
|
||||||
description='FDFD Electromagnetic simulation tools',
|
|
||||||
|
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',
|
||||||
author='Jan Petykiewicz',
|
author='Jan Petykiewicz',
|
||||||
author_email='anewusername@gmail.com',
|
author_email='anewusername@gmail.com',
|
||||||
url='https://mpxd.net/gogs/jan/fdfd_tools',
|
url='https://mpxd.net/code/jan/fdfd_tools',
|
||||||
packages=find_packages(),
|
packages=find_packages(),
|
||||||
|
package_data={
|
||||||
|
'meanas': ['VERSION']
|
||||||
|
},
|
||||||
install_requires=[
|
install_requires=[
|
||||||
'numpy',
|
'numpy',
|
||||||
'scipy',
|
'scipy',
|
||||||
],
|
],
|
||||||
extras_require={
|
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',
|
||||||
|
],
|
||||||
)
|
)
|
||||||
|
Loading…
Reference in New Issue