From f61bcf3dfa6874d09b5c4eecc9774b3cd6aaeb7c Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 4 Aug 2019 13:48:41 -0700 Subject: [PATCH] rename to meanas and split fdtd/fdfd --- README.md | 50 +-- fdfd_tools/__init__.py | 26 -- fdfd_tools/fdtd.py | 339 ------------------ meanas/__init__.py | 48 +++ {fdfd_tools => meanas}/eigensolvers.py | 0 {fdfd_tools => meanas/fdfd}/bloch.py | 0 {fdfd_tools => meanas/fdfd}/farfield.py | 0 {fdfd_tools => meanas/fdfd}/functional.py | 16 +- {fdfd_tools => meanas/fdfd}/operators.py | 30 +- fdfd_tools/grid.py => meanas/fdfd/scpml.py | 1 - {fdfd_tools => meanas/fdfd}/solvers.py | 2 +- {fdfd_tools => meanas/fdfd}/waveguide.py | 24 +- {fdfd_tools => meanas/fdfd}/waveguide_mode.py | 10 +- meanas/fdtd/__init__.py | 9 + meanas/fdtd/base.py | 87 +++++ meanas/fdtd/boundaries.py | 68 ++++ meanas/fdtd/energy.py | 84 +++++ meanas/fdtd/pml.py | 122 +++++++ {fdfd_tools => meanas}/test/test_fdtd.py | 2 +- meanas/types.py | 22 ++ {fdfd_tools => meanas}/vectorization.py | 10 +- setup.py | 8 +- 22 files changed, 519 insertions(+), 439 deletions(-) delete mode 100644 fdfd_tools/__init__.py delete mode 100644 fdfd_tools/fdtd.py create mode 100644 meanas/__init__.py rename {fdfd_tools => meanas}/eigensolvers.py (100%) rename {fdfd_tools => meanas/fdfd}/bloch.py (100%) rename {fdfd_tools => meanas/fdfd}/farfield.py (100%) rename {fdfd_tools => meanas/fdfd}/functional.py (87%) rename {fdfd_tools => meanas/fdfd}/operators.py (92%) rename fdfd_tools/grid.py => meanas/fdfd/scpml.py (99%) rename {fdfd_tools => meanas/fdfd}/solvers.py (97%) rename {fdfd_tools => meanas/fdfd}/waveguide.py (92%) rename {fdfd_tools => meanas/fdfd}/waveguide_mode.py (97%) create mode 100644 meanas/fdtd/__init__.py create mode 100644 meanas/fdtd/base.py create mode 100644 meanas/fdtd/boundaries.py create mode 100644 meanas/fdtd/energy.py create mode 100644 meanas/fdtd/pml.py rename {fdfd_tools => meanas}/test/test_fdtd.py (99%) create mode 100644 meanas/types.py rename {fdfd_tools => meanas}/vectorization.py (87%) diff --git a/README.md b/README.md index 5a3f49c..3ead0ca 100644 --- a/README.md +++ b/README.md @@ -1,46 +1,56 @@ -# fdfd_tools +# meanas -**fdfd_tools** is a python package containing utilities for -creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD) -electromagnetic simulations. +**meanas** is a python package for 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** -* Library of sparse matrices for representing the electromagnetic wave - equation in 3D, as well as auxiliary matrices for conversion between fields -* 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) +- 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 -```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 +`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) problems, I recommend a GPU-based iterative -solver, such as [opencl_fdfd](https://mpxd.net/gogs/jan/opencl_fdfd) or +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. + ## Installation **Requirements:** -* python 3 (written and tested with 3.5) +* python 3 (tests require 3.7) * numpy * scipy Install with pip, via git: ```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 -See examples/test.py for some simple examples; you may need additional -packages such as [gridlock](https://mpxd.net/gogs/jan/gridlock) +See `examples/` for some simple examples; you may need additional +packages such as [gridlock](https://mpxd.net/code/jan/gridlock) to run the examples. diff --git a/fdfd_tools/__init__.py b/fdfd_tools/__init__.py deleted file mode 100644 index ecf4e15..0000000 --- a/fdfd_tools/__init__.py +++ /dev/null @@ -1,26 +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' -version = '0.5' diff --git a/fdfd_tools/fdtd.py b/fdfd_tools/fdtd.py deleted file mode 100644 index 15b5635..0000000 --- a/fdfd_tools/fdtd.py +++ /dev/null @@ -1,339 +0,0 @@ -from typing import List, Callable, Tuple, Dict -import numpy - -from . import dx_lists_t, field_t - -#TODO fix pmls - -__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: - 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) -> 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: - 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) -> functional_matrix: - 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) -> functional_matrix: - 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 - - -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, - 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 - - -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 - - - diff --git a/meanas/__init__.py b/meanas/__init__.py new file mode 100644 index 0000000..d4288a5 --- /dev/null +++ b/meanas/__init__.py @@ -0,0 +1,48 @@ +""" +Electromagnetic simulation tools + +This package is intended for building simulation inputs, analyzing +simulation outputs, and running short simulations on unspecialized hardware. +It is designed to provide tooling and a baseline for other, high-performance +purpose- and hardware-specific solvers. + + +**Contents** +- Finite difference frequency domain (FDFD) + * Library of sparse matrices for representing the electromagnetic wave + equation in 3D, as well as auxiliary matrices for conversion between fields + * Waveguide mode operators + * Waveguide mode eigensolver + * Stretched-coordinate PML boundaries (SCPML) + * Functional versions of most operators + * Anisotropic media (limited to diagonal elements eps_xx, eps_yy, eps_zz, mu_xx, ...) + * Arbitrary distributions of perfect electric and magnetic conductors (PEC / PMC) +- Finite difference time domain (FDTD) + * Basic Maxwell time-steps + * Poynting vector and energy calculation + * Convolutional PMLs + +This package does *not* provide a fast matrix solver, though by default +```meanas.fdfd.solvers.generic(...)``` will call +```scipy.sparse.linalg.qmr(...)``` to perform a solve. +For 2D FDFD problems this should be fine; likewise, the waveguide mode +solver uses scipy's eigenvalue solver, with reasonable results. + +For solving large (or 3D) FDFD problems, I recommend a GPU-based iterative +solver, such as [opencl_fdfd](https://mpxd.net/code/jan/opencl_fdfd) or +those included in [MAGMA](http://icl.cs.utk.edu/magma/index.html)). Your +solver will need the ability to solve complex symmetric (non-Hermitian) +linear systems, ideally with double precision. + + +Dependencies: +- numpy +- scipy + +""" + +from .types import dx_lists_t, field_t, vfield_t, field_updater +from .vectorization import vec, unvec + +__author__ = 'Jan Petykiewicz' +version = '0.5' diff --git a/fdfd_tools/eigensolvers.py b/meanas/eigensolvers.py similarity index 100% rename from fdfd_tools/eigensolvers.py rename to meanas/eigensolvers.py diff --git a/fdfd_tools/bloch.py b/meanas/fdfd/bloch.py similarity index 100% rename from fdfd_tools/bloch.py rename to meanas/fdfd/bloch.py diff --git a/fdfd_tools/farfield.py b/meanas/fdfd/farfield.py similarity index 100% rename from fdfd_tools/farfield.py rename to meanas/fdfd/farfield.py diff --git a/fdfd_tools/functional.py b/meanas/fdfd/functional.py similarity index 87% rename from fdfd_tools/functional.py rename to meanas/fdfd/functional.py index 1d39d84..e57fe88 100644 --- a/fdfd_tools/functional.py +++ b/meanas/fdfd/functional.py @@ -2,8 +2,8 @@ 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 inputs in the form E = [E_x, E_y, E_z], where each - component E_* is an ndarray of equal shape. +The functions generated here expect field inputs with shape (3, X, Y, Z), +e.g. E = [E_x, E_y, E_z] where each component has shape (X, Y, Z) """ from typing import List, Callable import numpy @@ -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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Function for taking the discretized curl of the H-field, F(H) -> curlH """ dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij') @@ -41,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Function for taking the discretized curl of the E-field, F(E) -> curlE """ dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij') @@ -69,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Dielectric constant :param mu: Magnetic permeability (default 1 everywhere) :return: Function implementing the wave operator A(E) -> E @@ -100,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Dielectric constant :param mu: Magnetic permeability (default 1 everywhere) :return: Function implementing the wave operator A(E, H) -> (E, H) @@ -131,7 +131,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param mu: Magnetic permeability (default 1 everywhere) :return: Function for converting E to H """ @@ -159,7 +159,7 @@ def m2j(omega: complex, 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 fdfd_tools.operators header + :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 """ diff --git a/fdfd_tools/operators.py b/meanas/fdfd/operators.py similarity index 92% rename from fdfd_tools/operators.py rename to meanas/fdfd/operators.py index 8809d09..3b2de68 100644 --- a/fdfd_tools/operators.py +++ b/meanas/fdfd/operators.py @@ -3,17 +3,13 @@ Sparse matrix operators for use with electromagnetic wave equations. These functions return sparse-matrix (scipy.sparse.spmatrix) representations of a variety of operators, intended for use with E and H fields vectorized using the - fdfd_tools.vec() and .unvec() functions (column-major/Fortran ordering). + meanas.vec() and .unvec() functions (column-major/Fortran ordering). E- and H-field values are defined on a Yee cell; epsilon values should be calculated for cells centered at each E component (mu at each H component). -Many of these functions require a 'dxes' parameter, of type 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. +Many of these functions require a 'dxes' parameter, of type meanas.dx_lists_type; see +the meanas.types submodule for details. The following operators are included: @@ -57,7 +53,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Vectorized dielectric constant :param mu: Vectorized magnetic permeability (default 1 everywhere). :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted @@ -101,7 +97,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Preconditioner matrices (Pl, Pr) """ p_squared = [dxes[0][0][:, None, None] * dxes[1][1][None, :, None] * dxes[1][2][None, None, :], @@ -127,7 +123,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Vectorized dielectric constant :param mu: Vectorized magnetic permeability (default 1 everywhere) :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted @@ -177,7 +173,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Vectorized dielectric constant :param mu: Vectorized magnetic permeability (default 1 everywhere) :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted @@ -216,7 +212,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Sparse matrix for taking the discretized curl of the H-field """ return cross(deriv_back(dxes[1])) @@ -226,7 +222,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Sparse matrix for taking the discretized curl of the E-field """ return cross(deriv_forward(dxes[0])) @@ -242,7 +238,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param mu: Vectorized magnetic permeability (default 1 everywhere) :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). @@ -270,7 +266,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param mu: Vectorized magnetic permeability (default 1 everywhere) :return: Sparse matrix for converting E to H """ @@ -454,7 +450,7 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: Operator for computing the Poynting vector, containing the (E x) portion of the Poynting vector. :param e: Vectorized E-field for the ExH cross product - :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Sparse matrix containing (E x) portion of Poynting cross product """ shape = [len(dx) for dx in dxes[0]] @@ -483,7 +479,7 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector. :param h: Vectorized H-field for the HxE cross product - :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :return: Sparse matrix containing (H x) portion of Poynting cross product """ shape = [len(dx) for dx in dxes[0]] diff --git a/fdfd_tools/grid.py b/meanas/fdfd/scpml.py similarity index 99% rename from fdfd_tools/grid.py rename to meanas/fdfd/scpml.py index 8ecb44f..c4091a0 100644 --- a/fdfd_tools/grid.py +++ b/meanas/fdfd/scpml.py @@ -8,7 +8,6 @@ import numpy __author__ = 'Jan Petykiewicz' -dx_lists_t = List[List[numpy.ndarray]] s_function_type = Callable[[float], float] diff --git a/fdfd_tools/solvers.py b/meanas/fdfd/solvers.py similarity index 97% rename from fdfd_tools/solvers.py rename to meanas/fdfd/solvers.py index 066725c..a0ce403 100644 --- a/fdfd_tools/solvers.py +++ b/meanas/fdfd/solvers.py @@ -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 fdfd_tools.vec(). + All ndarray arguments should be 1D array, as returned by meanas.vec(). :param omega: Complex frequency to solve at. :param dxes: [[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]] (complex cell sizes) diff --git a/fdfd_tools/waveguide.py b/meanas/fdfd/waveguide.py similarity index 92% rename from fdfd_tools/waveguide.py rename to meanas/fdfd/waveguide.py index 89e0d9c..48a1510 100644 --- a/fdfd_tools/waveguide.py +++ b/meanas/fdfd/waveguide.py @@ -51,7 +51,7 @@ def operator(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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Sparse matrix representation of the operator @@ -91,7 +91,7 @@ def normalized_fields(v: numpy.ndarray, :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 dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Normalized, vectorized (e, h) containing all vector components. @@ -120,6 +120,8 @@ def normalized_fields(v: numpy.ndarray, # Try to break symmetry to assign a consistent sign [experimental] E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) + logger.debug('norm_angle = {}'.format(norm_angle)) + logger.debug('norm_sign = {}'.format(sign) norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle) @@ -140,7 +142,7 @@ def v2h(v: numpy.ndarray, :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 dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Vectorized H field with all vector components """ @@ -172,7 +174,7 @@ def v2e(v: numpy.ndarray, :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 dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Vectorized E field with all vector components. @@ -192,7 +194,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Sparse matrix representation of the operator """ @@ -213,7 +215,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :return: Sparse matrix representation of the operator """ @@ -226,7 +228,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :return: Sparse matrix representation of the operator """ n = 1 @@ -243,7 +245,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :return: Sparse matrix representation of the operator """ n = 1 @@ -268,7 +270,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Relative error norm(OP @ h) / norm(h) @@ -299,7 +301,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param mu: Vectorized magnetic permeability grid (default 1 everywhere) :return: Relative error norm(OP @ e) / norm(e) @@ -335,7 +337,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 fdfd_tools.operators header (2D) + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) :param epsilon: Vectorized dielectric constant grid :param r0: Radius of curvature for the simulation. This should be the minimum value of r within the simulation domain. diff --git a/fdfd_tools/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py similarity index 97% rename from fdfd_tools/waveguide_mode.py rename to meanas/fdfd/waveguide_mode.py index 0b839d8..6e2b192 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/meanas/fdfd/waveguide_mode.py @@ -19,7 +19,7 @@ def solve_waveguide_mode_2d(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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Dielectric constant :param mu: Magnetic permeability (default 1 everywhere) :param wavenumber_correction: Whether to correct the wavenumber to @@ -87,7 +87,7 @@ 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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param axis: Propagation axis (0=x, 1=y, 2=z) :param polarity: Propagation direction (+1 for +ve, -1 for -ve) :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use @@ -167,7 +167,7 @@ def compute_source(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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param axis: Propagation axis (0=x, 1=y, 2=z) :param polarity: Propagation direction (+1 for +ve, -1 for -ve) :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use @@ -219,7 +219,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 fdfd_tools.operators header + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param axis: Propagation axis (0=x, 1=y, 2=z) :param polarity: Propagation direction (+1 for +ve, -1 for -ve) :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use @@ -283,7 +283,7 @@ def solve_waveguide_mode_cylindrical(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 fdfd_tools.operators header. + :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types. The first coordinate is assumed to be r, the second is y. :param epsilon: Dielectric constant :param r0: Radius of curvature for the simulation. This should be the minimum value of diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py new file mode 100644 index 0000000..a1d278a --- /dev/null +++ b/meanas/fdtd/__init__.py @@ -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 diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py new file mode 100644 index 0000000..8dd1df3 --- /dev/null +++ b/meanas/fdtd/base.py @@ -0,0 +1,87 @@ +""" +Basic FDTD field updates +""" +from typing import List, Callable, Tuple, Dict +import numpy + +from .. import dx_lists_t, field_t, field_updater + +__author__ = 'Jan Petykiewicz' + + +def curl_h(dxes: dx_lists_t = None) -> field_updater: + """ + Curl operator for use with the H field. + + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :return: Function for taking the discretized curl of the H-field, F(H) -> curlH + """ + if dxes: + dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij') + + def dh(f, ax): + return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax] + else: + def dh(f, ax): + return f - numpy.roll(f, 1, axis=ax) + + def ch_fun(h: field_t) -> field_t: + output = numpy.empty_like(h) + output[0] = dh(h[2], 1) + output[1] = dh(h[0], 2) + output[2] = dh(h[1], 0) + output[0] -= dh(h[1], 2) + output[1] -= dh(h[2], 0) + output[2] -= dh(h[0], 1) + return output + + return ch_fun + + +def curl_e(dxes: dx_lists_t = None) -> field_updater: + """ + Curl operator for use with the E field. + + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :return: Function for taking the discretized curl of the E-field, F(E) -> curlE + """ + if dxes is not None: + dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij') + + def de(f, ax): + return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax] + else: + def de(f, ax): + return numpy.roll(f, -1, axis=ax) - f + + def ce_fun(e: field_t) -> field_t: + output = numpy.empty_like(e) + output[0] = de(e[2], 1) + output[1] = de(e[0], 2) + output[2] = de(e[1], 0) + output[0] -= de(e[1], 2) + output[1] -= de(e[2], 0) + output[2] -= de(e[0], 1) + return output + + return ce_fun + + +def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater: + curl_h_fun = curl_h(dxes) + + def me_fun(e: field_t, h: field_t, epsilon: field_t): + e += dt * curl_h_fun(h) / epsilon + return e + + return me_fun + + +def maxwell_h(dt: float, dxes: dx_lists_t = None) -> field_updater: + curl_e_fun = curl_e(dxes) + + def mh_fun(e: field_t, h: field_t): + h -= dt * curl_e_fun(e) + return h + + return mh_fun diff --git a/meanas/fdtd/boundaries.py b/meanas/fdtd/boundaries.py new file mode 100644 index 0000000..34a8d4a --- /dev/null +++ b/meanas/fdtd/boundaries.py @@ -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)) + + diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py new file mode 100644 index 0000000..26fb036 --- /dev/null +++ b/meanas/fdtd/energy.py @@ -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 diff --git a/meanas/fdtd/pml.py b/meanas/fdtd/pml.py new file mode 100644 index 0000000..3e10aa6 --- /dev/null +++ b/meanas/fdtd/pml.py @@ -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 diff --git a/fdfd_tools/test/test_fdtd.py b/meanas/test/test_fdtd.py similarity index 99% rename from fdfd_tools/test/test_fdtd.py rename to meanas/test/test_fdtd.py index d74b8be..8443afc 100644 --- a/fdfd_tools/test/test_fdtd.py +++ b/meanas/test/test_fdtd.py @@ -4,7 +4,7 @@ import dataclasses from typing import List, Tuple from numpy.testing import assert_allclose, assert_array_equal -from fdfd_tools import fdtd +from meanas import fdtd prng = numpy.random.RandomState(12345) diff --git a/meanas/types.py b/meanas/types.py new file mode 100644 index 0000000..7dc5c1c --- /dev/null +++ b/meanas/types.py @@ -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] diff --git a/fdfd_tools/vectorization.py b/meanas/vectorization.py similarity index 87% rename from fdfd_tools/vectorization.py rename to meanas/vectorization.py index 57b58fb..fd6cdcf 100644 --- a/fdfd_tools/vectorization.py +++ b/meanas/vectorization.py @@ -4,15 +4,13 @@ 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 -__author__ = 'Jan Petykiewicz' +from .types import field_t, vfield_t -# Types -field_t = List[numpy.ndarray] # vector field (eg. [E_x, E_y, E_z] -vfield_t = numpy.ndarray # linearized vector field + +__author__ = 'Jan Petykiewicz' def vec(f: field_t) -> vfield_t: @@ -27,7 +25,7 @@ def vec(f: field_t) -> vfield_t: """ if numpy.any(numpy.equal(f, None)): return None - return numpy.hstack(tuple((fi.ravel(order='C') for fi in f))) + return numpy.ravel(f, order='C') def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t: diff --git a/setup.py b/setup.py index 2a64b37..8e817cb 100644 --- a/setup.py +++ b/setup.py @@ -1,14 +1,14 @@ #!/usr/bin/env python3 from setuptools import setup, find_packages -import fdfd_tools +import meanas with open('README.md', 'r') as f: long_description = f.read() -setup(name='fdfd_tools', - version=fdfd_tools.version, - description='FDFD Electromagnetic simulation tools', +setup(name='meanas', + version=meanas.version, + description='Electromagnetic simulation tools', long_description=long_description, long_description_content_type='text/markdown', author='Jan Petykiewicz',