From a956323b94f1cff5d896dd9156af679c6626dd90 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Wed, 27 Nov 2019 22:59:52 -0800 Subject: [PATCH] move stuff under fdmath --- examples/fdfd.py | 16 ++-- meanas/__init__.py | 3 - meanas/eigensolvers.py | 15 ++-- meanas/fdfd/bloch.py | 28 +++--- meanas/fdfd/farfield.py | 10 +-- meanas/fdfd/functional.py | 52 ++++++----- meanas/fdfd/operators.py | 70 +++++++-------- meanas/fdfd/scpml.py | 10 +-- meanas/fdfd/solvers.py | 17 ++-- meanas/fdfd/waveguide_2d.py | 126 ++++++++++++++------------- meanas/fdfd/waveguide_3d.py | 32 +++---- meanas/fdfd/waveguide_cyl.py | 12 +-- meanas/fdmath/__init__.py | 8 ++ meanas/fdmath/functional.py | 14 +-- meanas/fdmath/operators.py | 4 +- meanas/{ => fdmath}/types.py | 6 +- meanas/{ => fdmath}/vectorization.py | 6 +- meanas/fdtd/base.py | 20 +++-- meanas/fdtd/boundaries.py | 12 +-- meanas/fdtd/energy.py | 83 +++++++++--------- meanas/fdtd/pml.py | 12 +-- meanas/test/test_fdfd.py | 17 +++- meanas/test/test_fdfd_pml.py | 3 +- 23 files changed, 304 insertions(+), 272 deletions(-) rename meanas/{ => fdmath}/types.py (89%) rename meanas/{ => fdmath}/vectorization.py (90%) diff --git a/examples/fdfd.py b/examples/fdfd.py index aa684e8..3fe895f 100644 --- a/examples/fdfd.py +++ b/examples/fdfd.py @@ -3,14 +3,18 @@ import numpy from numpy.linalg import norm import meanas -from meanas import vec, unvec, fdtd -from meanas.fdfd import waveguide_mode, functional, scpml, operators +from meanas import fdtd +from meanas.fdmath import vec, unvec +from meanas.fdfd import waveguide_3d, functional, scpml, operators from meanas.fdfd.solvers import generic as generic_solver import gridlock from matplotlib import pyplot +import logging + +logging.basicConfig(level=logging.DEBUG) __author__ = 'Jan Petykiewicz' @@ -134,10 +138,10 @@ def test1(solver=generic_solver): 'polarity': +1, } - wg_results = waveguide_mode.solve_waveguide_mode(mode_number=0, omega=omega, epsilon=grid.grids, **wg_args) - J = waveguide_mode.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'], - omega=omega, epsilon=grid.grids, **wg_args) - e_overlap = waveguide_mode.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args) + wg_results = waveguide_3d.solve_mode(mode_number=0, omega=omega, epsilon=grid.grids, **wg_args) + J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'], + omega=omega, epsilon=grid.grids, **wg_args) + e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args) pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) # pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) diff --git a/meanas/__init__.py b/meanas/__init__.py index 8d30d47..770bdf4 100644 --- a/meanas/__init__.py +++ b/meanas/__init__.py @@ -6,9 +6,6 @@ See the readme or `import meanas; help(meanas)` for more info. 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: diff --git a/meanas/eigensolvers.py b/meanas/eigensolvers.py index 960df60..51aeb6e 100644 --- a/meanas/eigensolvers.py +++ b/meanas/eigensolvers.py @@ -37,20 +37,17 @@ def power_iteration(operator: sparse.spmatrix, def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperator, - guess_vectors: numpy.ndarray, + guess_vector: numpy.ndarray, iterations: int = 40, tolerance: float = 1e-13, - solver=None, + solver = None, ) -> Tuple[complex, numpy.ndarray]: """ Use Rayleigh quotient iteration to refine an eigenvector guess. - TODO: - Need to test this for more than one guess_vectors. - Args: operator: Matrix to analyze. - guess_vectors: Eigenvectors to refine. + guess_vector: Eigenvector to refine. iterations: Maximum number of iterations to perform. Default 40. tolerance: Stop iteration if `(A - I*eigenvalue) @ v < num_vectors * tolerance`, Default 1e-13. @@ -73,16 +70,16 @@ def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperato if solver is None: solver = lambda A, b: spalg.bicgstab(A, b)[0] - v = numpy.atleast_2d(guess_vectors) + v = numpy.squeeze(guess_vector) v /= norm(v) for _ in range(iterations): eigval = v.conj() @ (operator @ v) - if norm(operator @ v - eigval * v) < v.shape[1] * tolerance: + if norm(operator @ v - eigval * v) < tolerance: break shifted_operator = operator - shift(eigval) v = solver(shifted_operator, v) - v /= norm(v, axis=0) + v /= norm(v) return eigval, v diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py index 0ca64a7..28b82df 100644 --- a/meanas/fdfd/bloch.py +++ b/meanas/fdfd/bloch.py @@ -83,7 +83,7 @@ import scipy.optimize from scipy.linalg import norm import scipy.sparse.linalg as spalg -from .. import field_t +from ..fdmath import fdfield_t logger = logging.getLogger(__name__) @@ -154,8 +154,8 @@ def generate_kmn(k0: numpy.ndarray, def maxwell_operator(k0: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t, - mu: field_t = None + epsilon: fdfield_t, + mu: fdfield_t = None ) -> Callable[[numpy.ndarray], numpy.ndarray]: """ Generate the Maxwell operator @@ -227,8 +227,8 @@ def maxwell_operator(k0: numpy.ndarray, def hmn_2_exyz(k0: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t, - ) -> Callable[[numpy.ndarray], field_t]: + epsilon: fdfield_t, + ) -> Callable[[numpy.ndarray], fdfield_t]: """ Generate an operator which converts a vectorized spatial-frequency-space h_mn into an E-field distribution, i.e. @@ -249,7 +249,7 @@ def hmn_2_exyz(k0: numpy.ndarray, k_mag, m, n = generate_kmn(k0, G_matrix, shape) - def operator(h: numpy.ndarray) -> field_t: + def operator(h: numpy.ndarray) -> fdfield_t: hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] d_xyz = (n * hin_m - m * hin_n) * k_mag @@ -262,8 +262,8 @@ def hmn_2_exyz(k0: numpy.ndarray, def hmn_2_hxyz(k0: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t - ) -> Callable[[numpy.ndarray], field_t]: + epsilon: fdfield_t + ) -> Callable[[numpy.ndarray], fdfield_t]: """ Generate an operator which converts a vectorized spatial-frequency-space h_mn into an H-field distribution, i.e. @@ -293,8 +293,8 @@ def hmn_2_hxyz(k0: numpy.ndarray, def inverse_maxwell_operator_approx(k0: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t, - mu: field_t = None + epsilon: fdfield_t, + mu: fdfield_t = None ) -> Callable[[numpy.ndarray], numpy.ndarray]: """ Generate an approximate inverse of the Maxwell operator, @@ -366,8 +366,8 @@ def find_k(frequency: float, tolerance: float, direction: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t, - mu: field_t = None, + epsilon: fdfield_t, + mu: fdfield_t = None, band: int = 0, k_min: float = 0, k_max: float = 0.5, @@ -409,8 +409,8 @@ def find_k(frequency: float, def eigsolve(num_modes: int, k0: numpy.ndarray, G_matrix: numpy.ndarray, - epsilon: field_t, - mu: field_t = None, + epsilon: fdfield_t, + mu: fdfield_t = None, tolerance: float = 1e-20, max_iters: int = 10000, reset_iters: int = 100, diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py index 665e70f..b4b9dbe 100644 --- a/meanas/fdfd/farfield.py +++ b/meanas/fdfd/farfield.py @@ -6,11 +6,11 @@ import numpy from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift from numpy import pi -from .. import field_t +from .. import fdfield_t -def near_to_farfield(E_near: field_t, - H_near: field_t, +def near_to_farfield(E_near: fdfield_t, + H_near: fdfield_t, dx: float, dy: float, padded_size: List[int] = None @@ -117,8 +117,8 @@ def near_to_farfield(E_near: field_t, -def far_to_nearfield(E_far: field_t, - H_far: field_t, +def far_to_nearfield(E_far: fdfield_t, + H_far: fdfield_t, dkx: float, dky: float, padded_size: List[int] = None diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py index 74b66da..d995e7b 100644 --- a/meanas/fdfd/functional.py +++ b/meanas/fdfd/functional.py @@ -2,32 +2,30 @@ Functional versions of many FDFD operators. These can be useful for performing FDFD calculations without needing to construct large matrices in memory. -The functions generated here expect `field_t` inputs with shape (3, X, Y, Z), +The functions generated here expect `fdfield_t` 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, Tuple import numpy -from .. import dx_lists_t, field_t +from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t from ..fdmath.functional import curl_forward, curl_back + __author__ = 'Jan Petykiewicz' -field_transform_t = Callable[[field_t], field_t] - - def e_full(omega: complex, dxes: dx_lists_t, - epsilon: field_t, - mu: field_t = None - ) -> field_transform_t: + epsilon: fdfield_t, + mu: fdfield_t = None + ) -> fdfield_updater_t: """ Wave operator for use with E-field. See `operators.e_full` for details. Args: omega: Angular frequency of the simulation - dxes: Grid parameters [dx_e, dx_h] as described in meanas.types + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Dielectric constant mu: Magnetic permeability (default 1 everywhere) @@ -54,16 +52,16 @@ def e_full(omega: complex, def eh_full(omega: complex, dxes: dx_lists_t, - epsilon: field_t, - mu: field_t = None - ) -> Callable[[field_t, field_t], Tuple[field_t, field_t]]: + epsilon: fdfield_t, + mu: fdfield_t = None + ) -> Callable[[fdfield_t, fdfield_t], Tuple[fdfield_t, fdfield_t]]: """ Wave operator for full (both E and H) field representation. See `operators.eh_full`. Args: omega: Angular frequency of the simulation - dxes: Grid parameters [dx_e, dx_h] as described in meanas.types + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Dielectric constant mu: Magnetic permeability (default 1 everywhere) @@ -90,15 +88,15 @@ def eh_full(omega: complex, def e2h(omega: complex, dxes: dx_lists_t, - mu: field_t = None, - ) -> field_transform_t: + mu: fdfield_t = None, + ) -> fdfield_updater_t: """ Utility operator for converting the `E` field into the `H` field. For use with `e_full` -- assumes that there is no magnetic current `M`. Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` mu: Magnetic permeability (default 1 everywhere) Return: @@ -121,8 +119,8 @@ def e2h(omega: complex, def m2j(omega: complex, dxes: dx_lists_t, - mu: field_t = None, - ) -> field_transform_t: + mu: fdfield_t = None, + ) -> fdfield_updater_t: """ Utility operator for converting magnetic current `M` distribution into equivalent electric current distribution `J`. @@ -130,7 +128,7 @@ def m2j(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` mu: Magnetic permeability (default 1 everywhere) Returns: @@ -153,12 +151,12 @@ def m2j(omega: complex, return m2j_mu -def e_tfsf_source(TF_region: field_t, +def e_tfsf_source(TF_region: fdfield_t, omega: complex, dxes: dx_lists_t, - epsilon: field_t, - mu: field_t = None, - ) -> field_transform_t: + epsilon: fdfield_t, + mu: fdfield_t = None, + ) -> fdfield_updater_t: """ Operator that turns an E-field distribution into a total-field/scattered-field (TFSF) source. @@ -168,7 +166,7 @@ def e_tfsf_source(TF_region: field_t, (i.e. in the scattered-field region). Should have the same shape as the simulation grid, e.g. `epsilon[0].shape`. omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Dielectric constant distribution mu: Magnetic permeability (default 1 everywhere) @@ -184,7 +182,7 @@ def e_tfsf_source(TF_region: field_t, return neg_iwj / (-1j * omega) -def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[field_t, field_t], field_t]: +def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[fdfield_t, fdfield_t], fdfield_t]: """ Generates a function that takes the single-frequency `E` and `H` fields and calculates the cross product `E` x `H` = \\( E \\times H \\) as required @@ -201,12 +199,12 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[field_t, field_t], field_t instead. Args: - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: Function `f` that returns E x H as required for the poynting vector. """ - def exh(e: field_t, h: field_t): + def exh(e: fdfield_t, h: fdfield_t): s = numpy.empty_like(e) ex = e[0] * dxes[0][0][:, None, None] ey = e[1] * dxes[0][1][None, :, None] diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 20cf96a..b90ec67 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -9,7 +9,7 @@ E- and H-field values are defined on a Yee cell; `epsilon` values should be calc cells centered at each E component (`mu` at each H component). Many of these functions require a `dxes` parameter, of type `dx_lists_t`; see -the `meanas.types` submodule for details. +the `meanas.fdmath.types` submodule for details. The following operators are included: @@ -31,7 +31,7 @@ from typing import List, Tuple import numpy import scipy.sparse as sparse -from .. import vec, dx_lists_t, vfield_t +from ..fdmath import vec, dx_lists_t, vfdfield_t from ..fdmath.operators import shift_with_mirror, rotation, curl_forward, curl_back @@ -40,10 +40,10 @@ __author__ = 'Jan Petykiewicz' def e_full(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - pec: vfield_t = None, - pmc: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, + pec: vfdfield_t = None, + pmc: vfdfield_t = None, ) -> sparse.spmatrix: """ Wave operator @@ -60,7 +60,7 @@ def e_full(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere). pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted @@ -107,7 +107,7 @@ def e_full_preconditioners(dxes: dx_lists_t The preconditioner matrices are diagonal and complex, with `Pr = 1 / Pl` Args: - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: Preconditioner matrices `(Pl, Pr)`. @@ -124,10 +124,10 @@ def e_full_preconditioners(dxes: dx_lists_t def h_full(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - pec: vfield_t = None, - pmc: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, + pec: vfdfield_t = None, + pmc: vfdfield_t = None, ) -> sparse.spmatrix: """ Wave operator @@ -142,7 +142,7 @@ def h_full(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere) pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted @@ -180,10 +180,10 @@ def h_full(omega: complex, def eh_full(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - pec: vfield_t = None, - pmc: vfield_t = None + epsilon: vfdfield_t, + mu: vfdfield_t = None, + pec: vfdfield_t = None, + pmc: vfdfield_t = None ) -> sparse.spmatrix: """ Wave operator for `[E, H]` field representation. This operator implements Maxwell's @@ -210,7 +210,7 @@ def eh_full(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere) pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted @@ -249,8 +249,8 @@ def eh_full(omega: complex, def e2h(omega: complex, dxes: dx_lists_t, - mu: vfield_t = None, - pmc: vfield_t = None, + mu: vfdfield_t = None, + pmc: vfdfield_t = None, ) -> sparse.spmatrix: """ Utility operator for converting the E field into the H field. @@ -258,7 +258,7 @@ def e2h(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` mu: Vectorized magnetic permeability (default 1 everywhere) pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted as containing a perfect magnetic conductor (PMC). @@ -280,7 +280,7 @@ def e2h(omega: complex, def m2j(omega: complex, dxes: dx_lists_t, - mu: vfield_t = None + mu: vfdfield_t = None ) -> sparse.spmatrix: """ Operator for converting a magnetic current M into an electric current J. @@ -288,7 +288,7 @@ def m2j(omega: complex, Args: omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` mu: Vectorized magnetic permeability (default 1 everywhere) Returns: @@ -302,14 +302,14 @@ def m2j(omega: complex, return op -def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: +def poynting_e_cross(e: vfdfield_t, dxes: dx_lists_t) -> sparse.spmatrix: """ Operator for computing the Poynting vector, containing the (E x) portion of the Poynting vector. Args: e: Vectorized E-field for the ExH cross product - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: Sparse matrix containing (E x) portion of Poynting cross product. @@ -331,13 +331,13 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: return P -def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: +def poynting_h_cross(h: vfdfield_t, dxes: dx_lists_t) -> sparse.spmatrix: """ Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector. Args: h: Vectorized H-field for the HxE cross product - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` Returns: Sparse matrix containing (H x) portion of Poynting cross product. @@ -358,11 +358,11 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: return P -def e_tfsf_source(TF_region: vfield_t, +def e_tfsf_source(TF_region: vfdfield_t, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, ) -> sparse.spmatrix: """ Operator that turns a desired E-field distribution into a @@ -374,7 +374,7 @@ def e_tfsf_source(TF_region: vfield_t, TF_region: Mask, which is set to 1 inside the total-field region and 0 in the scattered-field region omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere). @@ -388,11 +388,11 @@ def e_tfsf_source(TF_region: vfield_t, return (A @ Q - Q @ A) / (-1j * omega) -def e_boundary_source(mask: vfield_t, +def e_boundary_source(mask: vfdfield_t, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, periodic_mask_edges: bool = False, ) -> sparse.spmatrix: """ @@ -405,7 +405,7 @@ def e_boundary_source(mask: vfield_t, i.e. any points where shifting the mask by one cell in any direction would change its value. omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere). diff --git a/meanas/fdfd/scpml.py b/meanas/fdfd/scpml.py index c9a93e6..7e7694d 100644 --- a/meanas/fdfd/scpml.py +++ b/meanas/fdfd/scpml.py @@ -5,14 +5,14 @@ Functions for creating stretched coordinate perfectly matched layer (PML) absorb from typing import List, Callable import numpy -from .. import dx_lists_t +from ..fdmath import dx_lists_t __author__ = 'Jan Petykiewicz' s_function_t = Callable[[float], float] -"""Typedef for s-functions""" +"""Typedef for s-functions, see `prepare_s_function()`""" def prepare_s_function(ln_R: float = -16, @@ -63,7 +63,7 @@ def uniform_grid_scpml(shape: numpy.ndarray or List[int], Default uses `prepare_s_function()` with no parameters. Returns: - Complex cell widths (dx_lists_t) as discussed in `meanas.types`. + Complex cell widths (dx_lists_t) as discussed in `meanas.fdmath.types`. """ if s_function is None: s_function = prepare_s_function() @@ -102,7 +102,7 @@ def stretch_with_scpml(dxes: dx_lists_t, Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis. Args: - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` axis: axis to stretch (0=x, 1=y, 2=z) polarity: direction to stretch (-1 for -ve, +1 for +ve) omega: Angular frequency for the simulation @@ -113,7 +113,7 @@ def stretch_with_scpml(dxes: dx_lists_t, of pml parameters. Default uses `prepare_s_function()` with no parameters. Returns: - Complex cell widths (dx_lists_t) as discussed in `meanas.types`. + Complex cell widths (dx_lists_t) as discussed in `meanas.fdmath.types`. Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed. """ if s_function is None: diff --git a/meanas/fdfd/solvers.py b/meanas/fdfd/solvers.py index aa9633a..ff5bfe3 100644 --- a/meanas/fdfd/solvers.py +++ b/meanas/fdfd/solvers.py @@ -9,6 +9,7 @@ import numpy from numpy.linalg import norm import scipy.sparse.linalg +from ..fdmath import dx_lists_t, vfdfield_t from . import operators @@ -60,16 +61,16 @@ def _scipy_qmr(A: scipy.sparse.csr_matrix, def generic(omega: complex, - dxes: List[List[numpy.ndarray]], - J: numpy.ndarray, - epsilon: numpy.ndarray, - mu: numpy.ndarray = None, - pec: numpy.ndarray = None, - pmc: numpy.ndarray = None, + dxes: dx_lists_t, + J: vfdfield_t, + epsilon: vfdfield_t, + mu: vfdfield_t = None, + pec: vfdfield_t = None, + pmc: vfdfield_t = None, adjoint: bool = False, matrix_solver: Callable[..., numpy.ndarray] = _scipy_qmr, matrix_solver_opts: Dict[str, Any] = None, - ) -> numpy.ndarray: + ) -> vfdfield_t: """ Conjugate gradient FDFD solver using CSR sparse matrices. @@ -78,7 +79,7 @@ def generic(omega: complex, Args: omega: Complex frequency to solve at. dxes: `[[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]]` (complex cell sizes) as - discussed in `meanas.types` + discussed in `meanas.fdmath.types` J: Electric current distribution (at E-field locations) epsilon: Dielectric constant distribution (at E-field locations) mu: Magnetic permeability distribution (at H-field locations) diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index 6a7f24f..ecd40de 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -14,7 +14,8 @@ import numpy from numpy.linalg import norm import scipy.sparse as sparse -from .. import vec, unvec, dx_lists_t, field_t, vfield_t +from ..fdmath.operators import deriv_forward, deriv_back, curl_forward, curl_back, cross +from ..fdmath import vec, unvec, dx_lists_t, fdfield_t, vfdfield_t from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration from . import operators @@ -24,8 +25,8 @@ __author__ = 'Jan Petykiewicz' def operator_e(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, ) -> sparse.spmatrix: """ Waveguide operator of the form @@ -66,7 +67,7 @@ def operator_e(omega: complex, Args: omega: The angular frequency of the system. - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -76,8 +77,8 @@ def operator_e(omega: complex, if numpy.any(numpy.equal(mu, None)): mu = numpy.ones_like(epsilon) - Dfx, Dfy = operators.deriv_forward(dxes[0]) - Dbx, Dby = operators.deriv_back(dxes[1]) + Dfx, Dfy = deriv_forward(dxes[0]) + Dbx, Dby = deriv_back(dxes[1]) eps_parts = numpy.split(epsilon, 3) eps_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1]))) @@ -95,8 +96,8 @@ def operator_e(omega: complex, def operator_h(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, ) -> sparse.spmatrix: """ Waveguide operator of the form @@ -137,7 +138,7 @@ def operator_h(omega: complex, Args: omega: The angular frequency of the system. - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -169,10 +170,10 @@ def normalized_fields_e(e_xy: numpy.ndarray, wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: + ) -> Tuple[vfdfield_t, vfdfield_t]: """ Given a vector `e_xy` containing the vectorized E_x and E_y fields, returns normalized, vectorized E and H fields for the system. @@ -182,7 +183,7 @@ def normalized_fields_e(e_xy: numpy.ndarray, wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`. It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction. @@ -203,10 +204,10 @@ def normalized_fields_h(h_xy: numpy.ndarray, wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: + ) -> Tuple[vfdfield_t, vfdfield_t]: """ Given a vector `h_xy` containing the vectorized H_x and H_y fields, returns normalized, vectorized E and H fields for the system. @@ -216,7 +217,7 @@ def normalized_fields_h(h_xy: numpy.ndarray, wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`. It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction. @@ -237,10 +238,10 @@ def _normalized_fields(e: numpy.ndarray, h: numpy.ndarray, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: + ) -> Tuple[vfdfield_t, vfdfield_t]: # TODO documentation shape = [s.size for s in dxes[0]] dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)] @@ -276,8 +277,8 @@ def _normalized_fields(e: numpy.ndarray, def exy2h(wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None + epsilon: vfdfield_t, + mu: vfdfield_t = None ) -> sparse.spmatrix: """ Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, @@ -287,7 +288,7 @@ def exy2h(wavenumber: complex, wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`. It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -301,8 +302,8 @@ def exy2h(wavenumber: complex, def hxy2e(wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None + epsilon: vfdfield_t, + mu: vfdfield_t = None ) -> sparse.spmatrix: """ Operator which transforms the vector `h_xy` containing the vectorized H_x and H_y fields, @@ -312,7 +313,7 @@ def hxy2e(wavenumber: complex, wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`. It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -325,7 +326,7 @@ def hxy2e(wavenumber: complex, def hxy2h(wavenumber: complex, dxes: dx_lists_t, - mu: vfield_t = None + mu: vfdfield_t = None ) -> sparse.spmatrix: """ Operator which transforms the vector `h_xy` containing the vectorized H_x and H_y fields, @@ -334,13 +335,13 @@ def hxy2h(wavenumber: complex, Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`. It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy` - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) mu: Vectorized magnetic permeability grid (default 1 everywhere) Returns: Sparse matrix representing the operator. """ - Dfx, Dfy = operators.deriv_forward(dxes[0]) + Dfx, Dfy = deriv_forward(dxes[0]) hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber) if not numpy.any(numpy.equal(mu, None)): @@ -358,7 +359,7 @@ def hxy2h(wavenumber: complex, def exy2e(wavenumber: complex, dxes: dx_lists_t, - epsilon: vfield_t, + epsilon: vfdfield_t, ) -> sparse.spmatrix: """ Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, @@ -367,13 +368,13 @@ def exy2e(wavenumber: complex, Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy` - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid Returns: Sparse matrix representing the operator. """ - Dbx, Dby = operators.deriv_back(dxes[1]) + Dbx, Dby = deriv_back(dxes[1]) exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber) if not numpy.any(numpy.equal(epsilon, None)): @@ -392,7 +393,7 @@ def exy2e(wavenumber: complex, def e2h(wavenumber: complex, omega: complex, dxes: dx_lists_t, - mu: vfield_t = None + mu: vfdfield_t = None ) -> sparse.spmatrix: """ Returns an operator which, when applied to a vectorized E eigenfield, produces @@ -401,7 +402,7 @@ def e2h(wavenumber: complex, Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) mu: Vectorized magnetic permeability grid (default 1 everywhere) Returns: @@ -416,7 +417,7 @@ def e2h(wavenumber: complex, def h2e(wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t + epsilon: vfdfield_t ) -> sparse.spmatrix: """ Returns an operator which, when applied to a vectorized H eigenfield, produces @@ -425,7 +426,7 @@ def h2e(wavenumber: complex, Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid Returns: @@ -441,7 +442,7 @@ def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) Return: Sparse matrix representation of the operator. @@ -450,9 +451,10 @@ def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: for d in dxes[0]: n *= len(d) + print(wavenumber, n) Bz = -1j * wavenumber * sparse.eye(n) - Dfx, Dfy = operators.deriv_forward(dxes[0]) - return operators.cross([Dfx, Dfy, Bz]) + Dfx, Dfy = deriv_forward(dxes[0]) + return cross([Dfx, Dfy, Bz]) def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: @@ -461,7 +463,7 @@ def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: Args: wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) Return: Sparse matrix representation of the operator. @@ -471,16 +473,16 @@ def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: n *= len(d) Bz = -1j * wavenumber * sparse.eye(n) - Dbx, Dby = operators.deriv_back(dxes[1]) - return operators.cross([Dbx, Dby, Bz]) + Dbx, Dby = deriv_back(dxes[1]) + return cross([Dbx, Dby, Bz]) -def h_err(h: vfield_t, +def h_err(h: vfdfield_t, wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None + epsilon: vfdfield_t, + mu: vfdfield_t = None ) -> float: """ Calculates the relative error in the H field @@ -489,7 +491,7 @@ def h_err(h: vfield_t, h: Vectorized H field wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -509,12 +511,12 @@ def h_err(h: vfield_t, return norm(op) / norm(h) -def e_err(e: vfield_t, +def e_err(e: vfdfield_t, wavenumber: complex, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None + epsilon: vfdfield_t, + mu: vfdfield_t = None ) -> float: """ Calculates the relative error in the E field @@ -523,7 +525,7 @@ def e_err(e: vfield_t, e: Vectorized E field wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)` omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid mu: Vectorized magnetic permeability grid (default 1 everywhere) @@ -545,17 +547,17 @@ def e_err(e: vfield_t, def solve_modes(mode_numbers: List[int], omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, + epsilon: vfdfield_t, + mu: vfdfield_t = None, mode_margin: int = 2, - ) -> Tuple[List[vfield_t], List[complex]]: + ) -> Tuple[List[vfdfield_t], List[complex]]: """ Given a 2D region, attempts to solve for the eigenmode with the specified mode numbers. Args: mode_numbers: List of 0-indexed mode numbers to solve for omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` epsilon: Dielectric constant mu: Magnetic permeability (default 1 everywhere) mode_margin: The eigensolver will actually solve for `(max(mode_number) + mode_margin)` @@ -570,17 +572,18 @@ def solve_modes(mode_numbers: List[int], Solve for the largest-magnitude eigenvalue of the real operator ''' dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] - A_r = waveguide.operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu)) + A_r = operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu)) - eigvals, eigvecs = signed_eigensolve(A_r, max(mode_number) + mode_margin) - e_xys = eigvecs[:, -(numpy.array(mode_number) + 1)] + eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin) + e_xys = eigvecs[:, -(numpy.array(mode_numbers) + 1)] ''' Now solve for the eigenvector of the full operator, using the real operator's eigenvector as an initial guess for Rayleigh quotient iteration. ''' - A = waveguide.operator_e(omega, dxes, epsilon, mu) - eigvals, e_xys = rayleigh_quotient_iteration(A, e_xys) + A = operator_e(omega, dxes, epsilon, mu) + for nn in range(len(mode_numbers)): + eigvals[nn], e_xys[:, nn] = rayleigh_quotient_iteration(A, e_xys[:, nn]) # Calculate the wave-vector (force the real part to be positive) wavenumbers = numpy.sqrt(eigvals) @@ -592,7 +595,7 @@ def solve_modes(mode_numbers: List[int], def solve_mode(mode_number: int, *args, **kwargs - ) -> Tuple[vfield_t, complex]: + ) -> Tuple[vfdfield_t, complex]: """ Wrapper around `solve_modes()` that solves for a single mode. @@ -604,4 +607,5 @@ def solve_mode(mode_number: int, Returns: (e_xy, wavenumber) """ - return solve_modes(mode_numbers=[mode_number], *args, **kwargs) + e_xys, wavenumbers = solve_modes(mode_numbers=[mode_number], *args, **kwargs) + return e_xys[:, 0], wavenumbers[0] diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index 18727ce..02fb7fd 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -8,7 +8,7 @@ from typing import Dict, List, Tuple import numpy import scipy.sparse as sparse -from .. import vec, unvec, dx_lists_t, vfield_t, field_t +from ..fdmath import vec, unvec, dx_lists_t, vfdfield_t, fdfield_t from . import operators, waveguide_2d, functional @@ -18,8 +18,8 @@ def solve_mode(mode_number: int, axis: int, polarity: int, slices: List[slice], - epsilon: field_t, - mu: field_t = None, + epsilon: fdfield_t, + mu: fdfield_t = None, ) -> Dict[str, complex or numpy.ndarray]: """ Given a 3D grid, selects a slice from the grid and attempts to @@ -28,7 +28,7 @@ def solve_mode(mode_number: int, Args: mode_number: Number of the mode, 0-indexed omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` axis: Propagation axis (0=x, 1=y, 2=z) polarity: Propagation direction (+1 for +ve, -1 for -ve) slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use @@ -71,7 +71,7 @@ def solve_mode(mode_number: int, wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2) shape = [d.size for d in args_2d['dxes'][0]] - ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber) + ve, vh = waveguide_2d.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber) e = unvec(ve, shape) h = unvec(vh, shape) @@ -98,16 +98,16 @@ def solve_mode(mode_number: int, return results -def compute_source(E: field_t, +def compute_source(E: fdfield_t, wavenumber: complex, omega: complex, dxes: dx_lists_t, axis: int, polarity: int, slices: List[slice], - epsilon: field_t, - mu: field_t = None, - ) -> field_t: + epsilon: fdfield_t, + mu: fdfield_t = None, + ) -> fdfield_t: """ Given an eigenmode obtained by `solve_mode`, returns the current source distribution necessary to position a unidirectional source at the slice location. @@ -116,7 +116,7 @@ def compute_source(E: field_t, E: E-field of the mode wavenumber: Wavenumber of the mode omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` axis: Propagation axis (0=x, 1=y, 2=z) polarity: Propagation direction (+1 for +ve, -1 for -ve) slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use @@ -143,13 +143,13 @@ def compute_source(E: field_t, return J -def compute_overlap_e(E: field_t, +def compute_overlap_e(E: fdfield_t, wavenumber: complex, dxes: dx_lists_t, axis: int, polarity: int, slices: List[slice], - ) -> field_t: # TODO DOCS + ) -> fdfield_t: # TODO DOCS """ Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn) @@ -160,7 +160,7 @@ def compute_overlap_e(E: field_t, H: H-field of the mode (advanced by half of a Yee cell from E) wavenumber: Wavenumber of the mode omega: Angular frequency of the simulation - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` axis: Propagation axis (0=x, 1=y, 2=z) polarity: Propagation direction (+1 for +ve, -1 for -ve) slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use @@ -188,13 +188,13 @@ def compute_overlap_e(E: field_t, return Etgt -def expand_e(E: field_t, +def expand_e(E: fdfield_t, wavenumber: complex, dxes: dx_lists_t, axis: int, polarity: int, slices: List[slice], - ) -> field_t: + ) -> fdfield_t: """ Given an eigenmode obtained by `solve_mode`, expands the E-field from the 2D slice where the mode was calculated to the entire domain (along the propagation @@ -205,7 +205,7 @@ def expand_e(E: field_t, Args: E: E-field of the mode wavenumber: Wavenumber of the mode - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` axis: Propagation axis (0=x, 1=y, 2=z) polarity: Propagation direction (+1 for +ve, -1 for -ve) slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index ebfb41d..0995984 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -13,7 +13,7 @@ import numpy from numpy.linalg import norm import scipy.sparse as sparse -from .. import vec, unvec, dx_lists_t, field_t, vfield_t +from ..fdmath import vec, unvec, dx_lists_t, fdfield_t, vfdfield_t from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration from . import operators @@ -23,7 +23,7 @@ __author__ = 'Jan Petykiewicz' def cylindrical_operator(omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, + epsilon: vfdfield_t, r0: float, ) -> sparse.spmatrix: """ @@ -41,7 +41,7 @@ def cylindrical_operator(omega: complex, Args: omega: The angular frequency of the system - dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) epsilon: Vectorized dielectric constant grid r0: Radius of curvature for the simulation. This should be the minimum value of r within the simulation domain. @@ -83,9 +83,9 @@ def cylindrical_operator(omega: complex, def solve_mode(mode_number: int, omega: complex, dxes: dx_lists_t, - epsilon: vfield_t, + epsilon: vfdfield_t, r0: float, - ) -> Dict[str, complex or field_t]: + ) -> Dict[str, complex or fdfield_t]: """ TODO: fixup Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode @@ -94,7 +94,7 @@ def solve_mode(mode_number: int, Args: mode_number: Number of the mode, 0-indexed omega: Angular frequency of the simulation - dxes: Grid parameters [dx_e, dx_h] as described in meanas.types. + dxes: Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types. The first coordinate is assumed to be r, the second is y. epsilon: Dielectric constant r0: Radius of curvature for the simulation. This should be the minimum value of diff --git a/meanas/fdmath/__init__.py b/meanas/fdmath/__init__.py index 74deff1..f70634f 100644 --- a/meanas/fdmath/__init__.py +++ b/meanas/fdmath/__init__.py @@ -91,4 +91,12 @@ and while \\( \\hat{h} \\) and \\( \\tilde{h} \\) are located at \\((m \pm \\frac{1}{2}, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})\\) with components at \\((m, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})\\) etc. +TODO: Explain fdfield_t vs vfdfield_t / operators vs functional +TODO: explain dxes + """ + +from .types import fdfield_t, vfdfield_t, dx_lists_t, fdfield_updater_t +from .vectorization import vec, unvec +from . import operators, functional, types, vectorization + diff --git a/meanas/fdmath/functional.py b/meanas/fdmath/functional.py index 8593d4a..fc5b0ca 100644 --- a/meanas/fdmath/functional.py +++ b/meanas/fdmath/functional.py @@ -6,10 +6,10 @@ Basic discrete calculus etc. from typing import List, Callable, Tuple, Dict import numpy -from .. import field_t, field_updater +from .types import fdfield_t, fdfield_updater_t -def deriv_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: +def deriv_forward(dx_e: List[numpy.ndarray] = None) -> fdfield_updater_t: """ Utility operators for taking discretized derivatives (backward variant). @@ -31,7 +31,7 @@ def deriv_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: return derivs -def deriv_back(dx_h: List[numpy.ndarray] = None) -> field_updater: +def deriv_back(dx_h: List[numpy.ndarray] = None) -> fdfield_updater_t: """ Utility operators for taking discretized derivatives (forward variant). @@ -53,7 +53,7 @@ def deriv_back(dx_h: List[numpy.ndarray] = None) -> field_updater: return derivs -def curl_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: +def curl_forward(dx_e: List[numpy.ndarray] = None) -> fdfield_updater_t: """ Curl operator for use with the E field. @@ -67,7 +67,7 @@ def curl_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: """ Dx, Dy, Dz = deriv_forward(dx_e) - def ce_fun(e: field_t) -> field_t: + def ce_fun(e: fdfield_t) -> fdfield_t: output = numpy.empty_like(e) output[0] = Dy(e[2]) output[1] = Dz(e[0]) @@ -80,7 +80,7 @@ def curl_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: return ce_fun -def curl_back(dx_h: List[numpy.ndarray] = None) -> field_updater: +def curl_back(dx_h: List[numpy.ndarray] = None) -> fdfield_updater_t: """ Create a function which takes the backward curl of a field. @@ -94,7 +94,7 @@ def curl_back(dx_h: List[numpy.ndarray] = None) -> field_updater: """ Dx, Dy, Dz = deriv_back(dx_h) - def ch_fun(h: field_t) -> field_t: + def ch_fun(h: fdfield_t) -> fdfield_t: output = numpy.empty_like(h) output[0] = Dy(h[2]) output[1] = Dz(h[0]) diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index 64f04aa..3f83c8c 100644 --- a/meanas/fdmath/operators.py +++ b/meanas/fdmath/operators.py @@ -7,7 +7,7 @@ from typing import List, Callable, Tuple, Dict import numpy import scipy.sparse as sparse -from .. import field_t, vfield_t +from .types import fdfield_t, vfdfield_t def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: @@ -155,7 +155,7 @@ def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix: [-B[1], B[0], zero]]) -def vec_cross(b: vfield_t) -> sparse.spmatrix: +def vec_cross(b: vfdfield_t) -> sparse.spmatrix: """ Vector cross product operator diff --git a/meanas/types.py b/meanas/fdmath/types.py similarity index 89% rename from meanas/types.py rename to meanas/fdmath/types.py index 2a1351f..6e2fd63 100644 --- a/meanas/types.py +++ b/meanas/fdmath/types.py @@ -8,7 +8,7 @@ from typing import List, Callable # Field types # TODO: figure out a better way to set the docstrings without creating actual subclasses? # Probably not a big issue since they're only used for type hinting -class field_t(numpy.ndarray): +class fdfield_t(numpy.ndarray): """ Vector field with shape (3, X, Y, Z) (e.g. `[E_x, E_y, E_z]`) @@ -16,7 +16,7 @@ class field_t(numpy.ndarray): """ pass -class vfield_t(numpy.ndarray): +class vfdfield_t(numpy.ndarray): """ Linearized vector field (single vector of length 3*X*Y*Z) @@ -37,4 +37,4 @@ dx_lists_t = List[List[numpy.ndarray]] ''' -field_updater = Callable[[field_t], field_t] +fdfield_updater_t = Callable[[fdfield_t], fdfield_t] diff --git a/meanas/vectorization.py b/meanas/fdmath/vectorization.py similarity index 90% rename from meanas/vectorization.py rename to meanas/fdmath/vectorization.py index fd6cdcf..8e8099b 100644 --- a/meanas/vectorization.py +++ b/meanas/fdmath/vectorization.py @@ -7,13 +7,13 @@ Vectorized versions of the field use row-major (ie., C-style) ordering. from typing import List import numpy -from .types import field_t, vfield_t +from .types import fdfield_t, vfdfield_t __author__ = 'Jan Petykiewicz' -def vec(f: field_t) -> vfield_t: +def vec(f: fdfield_t) -> vfdfield_t: """ Create a 1D ndarray from a 3D vector field which spans a 1-3D region. @@ -28,7 +28,7 @@ def vec(f: field_t) -> vfield_t: return numpy.ravel(f, order='C') -def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t: +def unvec(v: vfdfield_t, shape: numpy.ndarray) -> fdfield_t: """ Perform the inverse of vec(): take a 1D ndarray and output a 3D field of form [f_x, f_y, f_z] where each of f_* is a len(shape)-dimensional diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index 87a234b..fa478ba 100644 --- a/meanas/fdtd/base.py +++ b/meanas/fdtd/base.py @@ -4,27 +4,33 @@ Basic FDTD field updates from typing import List, Callable, Tuple, Dict import numpy -from .. import dx_lists_t, field_t, field_updater +from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t from ..fdmath.functional import curl_forward, curl_back __author__ = 'Jan Petykiewicz' -def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater: - curl_h_fun = curl_back(dxes[1]) +def maxwell_e(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: + if dxes is not None: + curl_h_fun = curl_back(dxes[1]) + else: + curl_h_fun = curl_back() - def me_fun(e: field_t, h: field_t, epsilon: field_t): + def me_fun(e: fdfield_t, h: fdfield_t, epsilon: fdfield_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_forward(dxes[0]) +def maxwell_h(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t: + if dxes is not None: + curl_e_fun = curl_forward(dxes[0]) + else: + curl_e_fun = curl_forward() - def mh_fun(e: field_t, h: field_t): + def mh_fun(e: fdfield_t, h: fdfield_t): h -= dt * curl_e_fun(e) return h diff --git a/meanas/fdtd/boundaries.py b/meanas/fdtd/boundaries.py index cba1797..0f0b1a6 100644 --- a/meanas/fdtd/boundaries.py +++ b/meanas/fdtd/boundaries.py @@ -5,12 +5,12 @@ Boundary conditions from typing import List, Callable, Tuple, Dict import numpy -from .. import dx_lists_t, field_t, field_updater +from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t def conducting_boundary(direction: int, polarity: int - ) -> Tuple[field_updater, field_updater]: + ) -> Tuple[fdfield_updater_t, fdfield_updater_t]: dirs = [0, 1, 2] if direction not in dirs: raise Exception('Invalid direction: {}'.format(direction)) @@ -23,13 +23,13 @@ def conducting_boundary(direction: int, boundary_slice[direction] = 0 shifted1_slice[direction] = 1 - def en(e: field_t): + def en(e: fdfield_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): + def hn(h: fdfield_t): h[direction][boundary_slice] = h[direction][shifted1_slice] h[u][boundary_slice] = 0 h[v][boundary_slice] = 0 @@ -45,14 +45,14 @@ def conducting_boundary(direction: int, shifted1_slice[direction] = -2 shifted2_slice[direction] = -3 - def ep(e: field_t): + def ep(e: fdfield_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): + def hp(h: fdfield_t): h[direction][boundary_slice] = h[direction][shifted1_slice] h[u][boundary_slice] = -h[u][shifted2_slice] h[u][shifted1_slice] = 0 diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py index 608137f..41268b7 100644 --- a/meanas/fdtd/energy.py +++ b/meanas/fdtd/energy.py @@ -2,12 +2,13 @@ from typing import List, Callable, Tuple, Dict import numpy -from .. import dx_lists_t, field_t, field_updater, fdmath +from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t +from ..fdmath.functional import deriv_back, deriv_forward -def poynting(e: field_t, - h: field_t, +def poynting(e: fdfield_t, + h: fdfield_t, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: if dxes is None: dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) @@ -25,51 +26,51 @@ def poynting(e: field_t, return s -def poynting_divergence(s: field_t = None, +def poynting_divergence(s: fdfield_t = None, *, - e: field_t = None, - h: field_t = None, + e: fdfield_t = None, + h: fdfield_t = None, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: if s is None: s = poynting(e, h, dxes=dxes) - Dx, Dy, Dz = fdmath.functional.deriv_back() + Dx, Dy, Dz = deriv_back() ds = Dx(s[0]) + Dy(s[1]) + Dz(s[2]) return ds -def energy_hstep(e0: field_t, - h1: field_t, - e2: field_t, - epsilon: field_t = None, - mu: field_t = None, +def energy_hstep(e0: fdfield_t, + h1: fdfield_t, + e2: fdfield_t, + epsilon: fdfield_t = None, + mu: fdfield_t = None, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: u = dxmul(e0 * e2, h1 * h1, epsilon, mu, dxes) return u -def energy_estep(h0: field_t, - e1: field_t, - h2: field_t, - epsilon: field_t = None, - mu: field_t = None, +def energy_estep(h0: fdfield_t, + e1: fdfield_t, + h2: fdfield_t, + epsilon: fdfield_t = None, + mu: fdfield_t = None, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: u = dxmul(e1 * e1, h0 * h2, epsilon, mu, dxes) return u def delta_energy_h2e(dt: float, - e0: field_t, - h1: field_t, - e2: field_t, - h3: field_t, - epsilon: field_t = None, - mu: field_t = None, + e0: fdfield_t, + h1: fdfield_t, + e2: fdfield_t, + h3: fdfield_t, + epsilon: fdfield_t = None, + mu: fdfield_t = None, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: """ This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2) """ @@ -80,14 +81,14 @@ def delta_energy_h2e(dt: float, def delta_energy_e2h(dt: float, - h0: field_t, - e1: field_t, - h2: field_t, - e3: field_t, - epsilon: field_t = None, - mu: field_t = None, + h0: fdfield_t, + e1: fdfield_t, + h2: fdfield_t, + e3: fdfield_t, + epsilon: fdfield_t = None, + mu: fdfield_t = None, dxes: dx_lists_t = None, - ) -> field_t: + ) -> fdfield_t: """ This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2) """ @@ -97,7 +98,7 @@ def delta_energy_e2h(dt: float, return du -def delta_energy_j(j0: field_t, e1: field_t, dxes: dx_lists_t = None) -> field_t: +def delta_energy_j(j0: fdfield_t, e1: fdfield_t, dxes: dx_lists_t = None) -> fdfield_t: if dxes is None: dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2)) @@ -108,12 +109,12 @@ def delta_energy_j(j0: field_t, e1: field_t, dxes: dx_lists_t = None) -> field_t return du -def dxmul(ee: field_t, - hh: field_t, - epsilon: field_t = None, - mu: field_t = None, +def dxmul(ee: fdfield_t, + hh: fdfield_t, + epsilon: fdfield_t = None, + mu: fdfield_t = None, dxes: dx_lists_t = None - ) -> field_t: + ) -> fdfield_t: if epsilon is None: epsilon = 1 if mu is None: diff --git a/meanas/fdtd/pml.py b/meanas/fdtd/pml.py index 7d73e38..af15a5a 100644 --- a/meanas/fdtd/pml.py +++ b/meanas/fdtd/pml.py @@ -7,7 +7,7 @@ PML implementations from typing import List, Callable, Tuple, Dict import numpy -from .. import dx_lists_t, field_t, field_updater +from ..fdmath import dx_lists_t, fdfield_t, fdfield_updater_t __author__ = 'Jan Petykiewicz' @@ -16,7 +16,7 @@ __author__ = 'Jan Petykiewicz' def cpml(direction: int, polarity: int, dt: float, - epsilon: field_t, + epsilon: fdfield_t, thickness: int = 8, ln_R_per_layer: float = -1.6, epsilon_eff: float = 1, @@ -25,7 +25,7 @@ def cpml(direction: int, ma: float = 1, cfs_alpha: float = 0, dtype: numpy.dtype = numpy.float32, - ) -> Tuple[Callable, Callable, Dict[str, field_t]]: + ) -> Tuple[Callable, Callable, Dict[str, fdfield_t]]: if direction not in range(3): raise Exception('Invalid direction: {}'.format(direction)) @@ -58,6 +58,7 @@ def cpml(direction: int, expand_slice = [None] * 3 expand_slice[direction] = slice(None) + expand_slice = tuple(expand_slice) def par(x): scaling = (x / thickness) ** m @@ -79,6 +80,7 @@ def cpml(direction: int, region[direction] = slice(-thickness, None) else: raise Exception('Bad polarity!') + region = tuple(region) se = 1 if direction == 1 else -1 @@ -97,7 +99,7 @@ def cpml(direction: int, # 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]: + def pml_e(e: fdfield_t, h: fdfield_t, epsilon: fdfield_t) -> Tuple[fdfield_t, fdfield_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 @@ -108,7 +110,7 @@ def cpml(direction: int, 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]: + def pml_h(e: fdfield_t, h: fdfield_t) -> Tuple[fdfield_t, fdfield_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 diff --git a/meanas/test/test_fdfd.py b/meanas/test/test_fdfd.py index ea6f417..2cf0c44 100644 --- a/meanas/test/test_fdfd.py +++ b/meanas/test/test_fdfd.py @@ -5,7 +5,8 @@ import pytest import numpy #from numpy.testing import assert_allclose, assert_array_equal -from .. import fdfd, vec, unvec +from .. import fdfd +from ..fdmath import vec, unvec from .utils import assert_close, assert_fields_close @@ -20,6 +21,10 @@ def test_poynting_planes(sim): mask = (sim.j != 0).any(axis=0) if mask.sum() != 2: pytest.skip(f'test_poynting_planes will only test 2-point sources, got {mask.sum()}') +# for dxg in sim.dxes: +# for dxa in dxg: +# if not (dxa == sim.dxes[0][0][0]).all(): +# pytest.skip('test_poynting_planes skips nonuniform dxes') points = numpy.where(mask) mask[points[0][0], points[1][0], points[2][0]] = 0 @@ -43,7 +48,6 @@ def test_poynting_planes(sim): assert_close(sum(planes), src_energy.sum()) - ##################################### # Test fixtures ##################################### @@ -102,6 +106,15 @@ def sim(request, shape, epsilon, dxes, j_distribution, omega, pec, pmc): # if is3d: # pytest.skip('Skipping dt != 0.3 because test is 3D (for speed)') +# # If no edge currents, add pmls +# src_mask = j_distribution.any(axis=0) +# th = 10 +# #if src_mask.sum() - src_mask[th:-th, th:-th, th:-th].sum() == 0: +# if src_mask.sum() - src_mask[th:-th, :, :].sum() == 0: +# for axis in (0,): +# for polarity in (-1, 1): +# dxes = fdfd.scpml.stretch_with_scpml(dxes, axis=axis, polarity=polarity, + j_vec = vec(j_distribution) eps_vec = vec(epsilon) e_vec = fdfd.solvers.generic(J=j_vec, omega=omega, dxes=dxes, epsilon=eps_vec, diff --git a/meanas/test/test_fdfd_pml.py b/meanas/test/test_fdfd_pml.py index eafa14f..d13d753 100644 --- a/meanas/test/test_fdfd_pml.py +++ b/meanas/test/test_fdfd_pml.py @@ -6,7 +6,8 @@ import pytest import numpy from numpy.testing import assert_allclose, assert_array_equal -from .. import fdfd, vec, unvec +from .. import fdfd +from ..fdmath import vec, unvec from .utils import assert_close, assert_fields_close from .test_fdfd import FDResult