From d6e7e3dee166b17d1d28b03a82f87aa032199f5a Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 24 Nov 2019 23:47:31 -0800 Subject: [PATCH] Big documentation and structure updates - Split math into fdmath package - Rename waveguide into _2d _3d and _cyl variants - pdoc-based documentation --- make_docs.sh | 3 + meanas/eigensolvers.py | 26 +- meanas/fdfd/__init__.py | 17 +- meanas/fdfd/bloch.py | 2 +- meanas/fdfd/farfield.py | 6 +- meanas/fdfd/functional.py | 192 +++++------ meanas/fdfd/operators.py | 444 ++++++++++--------------- meanas/fdfd/scpml.py | 102 +++--- meanas/fdfd/solvers.py | 57 ++-- meanas/fdfd/waveguide.py | 492 --------------------------- meanas/fdfd/waveguide_2d.py | 607 ++++++++++++++++++++++++++++++++++ meanas/fdfd/waveguide_3d.py | 236 +++++++++++++ meanas/fdfd/waveguide_cyl.py | 138 ++++++++ meanas/fdfd/waveguide_mode.py | 307 ----------------- meanas/fdmath/functional.py | 109 ++++++ meanas/fdmath/operators.py | 231 +++++++++++++ meanas/fdtd/__init__.py | 2 +- meanas/fdtd/base.py | 64 +--- meanas/fdtd/energy.py | 1 + meanas/test/__init__.py | 3 + meanas/test/test_fdfd_pml.py | 5 +- meanas/types.py | 17 +- pdoc_templates/config.mako | 46 +++ pdoc_templates/css.mako | 389 ++++++++++++++++++++++ pdoc_templates/html.mako | 421 +++++++++++++++++++++++ 25 files changed, 2579 insertions(+), 1338 deletions(-) create mode 100755 make_docs.sh delete mode 100644 meanas/fdfd/waveguide.py create mode 100644 meanas/fdfd/waveguide_2d.py create mode 100644 meanas/fdfd/waveguide_3d.py create mode 100644 meanas/fdfd/waveguide_cyl.py delete mode 100644 meanas/fdfd/waveguide_mode.py create mode 100644 meanas/fdmath/functional.py create mode 100644 meanas/fdmath/operators.py create mode 100644 pdoc_templates/config.mako create mode 100644 pdoc_templates/css.mako create mode 100644 pdoc_templates/html.mako diff --git a/make_docs.sh b/make_docs.sh new file mode 100755 index 0000000..f2f2fe2 --- /dev/null +++ b/make_docs.sh @@ -0,0 +1,3 @@ +#!/bin/bash +cd ~/projects/meanas +pdoc3 --html --force --template-dir pdoc_templates -o doc . diff --git a/meanas/eigensolvers.py b/meanas/eigensolvers.py index 5ca9962..960df60 100644 --- a/meanas/eigensolvers.py +++ b/meanas/eigensolvers.py @@ -15,10 +15,13 @@ def power_iteration(operator: sparse.spmatrix, """ Use power iteration to estimate the dominant eigenvector of a matrix. - :param operator: Matrix to analyze. - :param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector. - :param iterations: Number of iterations to perform. Default 20. - :return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate) + Args: + operator: Matrix to analyze. + guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector. + iterations: Number of iterations to perform. Default 20. + + Returns: + (Largest-magnitude eigenvalue, Corresponding eigenvector estimate) """ if numpy.any(numpy.equal(guess_vector, None)): v = numpy.random.rand(operator.shape[0]) @@ -91,12 +94,15 @@ def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator, Find the largest-magnitude positive-only (or negative-only) eigenvalues and eigenvectors of the provided matrix. - :param operator: Matrix to analyze. - :param how_many: How many eigenvalues to find. - :param negative: Whether to find negative-only eigenvalues. - Default False (positive only). - :return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors) - eigenvectors[:, k] corresponds to the k-th eigenvalue + Args: + operator: Matrix to analyze. + how_many: How many eigenvalues to find. + negative: Whether to find negative-only eigenvalues. + Default False (positive only). + + Returns: + (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors) + `eigenvectors[:, k]` corresponds to the k-th eigenvalue """ # Use power iteration to estimate the dominant eigenvector lm_eigval, _ = power_iteration(operator) diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py index 80ba55e..d0d58a9 100644 --- a/meanas/fdfd/__init__.py +++ b/meanas/fdfd/__init__.py @@ -1,2 +1,17 @@ -from . import solvers, operators, functional, scpml, waveguide, waveguide_mode +""" +Tools for finite difference frequency-domain (FDFD) simulations and calculations. + +These mostly involve picking a single frequency, then setting up and solving a +matrix equation (Ax=b) or eigenvalue problem. + + +Submodules: + +- `operators`, `functional`: General FDFD problem setup. +- `solvers`: Solver interface and reference implementation. +- `scpml`: Stretched-coordinate perfectly matched layer (scpml) boundary conditions +- `waveguide_2d`: Operators and mode-solver for waveguides with constant cross-section. +- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D. +""" +from . import solvers, operators, functional, scpml, waveguide_2d, waveguide_3d # from . import farfield, bloch TODO diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py index 9c252e7..0ca64a7 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 .. import field_t logger = logging.getLogger(__name__) diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py index faa25b5..665e70f 100644 --- a/meanas/fdfd/farfield.py +++ b/meanas/fdfd/farfield.py @@ -1,7 +1,7 @@ """ Functions for performing near-to-farfield transformation (and the reverse). """ -from typing import Dict, List +from typing import Dict, List, Any import numpy from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift from numpy import pi @@ -14,7 +14,7 @@ def near_to_farfield(E_near: field_t, dx: float, dy: float, padded_size: List[int] = None - ) -> Dict[str]: + ) -> Dict[str, Any]: """ Compute the farfield, i.e. the distribution of the fields after propagation through several wavelengths of uniform medium. @@ -122,7 +122,7 @@ def far_to_nearfield(E_far: field_t, dkx: float, dky: float, padded_size: List[int] = None - ) -> Dict[str]: + ) -> Dict[str, Any]: """ Compute the farfield, i.e. the distribution of the fields after propagation through several wavelengths of uniform medium. diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py index bd31192..74b66da 100644 --- a/meanas/fdfd/functional.py +++ b/meanas/fdfd/functional.py @@ -2,88 +2,41 @@ Functional versions of many FDFD operators. These can be useful for performing FDFD calculations without needing to construct large matrices in memory. -The functions generated here expect field inputs with shape (3, X, Y, Z), +The functions generated here expect `field_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 +from typing import List, Callable, Tuple import numpy from .. import dx_lists_t, field_t +from ..fdmath.functional import curl_forward, curl_back __author__ = 'Jan Petykiewicz' -functional_matrix = Callable[[field_t], field_t] - - -def curl_h(dxes: dx_lists_t) -> functional_matrix: - """ - Curl operator for use with the H field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Function for taking the discretized curl of the H-field, F(H) -> curlH - """ - dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij') - - def dh(f, ax): - return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax] - - def ch_fun(h: field_t) -> field_t: - e = numpy.empty_like(h) - e[0] = dh(h[2], 1) - e[0] -= dh(h[1], 2) - e[1] = dh(h[0], 2) - e[1] -= dh(h[2], 0) - e[2] = dh(h[1], 0) - e[2] -= dh(h[0], 1) - return e - - return ch_fun - - -def curl_e(dxes: dx_lists_t) -> functional_matrix: - """ - Curl operator for use with the E field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Function for taking the discretized curl of the E-field, F(E) -> curlE - """ - dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij') - - def de(f, ax): - return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax] - - def ce_fun(e: field_t) -> field_t: - h = numpy.empty_like(e) - h[0] = de(e[2], 1) - h[0] -= de(e[1], 2) - h[1] = de(e[0], 2) - h[1] -= de(e[2], 0) - h[2] = de(e[1], 0) - h[2] -= de(e[0], 1) - return h - - return ce_fun +field_transform_t = Callable[[field_t], field_t] def e_full(omega: complex, dxes: dx_lists_t, epsilon: field_t, mu: field_t = None - ) -> functional_matrix: + ) -> field_transform_t: """ - Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field, - with wave equation - (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J + Wave operator for use with E-field. See `operators.e_full` for details. - :param omega: Angular frequency of the simulation - :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 + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters [dx_e, dx_h] as described in meanas.types + epsilon: Dielectric constant + mu: Magnetic permeability (default 1 everywhere) + + Return: + Function `f` implementing the wave operator + `f(E)` -> `-i * omega * J` """ - ch = curl_h(dxes) - ce = curl_e(dxes) + ch = curl_back(dxes[1]) + ce = curl_forward(dxes[0]) def op_1(e): curls = ch(ce(e)) @@ -103,18 +56,23 @@ def eh_full(omega: complex, dxes: dx_lists_t, epsilon: field_t, mu: field_t = None - ) -> functional_matrix: + ) -> Callable[[field_t, field_t], Tuple[field_t, field_t]]: """ Wave operator for full (both E and H) field representation. + See `operators.eh_full`. - :param omega: Angular frequency of the simulation - :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) + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters [dx_e, dx_h] as described in meanas.types + epsilon: Dielectric constant + mu: Magnetic permeability (default 1 everywhere) + + Returns: + Function `f` implementing the wave operator + `f(E, H)` -> `(J, -M)` """ - ch = curl_h(dxes) - ce = curl_e(dxes) + ch = curl_back(dxes[1]) + ce = curl_forward(dxes[0]) def op_1(e, h): return (ch(h) - 1j * omega * epsilon * e, @@ -133,23 +91,27 @@ def eh_full(omega: complex, def e2h(omega: complex, dxes: dx_lists_t, mu: field_t = None, - ) -> functional_matrix: + ) -> field_transform_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. + Utility operator for converting the `E` field into the `H` field. + For use with `e_full` -- assumes that there is no magnetic current `M`. - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param mu: Magnetic permeability (default 1 everywhere) - :return: Function for converting E to H - """ - A2 = curl_e(dxes) + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + mu: Magnetic permeability (default 1 everywhere) + + Return: + Function `f` for converting `E` to `H`, + `f(E)` -> `H` + """ + ce = curl_forward(dxes[0]) def e2h_1_1(e): - return A2(e) / (-1j * omega) + return ce(e) / (-1j * omega) def e2h_mu(e): - return A2(e) / (-1j * omega * mu) + return ce(e) / (-1j * omega * mu) if numpy.any(numpy.equal(mu, None)): return e2h_1_1 @@ -160,18 +122,22 @@ def e2h(omega: complex, def m2j(omega: complex, dxes: dx_lists_t, mu: field_t = None, - ) -> functional_matrix: + ) -> field_transform_t: """ - Utility operator for converting magnetic current (M) distribution - into equivalent electric current distribution (J). - For use with e.g. e_full(). + Utility operator for converting magnetic current `M` distribution + into equivalent electric current distribution `J`. + For use with e.g. `e_full`. - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param mu: Magnetic permeability (default 1 everywhere) - :return: Function for converting M to J - """ - ch = curl_h(dxes) + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + mu: Magnetic permeability (default 1 everywhere) + + Returns: + Function `f` for converting `M` to `J`, + `f(M)` -> `J` + """ + ch = curl_back(dxes[1]) def m2j_mu(m): J = ch(m / mu) / (-1j * omega) @@ -192,10 +158,23 @@ def e_tfsf_source(TF_region: field_t, dxes: dx_lists_t, epsilon: field_t, mu: field_t = None, - ) -> functional_matrix: + ) -> field_transform_t: """ - Operator that turuns an E-field distribution into a total-field/scattered-field + Operator that turns an E-field distribution into a total-field/scattered-field (TFSF) source. + + Args: + TF_region: mask which is set to 1 in the total-field region, and 0 elsewhere + (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` + epsilon: Dielectric constant distribution + mu: Magnetic permeability (default 1 everywhere) + + Returns: + Function `f` which takes an E field and returns a current distribution, + `f(E)` -> `J` """ # TODO documentation A = e_full(omega, dxes, epsilon, mu) @@ -205,7 +184,28 @@ def e_tfsf_source(TF_region: field_t, return neg_iwj / (-1j * omega) -def poynting_e_cross_h(dxes: dx_lists_t): +def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[field_t, field_t], field_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 + for the Poynting vector, \\( S = E \\times H \\) + + Note: + This function also shifts the input `E` field by one cell as required + for computing the Poynting cross product (see `meanas.fdfd` module docs). + + Note: + If `E` and `H` are peak amplitudes as assumed elsewhere in this code, + the time-average of the poynting vector is ` = Re(S)/2 = Re(E x H) / 2`. + The factor of `1/2` can be omitted if root-mean-square quantities are used + instead. + + Args: + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + + Returns: + Function `f` that returns E x H as required for the poynting vector. + """ def exh(e: field_t, h: field_t): s = numpy.empty_like(e) ex = e[0] * dxes[0][0][:, None, None] diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 905c23e..20cf96a 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -1,18 +1,19 @@ """ Sparse matrix operators for use with electromagnetic wave equations. -These functions return sparse-matrix (scipy.sparse.spmatrix) representations of +These functions return sparse-matrix (`scipy.sparse.spmatrix`) representations of a variety of operators, intended for use with E and H fields vectorized using the - meanas.vec() and .unvec() functions (column-major/Fortran ordering). + `meanas.vec()` and `meanas.unvec()` functions. -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). +E- and H-field values are defined on a Yee cell; `epsilon` values should be calculated for + cells centered at each E component (`mu` at each H component). -Many of these functions require a 'dxes' parameter, of type meanas.dx_lists_type; see -the meanas.types submodule for details. +Many of these functions require a `dxes` parameter, of type `dx_lists_t`; see +the `meanas.types` submodule for details. The following operators are included: + - E-only wave operator - H-only wave operator - EH wave operator @@ -20,8 +21,6 @@ The following operators are included: - E to H conversion - M to J conversion - Poynting cross products - -Also available: - Circular shifts - Discrete derivatives - Averaging operators @@ -33,6 +32,7 @@ import numpy import scipy.sparse as sparse from .. import vec, dx_lists_t, vfield_t +from ..fdmath.operators import shift_with_mirror, rotation, curl_forward, curl_back __author__ = 'Jan Petykiewicz' @@ -46,26 +46,35 @@ def e_full(omega: complex, pmc: vfield_t = None, ) -> sparse.spmatrix: """ - Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field, - with wave equation - (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J + Wave operator + $$ \\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\omega^2 \\epsilon $$ - To make this matrix symmetric, use the preconditions from e_full_preconditioners(). + del x (1/mu * del x) - omega**2 * epsilon - :param omega: Angular frequency of the simulation - :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 - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (ie, pec.size == epsilon.size) - :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (ie, pmc.size == epsilon.size) - :return: Sparse matrix containing the wave operator + for use with the E-field, with wave equation + $$ (\\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\omega^2 \\epsilon) E = -\\imath \\omega J $$ + + (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J + + To make this matrix symmetric, use the preconditioners from `e_full_preconditioners()`. + + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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 + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + + Returns: + Sparse matrix containing the wave operator. """ - ce = curl_e(dxes) - ch = curl_h(dxes) + ch = curl_back(dxes[1]) + ce = curl_forward(dxes[0]) if numpy.any(numpy.equal(pec, None)): pe = sparse.eye(epsilon.size) @@ -90,15 +99,18 @@ def e_full(omega: complex, def e_full_preconditioners(dxes: dx_lists_t ) -> Tuple[sparse.spmatrix, sparse.spmatrix]: """ - Left and right preconditioners (Pl, Pr) for symmetrizing the e_full wave operator. + Left and right preconditioners `(Pl, Pr)` for symmetrizing the `e_full` wave operator. - The preconditioned matrix A_symm = (Pl @ A @ Pr) is complex-symmetric + The preconditioned matrix `A_symm = (Pl @ A @ Pr)` is complex-symmetric (non-Hermitian unless there is no loss or PMLs). - The preconditioner matrices are diagonal and complex, with Pr = 1 / Pl + The preconditioner matrices are diagonal and complex, with `Pr = 1 / Pl` - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Preconditioner matrices (Pl, Pr) + Args: + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + + Returns: + Preconditioner matrices `(Pl, Pr)`. """ p_squared = [dxes[0][0][:, None, None] * dxes[1][1][None, :, None] * dxes[1][2][None, None, :], dxes[1][0][:, None, None] * dxes[0][1][None, :, None] * dxes[1][2][None, None, :], @@ -118,24 +130,33 @@ def h_full(omega: complex, pmc: vfield_t = None, ) -> sparse.spmatrix: """ - Wave operator del x (1/epsilon * del x) - omega**2 * mu, for use with H-field, - with wave equation - (del x (1/epsilon * del x) - omega**2 * mu) H = i * omega * M + Wave operator + $$ \\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu $$ - :param omega: Angular frequency of the simulation - :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 - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (ie, pec.size == epsilon.size) - :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (ie, pmc.size == epsilon.size) - :return: Sparse matrix containing the wave operator + del x (1/epsilon * del x) - omega**2 * mu + + for use with the H-field, with wave equation + $$ (\\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu) E = \\imath \\omega M $$ + + (del x (1/epsilon * del x) - omega**2 * mu) E = i * omega * M + + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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 + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + + Returns: + Sparse matrix containing the wave operator. """ - ec = curl_e(dxes) - hc = curl_h(dxes) + ch = curl_back(dxes[1]) + ce = curl_forward(dxes[0]) if numpy.any(numpy.equal(pec, None)): pe = sparse.eye(epsilon.size) @@ -153,7 +174,7 @@ def h_full(omega: complex, else: m = sparse.diags(mu) - A = pm @ (ec @ pe @ e_div @ hc - omega**2 * m) @ pm + A = pm @ (ce @ pe @ e_div @ ch - omega**2 * m) @ pm return A @@ -165,24 +186,42 @@ def eh_full(omega: complex, pmc: vfield_t = None ) -> sparse.spmatrix: """ - Wave operator for [E, H] field representation. This operator implements Maxwell's + Wave operator for `[E, H]` field representation. This operator implements Maxwell's equations without cancelling out either E or H. The operator is - [[-i * omega * epsilon, del x], - [del x, i * omega * mu]] + $$ \\begin{bmatrix} + -\\imath \\omega \\epsilon & \\nabla \\times \\\\ + \\nabla \\times & \\imath \\omega \\mu + \\end{bmatrix} $$ - for use with a field vector of the form hstack(vec(E), vec(H)). + [[-i * omega * epsilon, del x ], + [del x, i * omega * mu]] - :param omega: Angular frequency of the simulation - :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 - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (i.e., pec.size == epsilon.size) - :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (i.e., pmc.size == epsilon.size) - :return: Sparse matrix containing the wave operator + for use with a field vector of the form `cat(vec(E), vec(H))`: + $$ \\begin{bmatrix} + -\\imath \\omega \\epsilon & \\nabla \\times \\\\ + \\nabla \\times & \\imath \\omega \\mu + \\end{bmatrix} + \\begin{bmatrix} E \\\\ + H + \\end{bmatrix} + = \\begin{bmatrix} J \\\\ + -M + \\end{bmatrix} $$ + + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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 + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + + Returns: + Sparse matrix containing the wave operator. """ if numpy.any(numpy.equal(pec, None)): pe = sparse.eye(epsilon.size) @@ -200,34 +239,14 @@ def eh_full(omega: complex, iwm *= sparse.diags(mu) iwm = pm @ iwm @ pm - A1 = pe @ curl_h(dxes) @ pm - A2 = pm @ curl_e(dxes) @ pe + A1 = pe @ curl_back(dxes[1]) @ pm + A2 = pm @ curl_forward(dxes[0]) @ pe A = sparse.bmat([[-iwe, A1], [A2, iwm]]) return A -def curl_h(dxes: dx_lists_t) -> sparse.spmatrix: - """ - Curl operator for use with the H field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Sparse matrix for taking the discretized curl of the H-field - """ - return cross(deriv_back(dxes[1])) - - -def curl_e(dxes: dx_lists_t) -> sparse.spmatrix: - """ - Curl operator for use with the E field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Sparse matrix for taking the discretized curl of the E-field - """ - return cross(deriv_forward(dxes[0])) - - def e2h(omega: complex, dxes: dx_lists_t, mu: vfield_t = None, @@ -235,17 +254,20 @@ def e2h(omega: complex, ) -> sparse.spmatrix: """ Utility operator for converting the E field into the H field. - For use with e_full -- assumes that there is no magnetic current M. + For use with `e_full()` -- assumes that there is no magnetic current M. - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param 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). - The PMC is applied per-field-component (ie, pmc.size == epsilon.size) - :return: Sparse matrix for converting E to H + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + + Returns: + Sparse matrix for converting E to H. """ - op = curl_e(dxes) / (-1j * omega) + op = curl_forward(dxes[0]) / (-1j * omega) if not numpy.any(numpy.equal(mu, None)): op = sparse.diags(1 / mu) @ op @@ -261,16 +283,18 @@ def m2j(omega: complex, mu: vfield_t = None ) -> sparse.spmatrix: """ - Utility operator for converting M field into J. - Converts a magnetic current M into an electric current J. - For use with eg. e_full. + Operator for converting a magnetic current M into an electric current J. + For use with eg. `e_full()`. - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param mu: Vectorized magnetic permeability (default 1 everywhere) - :return: Sparse matrix for converting E to H + Args: + omega: Angular frequency of the simulation + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + mu: Vectorized magnetic permeability (default 1 everywhere) + + Returns: + Sparse matrix for converting M to J. """ - op = curl_h(dxes) / (1j * omega) + op = curl_back(dxes[1]) / (1j * omega) if not numpy.any(numpy.equal(mu, None)): op = op @ sparse.diags(1 / mu) @@ -278,178 +302,17 @@ def m2j(omega: complex, return op -def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: - """ - Utility operator for performing a circular shift along a specified axis by a - specified number of elements. - - :param axis: Axis to shift along. x=0, y=1, z=2 - :param shape: Shape of the grid being shifted - :param shift_distance: Number of cells to shift by. May be negative. Default 1. - :return: Sparse matrix for performing the circular shift - """ - if len(shape) not in (2, 3): - raise Exception('Invalid shape: {}'.format(shape)) - if axis not in range(len(shape)): - raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape)) - - shifts = [abs(shift_distance) if a == axis else 0 for a in range(3)] - shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)] - ijk = numpy.meshgrid(*shifted_diags, indexing='ij') - - n = numpy.prod(shape) - i_ind = numpy.arange(n) - j_ind = numpy.ravel_multi_index(ijk, shape, order='C') - - vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) - - d = sparse.csr_matrix(vij, shape=(n, n)) - - if shift_distance < 0: - d = d.T - - return d - - -def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: - """ - Utility operator for performing an n-element shift along a specified axis, with mirror - boundary conditions applied to the cells beyond the receding edge. - - :param axis: Axis to shift along. x=0, y=1, z=2 - :param shape: Shape of the grid being shifted - :param shift_distance: Number of cells to shift by. May be negative. Default 1. - :return: Sparse matrix for performing the circular shift - """ - if len(shape) not in (2, 3): - raise Exception('Invalid shape: {}'.format(shape)) - if axis not in range(len(shape)): - raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape)) - if shift_distance >= shape[axis]: - raise Exception('Shift ({}) is too large for axis {} of size {}'.format( - shift_distance, axis, shape[axis])) - - def mirrored_range(n, s): - v = numpy.arange(n) + s - v = numpy.where(v >= n, 2 * n - v - 1, v) - v = numpy.where(v < 0, - 1 - v, v) - return v - - shifts = [shift_distance if a == axis else 0 for a in range(3)] - shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts)] - ijk = numpy.meshgrid(*shifted_diags, indexing='ij') - - n = numpy.prod(shape) - i_ind = numpy.arange(n) - j_ind = numpy.ravel_multi_index(ijk, shape, order='C') - - vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) - - d = sparse.csr_matrix(vij, shape=(n, n)) - return d - - -def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]: - """ - Utility operators for taking discretized derivatives (forward variant). - - :param dx_e: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...]. - :return: List of operators for taking forward derivatives along each axis. - """ - shape = [s.size for s in dx_e] - n = numpy.prod(shape) - - dx_e_expanded = numpy.meshgrid(*dx_e, indexing='ij') - - def deriv(axis): - return rotation(axis, shape, 1) - sparse.eye(n) - - Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a) - for a, dx in enumerate(dx_e_expanded)] - - return Ds - - -def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]: - """ - Utility operators for taking discretized derivatives (backward variant). - - :param dx_h: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...]. - :return: List of operators for taking forward derivatives along each axis. - """ - shape = [s.size for s in dx_h] - n = numpy.prod(shape) - - dx_h_expanded = numpy.meshgrid(*dx_h, indexing='ij') - - def deriv(axis): - return rotation(axis, shape, -1) - sparse.eye(n) - - Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a) - for a, dx in enumerate(dx_h_expanded)] - - return Ds - - -def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix: - """ - Cross product operator - - :param B: List [Bx, By, Bz] of sparse matrices corresponding to the x, y, z - portions of the operator on the left side of the cross product. - :return: Sparse matrix corresponding to (B x), where x is the cross product - """ - n = B[0].shape[0] - zero = sparse.csr_matrix((n, n)) - return sparse.bmat([[zero, -B[2], B[1]], - [B[2], zero, -B[0]], - [-B[1], B[0], zero]]) - - -def vec_cross(b: vfield_t) -> sparse.spmatrix: - """ - Vector cross product operator - - :param b: Vector on the left side of the cross product - :return: Sparse matrix corresponding to (b x), where x is the cross product - """ - B = [sparse.diags(c) for c in numpy.split(b, 3)] - return cross(B) - - -def avgf(axis: int, shape: List[int]) -> sparse.spmatrix: - """ - Forward average operator (x4 = (x4 + x5) / 2) - - :param axis: Axis to average along (x=0, y=1, z=2) - :param shape: Shape of the grid to average - :return: Sparse matrix for forward average operation - """ - if len(shape) not in (2, 3): - raise Exception('Invalid shape: {}'.format(shape)) - - n = numpy.prod(shape) - return 0.5 * (sparse.eye(n) + rotation(axis, shape)) - - -def avgb(axis: int, shape: List[int]) -> sparse.spmatrix: - """ - Backward average operator (x4 = (x4 + x3) / 2) - - :param axis: Axis to average along (x=0, y=1, z=2) - :param shape: Shape of the grid to average - :return: Sparse matrix for backward average operation - """ - return avgf(axis, shape).T - - def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: """ - Operator for computing the Poynting vector, containing the (E x) portion of the Poynting vector. + 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 meanas.types - :return: Sparse matrix containing (E x) portion of Poynting cross product + Args: + e: Vectorized E-field for the ExH cross product + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + + Returns: + Sparse matrix containing (E x) portion of Poynting cross product. """ shape = [len(dx) for dx in dxes[0]] @@ -472,9 +335,12 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: """ Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector. - :param h: Vectorized H-field for the HxE cross product - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Sparse matrix containing (H x) portion of Poynting cross product + Args: + h: Vectorized H-field for the HxE cross product + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` + + Returns: + Sparse matrix containing (H x) portion of Poynting cross product. """ shape = [len(dx) for dx in dxes[0]] @@ -499,8 +365,22 @@ def e_tfsf_source(TF_region: vfield_t, mu: vfield_t = None, ) -> sparse.spmatrix: """ - Operator that turns an E-field distribution into a total-field/scattered-field - (TFSF) source. + Operator that turns a desired E-field distribution into a + total-field/scattered-field (TFSF) source. + + TODO: Reference Rumpf paper + + Args: + 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` + epsilon: Vectorized dielectric constant + mu: Vectorized magnetic permeability (default 1 everywhere). + + Returns: + Sparse matrix that turns an E-field into a current (J) distribution. + """ # TODO documentation A = e_full(omega, dxes, epsilon, mu) @@ -518,7 +398,19 @@ def e_boundary_source(mask: vfield_t, """ Operator that turns an E-field distrubtion into a current (J) distribution along the edges (external and internal) of the provided mask. This is just an - e_tfsf_source with an additional masking step. + `e_tfsf_source()` with an additional masking step. + + Args: + mask: The current distribution is generated at the edges of the mask, + 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` + epsilon: Vectorized dielectric constant + mu: Vectorized magnetic permeability (default 1 everywhere). + + Returns: + Sparse matrix that turns an E-field into a current (J) distribution. """ full = e_tfsf_source(TF_region=mask, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) diff --git a/meanas/fdfd/scpml.py b/meanas/fdfd/scpml.py index 897d43a..c9a93e6 100644 --- a/meanas/fdfd/scpml.py +++ b/meanas/fdfd/scpml.py @@ -1,5 +1,5 @@ """ -Functions for creating stretched coordinate PMLs. +Functions for creating stretched coordinate perfectly matched layer (PML) absorbers. """ from typing import List, Callable @@ -7,24 +7,29 @@ import numpy from .. import dx_lists_t + __author__ = 'Jan Petykiewicz' -s_function_type = Callable[[float], float] +s_function_t = Callable[[float], float] +"""Typedef for s-functions""" def prepare_s_function(ln_R: float = -16, m: float = 4 - ) -> s_function_type: + ) -> s_function_t: """ Create an s_function to pass to the SCPML functions. This is used when you would like to customize the PML parameters. - :param ln_R: Natural logarithm of the desired reflectance - :param m: Polynomial order for the PML (imaginary part increases as distance ** m) - :return: An s_function, which takes an ndarray (distances) and returns an ndarray (complex part - of the cell width; needs to be divided by sqrt(epilon_effective) * real(omega)) - before use. + Args: + ln_R: Natural logarithm of the desired reflectance + m: Polynomial order for the PML (imaginary part increases as distance ** m) + + Returns: + An s_function, which takes an ndarray (distances) and returns an ndarray (complex part + of the cell width; needs to be divided by `sqrt(epilon_effective) * real(omega))` + before use. """ def s_factor(distance: numpy.ndarray) -> numpy.ndarray: s_max = (m + 1) * ln_R / 2 # / 2 because we assume periodic boundaries @@ -36,26 +41,29 @@ def uniform_grid_scpml(shape: numpy.ndarray or List[int], thicknesses: numpy.ndarray or List[int], omega: float, epsilon_effective: float = 1.0, - s_function: s_function_type = None, + s_function: s_function_t = None, ) -> dx_lists_t: """ Create dx arrays for a uniform grid with a cell width of 1 and a pml. - If you want something more fine-grained, check out stretch_with_scpml(...). + If you want something more fine-grained, check out `stretch_with_scpml(...)`. - :param shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick) - :param thicknesses: [th_x, th_y, th_z] Thickness of the PML in each direction. - Both polarities are added. - Each th_ of pml is applied twice, once on each edge of the grid along the given axis. - th_* may be zero, in which case no pml is added. - :param omega: Angular frequency for the simulation - :param epsilon_effective: Effective epsilon of the PML. Match this to the material - at the edge of your grid. - Default 1. - :param s_function: s_function created by prepare_s_function(...), allowing - customization of pml parameters. - Default uses prepare_s_function() with no parameters. - :return: Complex cell widths (dx_lists) + Args: + shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick) + thicknesses: `[th_x, th_y, th_z]` + Thickness of the PML in each direction. + Both polarities are added. + Each th_ of pml is applied twice, once on each edge of the grid along the given axis. + `th_*` may be zero, in which case no pml is added. + omega: Angular frequency for the simulation + epsilon_effective: Effective epsilon of the PML. Match this to the material + at the edge of your grid. + Default 1. + s_function: created by `prepare_s_function(...)`, allowing customization of pml parameters. + Default uses `prepare_s_function()` with no parameters. + + Returns: + Complex cell widths (dx_lists_t) as discussed in `meanas.types`. """ if s_function is None: s_function = prepare_s_function() @@ -88,21 +96,25 @@ def stretch_with_scpml(dxes: dx_lists_t, omega: float, epsilon_effective: float = 1.0, thickness: int = 10, - s_function: s_function_type = None, + s_function: s_function_t = None, ) -> dx_lists_t: """ Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis. - :param dxes: dx_tuple with coordinates to stretch - :param axis: axis to stretch (0=x, 1=y, 2=z) - :param polarity: direction to stretch (-1 for -ve, +1 for +ve) - :param omega: Angular frequency for the simulation - :param epsilon_effective: Effective epsilon of the PML. Match this to the material at the - edge of your grid. Default 1. - :param thickness: number of cells to use for pml (default 10) - :param s_function: s_function created by prepare_s_function(...), allowing customization - of pml parameters. Default uses prepare_s_function() with no parameters. - :return: Complex cell widths + Args: + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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 + epsilon_effective: Effective epsilon of the PML. Match this to the material at the + edge of your grid. Default 1. + thickness: number of cells to use for pml (default 10) + s_function: Created by `prepare_s_function(...)`, allowing customization + of pml parameters. Default uses `prepare_s_function()` with no parameters. + + Returns: + Complex cell widths (dx_lists_t) as discussed in `meanas.types`. + Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed. """ if s_function is None: s_function = prepare_s_function() @@ -147,25 +159,3 @@ def stretch_with_scpml(dxes: dx_lists_t, dxes[1][axis] = dx_bi return dxes - - -def generate_periodic_dx(pos: List[numpy.ndarray]) -> dx_lists_t: - """ - Given a list of 3 ndarrays cell centers, creates the cell width parameters for a periodic grid. - - :param pos: List of 3 ndarrays of cell centers - :return: (dx_a, dx_b) cell widths (no pml) - """ - if len(pos) != 3: - raise Exception('Must have len(pos) == 3') - - dx_a = [numpy.array(numpy.inf)] * 3 - dx_b = [numpy.array(numpy.inf)] * 3 - - for i, p_orig in enumerate(pos): - p = numpy.array(p_orig, dtype=float) - if p.size != 1: - p_shifted = numpy.hstack((p[1:], p[-1] + (p[1] - p[0]))) - dx_a[i] = numpy.diff(p) - dx_b[i] = numpy.diff((p + p_shifted) / 2) - return dx_a, dx_b diff --git a/meanas/fdfd/solvers.py b/meanas/fdfd/solvers.py index a0ce403..aa9633a 100644 --- a/meanas/fdfd/solvers.py +++ b/meanas/fdfd/solvers.py @@ -1,5 +1,5 @@ """ -Solvers for FDFD problems. +Solvers and solver interface for FDFD problems. """ from typing import List, Callable, Dict, Any @@ -22,10 +22,13 @@ def _scipy_qmr(A: scipy.sparse.csr_matrix, """ Wrapper for scipy.sparse.linalg.qmr - :param A: Sparse matrix - :param b: Right-hand-side vector - :param kwargs: Passed as **kwargs to the wrapped function - :return: Guess for solution (returned even if didn't converge) + Args: + A: Sparse matrix + b: Right-hand-side vector + kwargs: Passed as **kwargs to the wrapped function + + Returns: + Guess for solution (returned even if didn't converge) """ ''' @@ -70,27 +73,31 @@ def generic(omega: complex, """ Conjugate gradient FDFD solver using CSR sparse matrices. - All ndarray arguments should be 1D array, as returned by meanas.vec(). + All ndarray arguments should be 1D arrays, 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) - :param J: Electric current distribution (at E-field locations) - :param epsilon: Dielectric constant distribution (at E-field locations) - :param mu: Magnetic permeability distribution (at H-field locations) - :param pec: Perfect electric conductor distribution - (at E-field locations; non-zero value indicates PEC is present) - :param pmc: Perfect magnetic conductor distribution - (at H-field locations; non-zero value indicates PMC is present) - :param adjoint: If true, solves the adjoint problem. - :param matrix_solver: Called as matrix_solver(A, b, **matrix_solver_opts) -> x - Where A: scipy.sparse.csr_matrix - b: numpy.ndarray - x: numpy.ndarray - Default is a wrapped version of scipy.sparse.linalg.qmr() - which doesn't return convergence info and logs the residual - every 100 iterations. - :param matrix_solver_opts: Passed as kwargs to matrix_solver(...) - :return: E-field which solves the system. + 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` + J: Electric current distribution (at E-field locations) + epsilon: Dielectric constant distribution (at E-field locations) + mu: Magnetic permeability distribution (at H-field locations) + pec: Perfect electric conductor distribution + (at E-field locations; non-zero value indicates PEC is present) + pmc: Perfect magnetic conductor distribution + (at H-field locations; non-zero value indicates PMC is present) + adjoint: If true, solves the adjoint problem. + matrix_solver: Called as `matrix_solver(A, b, **matrix_solver_opts) -> x`, + where `A`: `scipy.sparse.csr_matrix`; + `b`: `numpy.ndarray`; + `x`: `numpy.ndarray`; + Default is a wrapped version of `scipy.sparse.linalg.qmr()` + which doesn't return convergence info and logs the residual + every 100 iterations. + matrix_solver_opts: Passed as kwargs to `matrix_solver(...)` + + Returns: + E-field which solves the system. """ if matrix_solver_opts is None: diff --git a/meanas/fdfd/waveguide.py b/meanas/fdfd/waveguide.py deleted file mode 100644 index c24a471..0000000 --- a/meanas/fdfd/waveguide.py +++ /dev/null @@ -1,492 +0,0 @@ -""" -Various operators and helper functions for solving for waveguide modes. - -Assuming a z-dependence of the from exp(-i * wavenumber * z), we can simplify Maxwell's - equations in the absence of sources to the form - -A @ [H_x, H_y] = wavenumber**2 * [H_x, H_y] - -with A = -omega**2 * epsilon * mu + -epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] + -[[Dx], [Dy]] / mu * [Dx, Dy] * mu - -which is the form used in this file. - -As the z-dependence is known, all the functions in this file assume a 2D grid - (ie. dxes = [[[dx_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dx_h_0, ...], [dy_h_0, ...]]]) - with propagation along the z axis. -""" -# TODO update module docs - -from typing import List, Tuple -import numpy -from numpy.linalg import norm -import scipy.sparse as sparse - -from .. import vec, unvec, dx_lists_t, field_t, vfield_t -from . import operators - - -__author__ = 'Jan Petykiewicz' - - -def operator_e(omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - ) -> sparse.spmatrix: - if numpy.any(numpy.equal(mu, None)): - mu = numpy.ones_like(epsilon) - - Dfx, Dfy = operators.deriv_forward(dxes[0]) - Dbx, Dby = operators.deriv_back(dxes[1]) - - eps_parts = numpy.split(epsilon, 3) - eps_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1]))) - eps_z_inv = sparse.diags(1 / eps_parts[2]) - - mu_parts = numpy.split(mu, 3) - mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0]))) - mu_z_inv = sparse.diags(1 / mu_parts[2]) - - op = omega * omega * mu_yx @ eps_xy + \ - mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \ - sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy - return op - - -def operator_h(omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - ) -> sparse.spmatrix: - """ - Waveguide operator of the form - - omega**2 * epsilon * mu + - epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] + - [[Dx], [Dy]] / mu * [Dx, Dy] * mu - - for use with a field vector of the form [H_x, H_y]. - - This operator can be used to form an eigenvalue problem of the form - A @ [H_x, H_y] = wavenumber**2 * [H_x, H_y] - - which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * z) - z-dependence is assumed for the fields). - - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Sparse matrix representation of the operator - """ - if numpy.any(numpy.equal(mu, None)): - mu = numpy.ones_like(epsilon) - - Dfx, Dfy = operators.deriv_forward(dxes[0]) - Dbx, Dby = operators.deriv_back(dxes[1]) - - eps_parts = numpy.split(epsilon, 3) - eps_yx = sparse.diags(numpy.hstack((eps_parts[1], eps_parts[0]))) - eps_z_inv = sparse.diags(1 / eps_parts[2]) - - mu_parts = numpy.split(mu, 3) - mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1]))) - mu_z_inv = sparse.diags(1 / mu_parts[2]) - - op = omega * omega * eps_yx @ mu_xy + \ - eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \ - sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy - - return op - - -def normalized_fields_e(e_xy: numpy.ndarray, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: - """ - Given a vector e_xy containing the vectorized E_x and E_y fields, - returns normalized, vectorized E and H fields for the system. - - :param e_xy: Vector containing E_x and E_y fields - :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy` - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :param prop_phase: Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction. - Default 0 (continuous propagation direction, i.e. dz->0). - :return: Normalized, vectorized (e, h) containing all vector components. - """ - e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy - h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy - e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, - mu=mu, prop_phase=prop_phase) - return e_norm, h_norm - - -def normalized_fields_h(h_xy: numpy.ndarray, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: - """ - Given a vector e_xy containing the vectorized E_x and E_y fields, - returns normalized, vectorized E and H fields for the system. - - :param e_xy: Vector containing E_x and E_y fields - :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy` - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :param dxes_prop: Grid cell width in the propagation direction. Default 0 (continuous). - :return: Normalized, vectorized (e, h) containing all vector components. - """ - e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy - h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy - e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, - mu=mu, prop_phase=prop_phase) - return e_norm, h_norm - - -def _normalized_fields(e: numpy.ndarray, - h: numpy.ndarray, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - prop_phase: float = 0, - ) -> Tuple[vfield_t, vfield_t]: - # TODO documentation - shape = [s.size for s in dxes[0]] - dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)] - - E = unvec(e, shape) - H = unvec(h, shape) - - # Find time-averaged Sz and normalize to it - # H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting - phase = numpy.exp(-1j * -prop_phase / 2) - Sz_a = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0] - Sz_b = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1] - Sz_tavg = numpy.real(Sz_a.sum() - Sz_b.sum()) * 0.5 # 0.5 since E, H are assumed to be peak (not RMS) amplitudes - assert Sz_tavg > 0, 'Found a mode propagating in the wrong direction! Sz_tavg={}'.format(Sz_tavg) - - energy = epsilon * e.conj() * e - - norm_amplitude = 1 / numpy.sqrt(Sz_tavg) - norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric - - # Try to break symmetry to assign a consistent sign [experimental TODO] - E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) - sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) - - norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle) - - e *= norm_factor - h *= norm_factor - - return e, h - - -def exy2h(wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None - ) -> sparse.spmatrix: - """ - Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields, - into a vectorized H containing all three H components - - :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy` - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Sparse matrix representing the operator - """ - e2hop = e2h(wavenumber=wavenumber, omega=omega, dxes=dxes, mu=mu) - return e2hop @ exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) - - -def hxy2e(wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None - ) -> sparse.spmatrix: - """ - Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields, - into a vectorized E containing all three E components - - :param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy` - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Sparse matrix representing the operator - """ - h2eop = h2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon) - return h2eop @ hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) - - -def hxy2h(wavenumber: complex, - dxes: dx_lists_t, - mu: vfield_t = None - ) -> sparse.spmatrix: - """ - Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields, - into a vectorized H containing all three H components - - :param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy` - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Sparse matrix representing the operator - """ - Dfx, Dfy = operators.deriv_forward(dxes[0]) - hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber) - - if not numpy.any(numpy.equal(mu, None)): - mu_parts = numpy.split(mu, 3) - mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1]))) - mu_z_inv = sparse.diags(1 / mu_parts[2]) - - hxy2hz = mu_z_inv @ hxy2hz @ mu_xy - - n_pts = dxes[1][0].size * dxes[1][1].size - op = sparse.vstack((sparse.eye(2 * n_pts), - hxy2hz)) - return op - - -def exy2e(wavenumber: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - ) -> sparse.spmatrix: - """ - Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields, - into a vectorized E containing all three E components - - :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy` - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :return: Sparse matrix representing the operator - """ - Dbx, Dby = operators.deriv_back(dxes[1]) - exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber) - - if not numpy.any(numpy.equal(epsilon, None)): - epsilon_parts = numpy.split(epsilon, 3) - epsilon_xy = sparse.diags(numpy.hstack((epsilon_parts[0], epsilon_parts[1]))) - epsilon_z_inv = sparse.diags(1 / epsilon_parts[2]) - - exy2ez = epsilon_z_inv @ exy2ez @ epsilon_xy - - n_pts = dxes[0][0].size * dxes[0][1].size - op = sparse.vstack((sparse.eye(2 * n_pts), - exy2ez)) - return op - - -def e2h(wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - mu: vfield_t = None - ) -> sparse.spmatrix: - """ - Returns an operator which, when applied to a vectorized E eigenfield, produces - the vectorized H eigenfield. - - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Sparse matrix representation of the operator - """ - op = curl_e(wavenumber, dxes) / (-1j * omega) - if not numpy.any(numpy.equal(mu, None)): - op = sparse.diags(1 / mu) @ op - return op - - -def h2e(wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t - ) -> sparse.spmatrix: - """ - Returns an operator which, when applied to a vectorized H eigenfield, produces - the vectorized E eigenfield. - - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :return: Sparse matrix representation of the operator - """ - op = sparse.diags(1 / (1j * omega * epsilon)) @ curl_h(wavenumber, dxes) - return op - - -def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: - """ - Discretized curl operator for use with the waveguide E field. - - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :return: Sparse matrix representation of the operator - """ - n = 1 - for d in dxes[0]: - n *= len(d) - - Bz = -1j * wavenumber * sparse.eye(n) - Dfx, Dfy = operators.deriv_forward(dxes[0]) - return operators.cross([Dfx, Dfy, Bz]) - - -def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: - """ - Discretized curl operator for use with the waveguide H field. - - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :return: Sparse matrix representation of the operator - """ - n = 1 - for d in dxes[1]: - n *= len(d) - - Bz = -1j * wavenumber * sparse.eye(n) - Dbx, Dby = operators.deriv_back(dxes[1]) - return operators.cross([Dbx, Dby, Bz]) - - -def h_err(h: vfield_t, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None - ) -> float: - """ - Calculates the relative error in the H field - - :param h: Vectorized H field - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Relative error norm(OP @ h) / norm(h) - """ - ce = curl_e(wavenumber, dxes) - ch = curl_h(wavenumber, dxes) - - eps_inv = sparse.diags(1 / epsilon) - - if numpy.any(numpy.equal(mu, None)): - op = ce @ eps_inv @ ch @ h - omega ** 2 * h - else: - op = ce @ eps_inv @ ch @ h - omega ** 2 * (mu * h) - - return norm(op) / norm(h) - - -def e_err(e: vfield_t, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None - ) -> float: - """ - Calculates the relative error in the E field - - :param e: Vectorized E field - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param epsilon: Vectorized dielectric constant grid - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Relative error norm(OP @ e) / norm(e) - """ - ce = curl_e(wavenumber, dxes) - ch = curl_h(wavenumber, dxes) - - if numpy.any(numpy.equal(mu, None)): - op = ch @ ce @ e - omega ** 2 * (epsilon * e) - else: - mu_inv = sparse.diags(1 / mu) - op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e) - - return norm(op) / norm(e) - - -def cylindrical_operator(omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - r0: float, - ) -> sparse.spmatrix: - """ - Cylindrical coordinate waveguide operator of the form - - TODO - - for use with a field vector of the form [E_r, E_y]. - - This operator can be used to form an eigenvalue problem of the form - A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y] - - which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * theta) - theta-dependence is assumed for the fields). - - :param omega: The angular frequency of the system - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D) - :param 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. - :return: Sparse matrix representation of the operator - """ - - Dfx, Dfy = operators.deriv_forward(dxes[0]) - Dbx, Dby = operators.deriv_back(dxes[1]) - - rx = r0 + numpy.cumsum(dxes[0][0]) - ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0]) - tx = rx/r0 - ty = ry/r0 - - Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1))) - Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1))) - - eps_parts = numpy.split(epsilon, 3) - eps_x = sparse.diags(eps_parts[0]) - eps_y = sparse.diags(eps_parts[1]) - eps_z_inv = sparse.diags(1 / eps_parts[2]) - - pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby)) - pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx)) - a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy - a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx - b0 = Dbx @ Ty @ Dfy - b1 = Dby @ Ty @ Dfx - - diag = sparse.block_diag - op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \ - - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1)) - - return op - diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py new file mode 100644 index 0000000..6a7f24f --- /dev/null +++ b/meanas/fdfd/waveguide_2d.py @@ -0,0 +1,607 @@ +""" +Operators and helper functions for waveguides with unchanging cross-section. + +The propagation direction is chosen to be along the z axis, and all fields +are given an implicit z-dependence of the form `exp(-1 * wavenumber * z)`. + +As the z-dependence is known, all the functions in this file assume a 2D grid + (i.e. `dxes = [[[dx_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dx_h_0, ...], [dy_h_0, ...]]]`). +""" +# TODO update module docs + +from typing import List, Tuple +import numpy +from numpy.linalg import norm +import scipy.sparse as sparse + +from .. import vec, unvec, dx_lists_t, field_t, vfield_t +from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration +from . import operators + + +__author__ = 'Jan Petykiewicz' + + +def operator_e(omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + ) -> sparse.spmatrix: + """ + Waveguide operator of the form + + omega**2 * mu * epsilon + + mu * [[-Dy], [Dx]] / mu * [-Dy, Dx] + + [[Dx], [Dy]] / epsilon * [Dx, Dy] * epsilon + + for use with a field vector of the form `cat([E_x, E_y])`. + + More precisely, the operator is + $$ \\omega^2 \\mu_{yx} \\epsilon_{xy} + + \\mu_{yx} \\begin{bmatrix} -D_{by} \\\\ + D_{bx} \\end{bmatrix} \\mu_z^{-1} + \\begin{bmatrix} -D_{fy} & D_{fx} \\end{bmatrix} + + \\begin{bmatrix} D_{fx} \\\\ + D_{fy} \\end{bmatrix} \\epsilon_z^{-1} + \\begin{bmatrix} D_{bx} & D_{by} \\end{bmatrix} \\epsilon_{xy} $$ + + where + \\( \\epsilon_{xy} = \\begin{bmatrix} + \\epsilon_x & 0 \\\\ + 0 & \\epsilon_y + \\end{bmatrix} \\), + \\( \\mu_{yx} = \\begin{bmatrix} + \\mu_y & 0 \\\\ + 0 & \\mu_x + \\end{bmatrix} \\), + \\( D_{fx} \\) and \\( D_{bx} \\) are the forward and backward derivatives along x, + and each \\( \\epsilon_x, \\mu_y, \\) etc. is a diagonal matrix representing + + + This operator can be used to form an eigenvalue problem of the form + `operator_e(...) @ [E_x, E_y] = wavenumber**2 * [E_x, E_y]` + + which can then be solved for the eigenmodes of the system (an `exp(-i * wavenumber * z)` + z-dependence is assumed for the fields). + + Args: + omega: The angular frequency of the system. + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representation of the operator. + """ + if numpy.any(numpy.equal(mu, None)): + mu = numpy.ones_like(epsilon) + + Dfx, Dfy = operators.deriv_forward(dxes[0]) + Dbx, Dby = operators.deriv_back(dxes[1]) + + eps_parts = numpy.split(epsilon, 3) + eps_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1]))) + eps_z_inv = sparse.diags(1 / eps_parts[2]) + + mu_parts = numpy.split(mu, 3) + mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0]))) + mu_z_inv = sparse.diags(1 / mu_parts[2]) + + op = omega * omega * mu_yx @ eps_xy + \ + mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \ + sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy + return op + + +def operator_h(omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + ) -> sparse.spmatrix: + """ + Waveguide operator of the form + + omega**2 * epsilon * mu + + epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] + + [[Dx], [Dy]] / mu * [Dx, Dy] * mu + + for use with a field vector of the form `cat([H_x, H_y])`. + + More precisely, the operator is + $$ \\omega^2 \\epsilon_{yx} \\mu_{xy} + + \\epsilon_{yx} \\begin{bmatrix} -D_{fy} \\\\ + D_{fx} \\end{bmatrix} \\epsilon_z^{-1} + \\begin{bmatrix} -D_{by} & D_{bx} \\end{bmatrix} + + \\begin{bmatrix} D_{bx} \\\\ + D_{by} \\end{bmatrix} \\mu_z^{-1} + \\begin{bmatrix} D_{fx} & D_{fy} \\end{bmatrix} \\mu_{xy} $$ + + where + \\( \\epsilon_{yx} = \\begin{bmatrix} + \\epsilon_y & 0 \\\\ + 0 & \\epsilon_x + \\end{bmatrix} \\), + \\( \\mu_{xy} = \\begin{bmatrix} + \\mu_x & 0 \\\\ + 0 & \\mu_y + \\end{bmatrix} \\), + \\( D_{fx} \\) and \\( D_{bx} \\) are the forward and backward derivatives along x, + and each \\( \\epsilon_x, \\mu_y, \\) etc. is a diagonal matrix. + + + This operator can be used to form an eigenvalue problem of the form + `operator_h(...) @ [H_x, H_y] = wavenumber**2 * [H_x, H_y]` + + which can then be solved for the eigenmodes of the system (an `exp(-i * wavenumber * z)` + z-dependence is assumed for the fields). + + Args: + omega: The angular frequency of the system. + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representation of the operator. + """ + if numpy.any(numpy.equal(mu, None)): + mu = numpy.ones_like(epsilon) + + Dfx, Dfy = operators.deriv_forward(dxes[0]) + Dbx, Dby = operators.deriv_back(dxes[1]) + + eps_parts = numpy.split(epsilon, 3) + eps_yx = sparse.diags(numpy.hstack((eps_parts[1], eps_parts[0]))) + eps_z_inv = sparse.diags(1 / eps_parts[2]) + + mu_parts = numpy.split(mu, 3) + mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1]))) + mu_z_inv = sparse.diags(1 / mu_parts[2]) + + op = omega * omega * eps_yx @ mu_xy + \ + eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \ + sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy + + return op + + +def normalized_fields_e(e_xy: numpy.ndarray, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + prop_phase: float = 0, + ) -> Tuple[vfield_t, vfield_t]: + """ + Given a vector `e_xy` containing the vectorized E_x and E_y fields, + returns normalized, vectorized E and H fields for the system. + + Args: + e_xy: Vector containing E_x and E_y fields + 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) + 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. + Default 0 (continuous propagation direction, i.e. dz->0). + + Returns: + `(e, h)`, where each field is vectorized, normalized, + and contains all three vector components. + """ + e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy + h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy + e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, + mu=mu, prop_phase=prop_phase) + return e_norm, h_norm + + +def normalized_fields_h(h_xy: numpy.ndarray, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + prop_phase: float = 0, + ) -> Tuple[vfield_t, vfield_t]: + """ + Given a vector `h_xy` containing the vectorized H_x and H_y fields, + returns normalized, vectorized E and H fields for the system. + + Args: + h_xy: Vector containing H_x and H_y fields + 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) + 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. + Default 0 (continuous propagation direction, i.e. dz->0). + + Returns: + `(e, h)`, where each field is vectorized, normalized, + and contains all three vector components. + """ + e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy + h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy + e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, + mu=mu, prop_phase=prop_phase) + return e_norm, h_norm + + +def _normalized_fields(e: numpy.ndarray, + h: numpy.ndarray, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + prop_phase: float = 0, + ) -> Tuple[vfield_t, vfield_t]: + # TODO documentation + shape = [s.size for s in dxes[0]] + dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)] + + E = unvec(e, shape) + H = unvec(h, shape) + + # Find time-averaged Sz and normalize to it + # H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting + phase = numpy.exp(-1j * -prop_phase / 2) + Sz_a = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0] + Sz_b = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1] + Sz_tavg = numpy.real(Sz_a.sum() - Sz_b.sum()) * 0.5 # 0.5 since E, H are assumed to be peak (not RMS) amplitudes + assert Sz_tavg > 0, 'Found a mode propagating in the wrong direction! Sz_tavg={}'.format(Sz_tavg) + + energy = epsilon * e.conj() * e + + norm_amplitude = 1 / numpy.sqrt(Sz_tavg) + norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric + + # Try to break symmetry to assign a consistent sign [experimental TODO] + E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) + sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) + + norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle) + + e *= norm_factor + h *= norm_factor + + return e, h + + +def exy2h(wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + into a vectorized H containing all three H components + + Args: + 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) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representing the operator. + """ + e2hop = e2h(wavenumber=wavenumber, omega=omega, dxes=dxes, mu=mu) + return e2hop @ exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) + + +def hxy2e(wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `h_xy` containing the vectorized H_x and H_y fields, + into a vectorized E containing all three E components + + Args: + 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) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representing the operator. + """ + h2eop = h2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon) + return h2eop @ hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) + + +def hxy2h(wavenumber: complex, + dxes: dx_lists_t, + mu: vfield_t = None + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `h_xy` containing the vectorized H_x and H_y fields, + into a vectorized H containing all three H components + + 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) + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representing the operator. + """ + Dfx, Dfy = operators.deriv_forward(dxes[0]) + hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber) + + if not numpy.any(numpy.equal(mu, None)): + mu_parts = numpy.split(mu, 3) + mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1]))) + mu_z_inv = sparse.diags(1 / mu_parts[2]) + + hxy2hz = mu_z_inv @ hxy2hz @ mu_xy + + n_pts = dxes[1][0].size * dxes[1][1].size + op = sparse.vstack((sparse.eye(2 * n_pts), + hxy2hz)) + return op + + +def exy2e(wavenumber: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + ) -> sparse.spmatrix: + """ + Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields, + into a vectorized E containing all three E components + + 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) + epsilon: Vectorized dielectric constant grid + + Returns: + Sparse matrix representing the operator. + """ + Dbx, Dby = operators.deriv_back(dxes[1]) + exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber) + + if not numpy.any(numpy.equal(epsilon, None)): + epsilon_parts = numpy.split(epsilon, 3) + epsilon_xy = sparse.diags(numpy.hstack((epsilon_parts[0], epsilon_parts[1]))) + epsilon_z_inv = sparse.diags(1 / epsilon_parts[2]) + + exy2ez = epsilon_z_inv @ exy2ez @ epsilon_xy + + n_pts = dxes[0][0].size * dxes[0][1].size + op = sparse.vstack((sparse.eye(2 * n_pts), + exy2ez)) + return op + + +def e2h(wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + mu: vfield_t = None + ) -> sparse.spmatrix: + """ + Returns an operator which, when applied to a vectorized E eigenfield, produces + the vectorized H eigenfield. + + 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) + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Sparse matrix representation of the operator. + """ + op = curl_e(wavenumber, dxes) / (-1j * omega) + if not numpy.any(numpy.equal(mu, None)): + op = sparse.diags(1 / mu) @ op + return op + + +def h2e(wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t + ) -> sparse.spmatrix: + """ + Returns an operator which, when applied to a vectorized H eigenfield, produces + the vectorized E eigenfield. + + 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) + epsilon: Vectorized dielectric constant grid + + Returns: + Sparse matrix representation of the operator. + """ + op = sparse.diags(1 / (1j * omega * epsilon)) @ curl_h(wavenumber, dxes) + return op + + +def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: + """ + Discretized curl operator for use with the waveguide E field. + + 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) + + Return: + Sparse matrix representation of the operator. + """ + n = 1 + for d in dxes[0]: + n *= len(d) + + Bz = -1j * wavenumber * sparse.eye(n) + Dfx, Dfy = operators.deriv_forward(dxes[0]) + return operators.cross([Dfx, Dfy, Bz]) + + +def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix: + """ + Discretized curl operator for use with the waveguide H field. + + 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) + + Return: + Sparse matrix representation of the operator. + """ + n = 1 + for d in dxes[1]: + n *= len(d) + + Bz = -1j * wavenumber * sparse.eye(n) + Dbx, Dby = operators.deriv_back(dxes[1]) + return operators.cross([Dbx, Dby, Bz]) + + +def h_err(h: vfield_t, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None + ) -> float: + """ + Calculates the relative error in the H field + + Args: + 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) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Relative error `norm(A_h @ h) / norm(h)`. + """ + ce = curl_e(wavenumber, dxes) + ch = curl_h(wavenumber, dxes) + + eps_inv = sparse.diags(1 / epsilon) + + if numpy.any(numpy.equal(mu, None)): + op = ce @ eps_inv @ ch @ h - omega ** 2 * h + else: + op = ce @ eps_inv @ ch @ h - omega ** 2 * (mu * h) + + return norm(op) / norm(h) + + +def e_err(e: vfield_t, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None + ) -> float: + """ + Calculates the relative error in the E field + + Args: + 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) + epsilon: Vectorized dielectric constant grid + mu: Vectorized magnetic permeability grid (default 1 everywhere) + + Returns: + Relative error `norm(A_e @ e) / norm(e)`. + """ + ce = curl_e(wavenumber, dxes) + ch = curl_h(wavenumber, dxes) + + if numpy.any(numpy.equal(mu, None)): + op = ch @ ce @ e - omega ** 2 * (epsilon * e) + else: + mu_inv = sparse.diags(1 / mu) + op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e) + + return norm(op) / norm(e) + + +def solve_modes(mode_numbers: List[int], + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + mode_margin: int = 2, + ) -> Tuple[List[vfield_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` + epsilon: Dielectric constant + mu: Magnetic permeability (default 1 everywhere) + mode_margin: The eigensolver will actually solve for `(max(mode_number) + mode_margin)` + modes, but only return the target mode. Increasing this value can improve the solver's + ability to find the correct mode. Default 2. + + Returns: + (e_xys, wavenumbers) + """ + + ''' + 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)) + + eigvals, eigvecs = signed_eigensolve(A_r, max(mode_number) + mode_margin) + e_xys = eigvecs[:, -(numpy.array(mode_number) + 1)] + + ''' + Now solve for the eigenvector of the full operator, using the real operator's + eigenvector as an initial guess for Rayleigh quotient iteration. + ''' + A = waveguide.operator_e(omega, dxes, epsilon, mu) + eigvals, e_xys = rayleigh_quotient_iteration(A, e_xys) + + # Calculate the wave-vector (force the real part to be positive) + wavenumbers = numpy.sqrt(eigvals) + wavenumbers *= numpy.sign(numpy.real(wavenumbers)) + + return e_xys, wavenumbers + + +def solve_mode(mode_number: int, + *args, + **kwargs + ) -> Tuple[vfield_t, complex]: + """ + Wrapper around `solve_modes()` that solves for a single mode. + + Args: + mode_number: 0-indexed mode number to solve for + *args: passed to `solve_modes()` + **kwargs: passed to `solve_modes()` + + Returns: + (e_xy, wavenumber) + """ + return solve_modes(mode_numbers=[mode_number], *args, **kwargs) diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py new file mode 100644 index 0000000..18727ce --- /dev/null +++ b/meanas/fdfd/waveguide_3d.py @@ -0,0 +1,236 @@ +""" +Tools for working with waveguide modes in 3D domains. + +This module relies heavily on `waveguide_2d` and mostly just transforms +its parameters into 2D equivalents and expands the results back into 3D. +""" +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 . import operators, waveguide_2d, functional + + +def solve_mode(mode_number: int, + omega: complex, + dxes: dx_lists_t, + axis: int, + polarity: int, + slices: List[slice], + epsilon: field_t, + mu: field_t = None, + ) -> Dict[str, complex or numpy.ndarray]: + """ + Given a 3D grid, selects a slice from the grid and attempts to + solve for an eigenmode propagating through that slice. + + 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` + 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 + as the waveguide cross-section. `slices[axis]` should select only one item. + epsilon: Dielectric constant + mu: Magnetic permeability (default 1 everywhere) + + Returns: + `{'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}` + """ + if mu is None: + mu = numpy.ones_like(epsilon) + + slices = tuple(slices) + + ''' + Solve the 2D problem in the specified plane + ''' + # Define rotation to set z as propagation direction + order = numpy.roll(range(3), 2 - axis) + reverse_order = numpy.roll(range(3), axis - 2) + + # Find dx in propagation direction + dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes]) + dx_prop = 0.5 * sum(dxab_forward)[0] + + # Reduce to 2D and solve the 2D problem + args_2d = { + 'omega': omega, + 'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes], + 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]), + 'mu': vec([mu[i][slices].transpose(order) for i in order]), + } + e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d) + + ''' + Apply corrections and expand to 3D + ''' + # Correct wavenumber to account for numerical dispersion. + wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2) + + shape = [d.size for d in args_2d['dxes'][0]] + ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber) + e = unvec(ve, shape) + h = unvec(vh, shape) + + # Adjust for propagation direction + h *= polarity + + # Apply phase shift to H-field + h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop) + e[2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop) + + # Expand E, H to full epsilon space we were given + E = numpy.zeros_like(epsilon, dtype=complex) + H = numpy.zeros_like(epsilon, dtype=complex) + for a, o in enumerate(reverse_order): + E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order) + H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order) + + results = { + 'wavenumber': wavenumber, + 'wavenumber_2d': wavenumber_2d, + 'H': H, + 'E': E, + } + return results + + +def compute_source(E: field_t, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + axis: int, + polarity: int, + slices: List[slice], + epsilon: field_t, + mu: field_t = None, + ) -> field_t: + """ + Given an eigenmode obtained by `solve_mode`, returns the current source distribution + necessary to position a unidirectional source at the slice location. + + Args: + 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` + 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 + as the waveguide cross-section. `slices[axis]` should select only one item. + mu: Magnetic permeability (default 1 everywhere) + + Returns: + J distribution for the unidirectional source + """ + E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis, + polarity=polarity, slices=slices) + + smask = [slice(None)] * 4 + if polarity > 0: + smask[axis + 1] = slice(slices[axis].start, None) + else: + smask[axis + 1] = slice(None, slices[axis].stop) + + mask = numpy.zeros_like(E_expanded, dtype=int) + mask[tuple(smask)] = 1 + + masked_e2j = operators.e_boundary_source(mask=vec(mask), omega=omega, dxes=dxes, epsilon=vec(epsilon), mu=vec(mu)) + J = unvec(masked_e2j @ vec(E_expanded), E.shape[1:]) + return J + + +def compute_overlap_e(E: field_t, + wavenumber: complex, + dxes: dx_lists_t, + axis: int, + polarity: int, + slices: List[slice], + ) -> field_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) + [assumes reflection symmetry]. + + Args: + E: E-field of the mode + 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` + 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 + as the waveguide cross-section. slices[axis] should select only one item. + mu: Magnetic permeability (default 1 everywhere) + + Returns: + overlap_e such that `numpy.sum(overlap_e * other_e)` computes the overlap integral + """ + slices = tuple(slices) + + Ee = expand_e(E=E, wavenumber=wavenumber, dxes=dxes, + axis=axis, polarity=polarity, slices=slices) + + start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity)) + + slices2 = list(slices) + slices2[axis] = slice(start, stop) + slices2 = (slice(None), *slices2) + + Etgt = numpy.zeros_like(Ee) + Etgt[slices2] = Ee[slices2] + + Etgt /= (Etgt.conj() * Etgt).sum() + return Etgt + + +def expand_e(E: field_t, + wavenumber: complex, + dxes: dx_lists_t, + axis: int, + polarity: int, + slices: List[slice], + ) -> field_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 + axis). This assumes the epsilon cross-section remains constant throughout the + entire domain; it is up to the caller to truncate the expansion to any regions + where it is valid. + + Args: + E: E-field of the mode + wavenumber: Wavenumber of the mode + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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 + as the waveguide cross-section. slices[axis] should select only one item. + + Returns: + `E`, with the original field expanded along the specified `axis`. + """ + slices = tuple(slices) + + # Determine phase factors for parallel slices + a_shape = numpy.roll([1, -1, 1, 1], axis) + a_E = numpy.real(dxes[0][axis]).cumsum() + r_E = a_E - a_E[slices[axis]] + iphi = polarity * -1j * wavenumber + phase_E = numpy.exp(iphi * r_E).reshape(a_shape) + + # Expand our slice to the entire grid using the phase factors + E_expanded = numpy.zeros_like(E) + + slices_exp = list(slices) + slices_exp[axis] = slice(E.shape[axis + 1]) + slices_exp = (slice(None), *slices_exp) + + slices_in = (slice(None), *slices) + + E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in] + return E_expanded diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py new file mode 100644 index 0000000..ebfb41d --- /dev/null +++ b/meanas/fdfd/waveguide_cyl.py @@ -0,0 +1,138 @@ +""" +Operators and helper functions for cylindrical waveguides with unchanging cross-section. + +WORK IN PROGRESS, CURRENTLY BROKEN + +As the z-dependence is known, all the functions in this file assume a 2D grid + (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`). +""" +# TODO update module docs + +from typing import List, Tuple, Dict +import numpy +from numpy.linalg import norm +import scipy.sparse as sparse + +from .. import vec, unvec, dx_lists_t, field_t, vfield_t +from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration +from . import operators + + +__author__ = 'Jan Petykiewicz' + + +def cylindrical_operator(omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + r0: float, + ) -> sparse.spmatrix: + """ + Cylindrical coordinate waveguide operator of the form + + TODO + + for use with a field vector of the form `[E_r, E_y]`. + + This operator can be used to form an eigenvalue problem of the form + A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y] + + which can then be solved for the eigenmodes of the system + (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields). + + Args: + omega: The angular frequency of the system + dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.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. + + Returns: + Sparse matrix representation of the operator + """ + + Dfx, Dfy = operators.deriv_forward(dxes[0]) + Dbx, Dby = operators.deriv_back(dxes[1]) + + rx = r0 + numpy.cumsum(dxes[0][0]) + ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0]) + tx = rx/r0 + ty = ry/r0 + + Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1))) + Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1))) + + eps_parts = numpy.split(epsilon, 3) + eps_x = sparse.diags(eps_parts[0]) + eps_y = sparse.diags(eps_parts[1]) + eps_z_inv = sparse.diags(1 / eps_parts[2]) + + pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby)) + pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx)) + a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy + a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx + b0 = Dbx @ Ty @ Dfy + b1 = Dby @ Ty @ Dfx + + diag = sparse.block_diag + op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \ + - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1)) + + return op + + +def solve_mode(mode_number: int, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + r0: float, + ) -> Dict[str, complex or field_t]: + """ + TODO: fixup + Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode + of the bent waveguide with the specified mode number. + + 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. + 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 + r within the simulation domain. + + Returns: + `{'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}` + """ + + ''' + Solve for the largest-magnitude eigenvalue of the real operator + ''' + dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] + + A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0) + eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3) + e_xy = eigvecs[:, -(mode_number+1)] + + ''' + Now solve for the eigenvector of the full operator, using the real operator's + eigenvector as an initial guess for Rayleigh quotient iteration. + ''' + A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0) + eigval, e_xy = rayleigh_quotient_iteration(A, e_xy) + + # Calculate the wave-vector (force the real part to be positive) + wavenumber = numpy.sqrt(eigval) + wavenumber *= numpy.sign(numpy.real(wavenumber)) + + # TODO: Perform correction on wavenumber to account for numerical dispersion. + + shape = [d.size for d in dxes[0]] + e_xy = numpy.hstack((e_xy, numpy.zeros(shape[0] * shape[1]))) + fields = { + 'wavenumber': wavenumber, + 'E': unvec(e_xy, shape), +# 'E': unvec(e, shape), +# 'H': unvec(h, shape), + } + + return fields diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py deleted file mode 100644 index 4318cc0..0000000 --- a/meanas/fdfd/waveguide_mode.py +++ /dev/null @@ -1,307 +0,0 @@ -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 . import operators, waveguide, functional -from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration - - -def vsolve_waveguide_mode_2d(mode_number: int, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - mode_margin: int = 2, - ) -> Tuple[vfield_t, complex]: - """ - Given a 2d region, attempts to solve for the eigenmode with the specified mode number. - - :param mode_number: Number of the mode, 0-indexed. - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param epsilon: Dielectric constant - :param mu: Magnetic permeability (default 1 everywhere) - :param mode_margin: The eigensolver will actually solve for (mode_number + mode_margin) - modes, but only return the target mode. Increasing this value can improve the solver's - ability to find the correct mode. Default 2. - :return: (e_xy, wavenumber) - """ - - ''' - 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)) - - eigvals, eigvecs = signed_eigensolve(A_r, mode_number + mode_margin) - e_xy = eigvecs[:, -(mode_number + 1)] - - ''' - Now solve for the eigenvector of the full operator, using the real operator's - eigenvector as an initial guess for Rayleigh quotient iteration. - ''' - A = waveguide.operator_e(omega, dxes, epsilon, mu) - eigval, e_xy = rayleigh_quotient_iteration(A, e_xy) - - # Calculate the wave-vector (force the real part to be positive) - wavenumber = numpy.sqrt(eigval) - wavenumber *= numpy.sign(numpy.real(wavenumber)) - - return e_xy, wavenumber - - - -def solve_waveguide_mode(mode_number: int, - omega: complex, - dxes: dx_lists_t, - axis: int, - polarity: int, - slices: List[slice], - epsilon: field_t, - mu: field_t = None, - ) -> Dict[str, complex or numpy.ndarray]: - """ - Given a 3D grid, selects a slice from the grid and attempts to - solve for an eigenmode propagating through that slice. - - :param mode_number: Number of the mode, 0-indexed - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param axis: Propagation axis (0=x, 1=y, 2=z) - :param polarity: Propagation direction (+1 for +ve, -1 for -ve) - :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use - as the waveguide cross-section. slices[axis] should select only one - :param epsilon: Dielectric constant - :param mu: Magnetic permeability (default 1 everywhere) - :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} - """ - if mu is None: - mu = numpy.ones_like(epsilon) - - slices = tuple(slices) - - ''' - Solve the 2D problem in the specified plane - ''' - # Define rotation to set z as propagation direction - order = numpy.roll(range(3), 2 - axis) - reverse_order = numpy.roll(range(3), axis - 2) - - # Find dx in propagation direction - dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes]) - dx_prop = 0.5 * sum(dxab_forward)[0] - - # Reduce to 2D and solve the 2D problem - args_2d = { - 'omega': omega, - 'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes], - 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]), - 'mu': vec([mu[i][slices].transpose(order) for i in order]), - } - e_xy, wavenumber_2d = vsolve_waveguide_mode_2d(mode_number, **args_2d) - - ''' - Apply corrections and expand to 3D - ''' - # Correct wavenumber to account for numerical dispersion. - wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2) - print(wavenumber_2d / wavenumber) - - shape = [d.size for d in args_2d['dxes'][0]] - ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber) - e = unvec(ve, shape) - h = unvec(vh, shape) - - # Adjust for propagation direction - h *= polarity - - # Apply phase shift to H-field - h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop) - e[2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop) - - # Expand E, H to full epsilon space we were given - E = numpy.zeros_like(epsilon, dtype=complex) - H = numpy.zeros_like(epsilon, dtype=complex) - for a, o in enumerate(reverse_order): - E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order) - H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order) - - results = { - 'wavenumber': wavenumber, - 'wavenumber_2d': wavenumber_2d, - 'H': H, - 'E': E, - } - return results - - -def compute_source(E: field_t, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - axis: int, - polarity: int, - slices: List[slice], - epsilon: field_t, - mu: field_t = None, - ) -> field_t: - """ - Given an eigenmode obtained by solve_waveguide_mode, returns the current source distribution - necessary to position a unidirectional source at the slice location. - - :param E: E-field of the mode - :param wavenumber: Wavenumber of the mode - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param axis: Propagation axis (0=x, 1=y, 2=z) - :param polarity: Propagation direction (+1 for +ve, -1 for -ve) - :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use - as the waveguide cross-section. slices[axis] should select only one - :param mu: Magnetic permeability (default 1 everywhere) - :return: J distribution for the unidirectional source - """ - E_expanded = expand_wgmode_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis, - polarity=polarity, slices=slices) - - smask = [slice(None)] * 4 - if polarity > 0: - smask[axis + 1] = slice(slices[axis].start, None) - else: - smask[axis + 1] = slice(None, slices[axis].stop) - - mask = numpy.zeros_like(E_expanded, dtype=int) - mask[tuple(smask)] = 1 - - masked_e2j = operators.e_boundary_source(mask=vec(mask), omega=omega, dxes=dxes, epsilon=vec(epsilon), mu=vec(mu)) - J = unvec(masked_e2j @ vec(E_expanded), E.shape[1:]) - return J - - -def compute_overlap_e(E: field_t, - wavenumber: complex, - dxes: dx_lists_t, - axis: int, - polarity: int, - slices: List[slice], - ) -> field_t: # TODO DOCS - """ - Given an eigenmode obtained by solve_waveguide_mode, calculates overlap_e for the - mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn) - [assumes reflection symmetry].i - - overlap_e makes use of the e2h operator to collapse the above expression into - (vec(E) @ vec(overlap_e)), allowing for simple calculation of the mode overlap. - - :param E: E-field of the mode - :param H: H-field of the mode (advanced by half of a Yee cell from E) - :param wavenumber: Wavenumber of the mode - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :param axis: Propagation axis (0=x, 1=y, 2=z) - :param polarity: Propagation direction (+1 for +ve, -1 for -ve) - :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use - as the waveguide cross-section. slices[axis] should select only one - :param mu: Magnetic permeability (default 1 everywhere) - :return: overlap_e for calculating the mode overlap - """ - slices = tuple(slices) - - Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes, - axis=axis, polarity=polarity, slices=slices) - - start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity)) - - slices2 = list(slices) - slices2[axis] = slice(start, stop) - slices2 = (slice(None), *slices2) - - Etgt = numpy.zeros_like(Ee) - Etgt[slices2] = Ee[slices2] - - Etgt /= (Etgt.conj() * Etgt).sum() - return Etgt - - -def solve_waveguide_mode_cylindrical(mode_number: int, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - r0: float, - ) -> Dict[str, complex or field_t]: - """ - TODO: fixup - Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode - of the bent waveguide with the specified mode number. - - :param mode_number: Number of the mode, 0-indexed - :param omega: Angular frequency of the simulation - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types. - The first coordinate is assumed to be r, the second is y. - :param epsilon: Dielectric constant - :param r0: Radius of curvature for the simulation. This should be the minimum value of - r within the simulation domain. - :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} - """ - - ''' - Solve for the largest-magnitude eigenvalue of the real operator - ''' - dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] - - A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0) - eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3) - e_xy = eigvecs[:, -(mode_number+1)] - - ''' - Now solve for the eigenvector of the full operator, using the real operator's - eigenvector as an initial guess for Rayleigh quotient iteration. - ''' - A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0) - eigval, e_xy = rayleigh_quotient_iteration(A, e_xy) - - # Calculate the wave-vector (force the real part to be positive) - wavenumber = numpy.sqrt(eigval) - wavenumber *= numpy.sign(numpy.real(wavenumber)) - - # TODO: Perform correction on wavenumber to account for numerical dispersion. - - shape = [d.size for d in dxes[0]] - e_xy = numpy.hstack((e_xy, numpy.zeros(shape[0] * shape[1]))) - fields = { - 'wavenumber': wavenumber, - 'E': unvec(e_xy, shape), -# 'E': unvec(e, shape), -# 'H': unvec(h, shape), - } - - return fields - - -def expand_wgmode_e(E: field_t, - wavenumber: complex, - dxes: dx_lists_t, - axis: int, - polarity: int, - slices: List[slice], - ) -> field_t: - slices = tuple(slices) - - # Determine phase factors for parallel slices - a_shape = numpy.roll([1, -1, 1, 1], axis) - a_E = numpy.real(dxes[0][axis]).cumsum() - r_E = a_E - a_E[slices[axis]] - iphi = polarity * -1j * wavenumber - phase_E = numpy.exp(iphi * r_E).reshape(a_shape) - - # Expand our slice to the entire grid using the phase factors - E_expanded = numpy.zeros_like(E) - - slices_exp = list(slices) - slices_exp[axis] = slice(E.shape[axis + 1]) - slices_exp = (slice(None), *slices_exp) - - slices_in = (slice(None), *slices) - - E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in] - return E_expanded diff --git a/meanas/fdmath/functional.py b/meanas/fdmath/functional.py new file mode 100644 index 0000000..8593d4a --- /dev/null +++ b/meanas/fdmath/functional.py @@ -0,0 +1,109 @@ +""" +Math functions for finite difference simulations + +Basic discrete calculus etc. +""" +from typing import List, Callable, Tuple, Dict +import numpy + +from .. import field_t, field_updater + + +def deriv_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: + """ + Utility operators for taking discretized derivatives (backward variant). + + Args: + dx_e: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + List of functions for taking forward derivatives along each axis. + """ + if dx_e: + derivs = [lambda f: (numpy.roll(f, -1, axis=0) - f) / dx_e[0][:, None, None], + lambda f: (numpy.roll(f, -1, axis=1) - f) / dx_e[1][None, :, None], + lambda f: (numpy.roll(f, -1, axis=2) - f) / dx_e[2][None, None, :]] + else: + derivs = [lambda f: numpy.roll(f, -1, axis=0) - f, + lambda f: numpy.roll(f, -1, axis=1) - f, + lambda f: numpy.roll(f, -1, axis=2) - f] + return derivs + + +def deriv_back(dx_h: List[numpy.ndarray] = None) -> field_updater: + """ + Utility operators for taking discretized derivatives (forward variant). + + Args: + dx_h: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + List of functions for taking forward derivatives along each axis. + """ + if dx_h: + derivs = [lambda f: (f - numpy.roll(f, 1, axis=0)) / dx_h[0][:, None, None], + lambda f: (f - numpy.roll(f, 1, axis=1)) / dx_h[1][None, :, None], + lambda f: (f - numpy.roll(f, 1, axis=2)) / dx_h[2][None, None, :]] + else: + derivs = [lambda f: f - numpy.roll(f, 1, axis=0), + lambda f: f - numpy.roll(f, 1, axis=1), + lambda f: f - numpy.roll(f, 1, axis=2)] + return derivs + + +def curl_forward(dx_e: List[numpy.ndarray] = None) -> field_updater: + """ + Curl operator for use with the E field. + + Args: + dx_e: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + Function `f` for taking the discrete forward curl of a field, + `f(E)` -> curlE \\( = \\nabla_f \\times E \\) + """ + Dx, Dy, Dz = deriv_forward(dx_e) + + def ce_fun(e: field_t) -> field_t: + output = numpy.empty_like(e) + output[0] = Dy(e[2]) + output[1] = Dz(e[0]) + output[2] = Dx(e[1]) + output[0] -= Dz(e[1]) + output[1] -= Dx(e[2]) + output[2] -= Dy(e[0]) + return output + + return ce_fun + + +def curl_back(dx_h: List[numpy.ndarray] = None) -> field_updater: + """ + Create a function which takes the backward curl of a field. + + Args: + dx_h: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + Function `f` for taking the discrete backward curl of a field, + `f(H)` -> curlH \\( = \\nabla_b \\times H \\) + """ + Dx, Dy, Dz = deriv_back(dx_h) + + def ch_fun(h: field_t) -> field_t: + output = numpy.empty_like(h) + output[0] = Dy(h[2]) + output[1] = Dz(h[0]) + output[2] = Dx(h[1]) + output[0] -= Dz(h[1]) + output[1] -= Dx(h[2]) + output[2] -= Dy(h[0]) + return output + + return ch_fun + + diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py new file mode 100644 index 0000000..64f04aa --- /dev/null +++ b/meanas/fdmath/operators.py @@ -0,0 +1,231 @@ +""" +Matrix operators for finite difference simulations + +Basic discrete calculus etc. +""" +from typing import List, Callable, Tuple, Dict +import numpy +import scipy.sparse as sparse + +from .. import field_t, vfield_t + + +def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: + """ + Utility operator for performing a circular shift along a specified axis by a + specified number of elements. + + Args: + axis: Axis to shift along. x=0, y=1, z=2 + shape: Shape of the grid being shifted + shift_distance: Number of cells to shift by. May be negative. Default 1. + + Returns: + Sparse matrix for performing the circular shift. + """ + if len(shape) not in (2, 3): + raise Exception('Invalid shape: {}'.format(shape)) + if axis not in range(len(shape)): + raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape)) + + shifts = [abs(shift_distance) if a == axis else 0 for a in range(3)] + shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)] + ijk = numpy.meshgrid(*shifted_diags, indexing='ij') + + n = numpy.prod(shape) + i_ind = numpy.arange(n) + j_ind = numpy.ravel_multi_index(ijk, shape, order='C') + + vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) + + d = sparse.csr_matrix(vij, shape=(n, n)) + + if shift_distance < 0: + d = d.T + + return d + + +def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: + """ + Utility operator for performing an n-element shift along a specified axis, with mirror + boundary conditions applied to the cells beyond the receding edge. + + Args: + axis: Axis to shift along. x=0, y=1, z=2 + shape: Shape of the grid being shifted + shift_distance: Number of cells to shift by. May be negative. Default 1. + + Returns: + Sparse matrix for performing the shift-with-mirror. + """ + if len(shape) not in (2, 3): + raise Exception('Invalid shape: {}'.format(shape)) + if axis not in range(len(shape)): + raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape)) + if shift_distance >= shape[axis]: + raise Exception('Shift ({}) is too large for axis {} of size {}'.format( + shift_distance, axis, shape[axis])) + + def mirrored_range(n, s): + v = numpy.arange(n) + s + v = numpy.where(v >= n, 2 * n - v - 1, v) + v = numpy.where(v < 0, - 1 - v, v) + return v + + shifts = [shift_distance if a == axis else 0 for a in range(3)] + shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts)] + ijk = numpy.meshgrid(*shifted_diags, indexing='ij') + + n = numpy.prod(shape) + i_ind = numpy.arange(n) + j_ind = numpy.ravel_multi_index(ijk, shape, order='C') + + vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) + + d = sparse.csr_matrix(vij, shape=(n, n)) + return d + + +def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]: + """ + Utility operators for taking discretized derivatives (forward variant). + + Args: + dx_e: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + List of operators for taking forward derivatives along each axis. + """ + shape = [s.size for s in dx_e] + n = numpy.prod(shape) + + dx_e_expanded = numpy.meshgrid(*dx_e, indexing='ij') + + def deriv(axis): + return rotation(axis, shape, 1) - sparse.eye(n) + + Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a) + for a, dx in enumerate(dx_e_expanded)] + + return Ds + + +def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]: + """ + Utility operators for taking discretized derivatives (backward variant). + + Args: + dx_h: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + List of operators for taking forward derivatives along each axis. + """ + shape = [s.size for s in dx_h] + n = numpy.prod(shape) + + dx_h_expanded = numpy.meshgrid(*dx_h, indexing='ij') + + def deriv(axis): + return rotation(axis, shape, -1) - sparse.eye(n) + + Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a) + for a, dx in enumerate(dx_h_expanded)] + + return Ds + + +def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix: + """ + Cross product operator + + Args: + B: List `[Bx, By, Bz]` of sparse matrices corresponding to the x, y, z + portions of the operator on the left side of the cross product. + + Returns: + Sparse matrix corresponding to (B x), where x is the cross product. + """ + n = B[0].shape[0] + zero = sparse.csr_matrix((n, n)) + return sparse.bmat([[zero, -B[2], B[1]], + [B[2], zero, -B[0]], + [-B[1], B[0], zero]]) + + +def vec_cross(b: vfield_t) -> sparse.spmatrix: + """ + Vector cross product operator + + Args: + b: Vector on the left side of the cross product. + + Returns: + + Sparse matrix corresponding to (b x), where x is the cross product. + + """ + B = [sparse.diags(c) for c in numpy.split(b, 3)] + return cross(B) + + +def avg_forward(axis: int, shape: List[int]) -> sparse.spmatrix: + """ + Forward average operator `(x4 = (x4 + x5) / 2)` + + Args: + axis: Axis to average along (x=0, y=1, z=2) + shape: Shape of the grid to average + + Returns: + Sparse matrix for forward average operation. + """ + if len(shape) not in (2, 3): + raise Exception('Invalid shape: {}'.format(shape)) + + n = numpy.prod(shape) + return 0.5 * (sparse.eye(n) + rotation(axis, shape)) + + +def avg_back(axis: int, shape: List[int]) -> sparse.spmatrix: + """ + Backward average operator `(x4 = (x4 + x3) / 2)` + + Args: + axis: Axis to average along (x=0, y=1, z=2) + shape: Shape of the grid to average + + Returns: + Sparse matrix for backward average operation. + """ + return avg_forward(axis, shape).T + + +def curl_forward(dx_e: List[numpy.ndarray]) -> sparse.spmatrix: + """ + Curl operator for use with the E field. + + Args: + dx_e: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + Sparse matrix for taking the discretized curl of the E-field + """ + return cross(deriv_forward(dx_e)) + + +def curl_back(dx_h: List[numpy.ndarray]) -> sparse.spmatrix: + """ + Curl operator for use with the H field. + + Args: + dx_h: Lists of cell sizes for all axes + `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`. + + Returns: + Sparse matrix for taking the discretized curl of the H-field + """ + return cross(deriv_back(dx_h)) diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py index a1d278a..02cfc05 100644 --- a/meanas/fdtd/__init__.py +++ b/meanas/fdtd/__init__.py @@ -1,5 +1,5 @@ """ -Basic FDTD functionality +Utilities for running finite-difference time-domain (FDTD) simulations """ from .base import maxwell_e, maxwell_h diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py index 389a7b5..87a234b 100644 --- a/meanas/fdtd/base.py +++ b/meanas/fdtd/base.py @@ -5,70 +5,14 @@ from typing import List, Callable, Tuple, Dict import numpy from .. import dx_lists_t, field_t, field_updater +from ..fdmath.functional import curl_forward, curl_back + __author__ = 'Jan Petykiewicz' -def curl_h(dxes: dx_lists_t = None) -> field_updater: - """ - Curl operator for use with the H field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Function for taking the discretized curl of the H-field, F(H) -> curlH - """ - if dxes: - dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij') - - def dh(f, ax): - return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax] - else: - def dh(f, ax): - return f - numpy.roll(f, 1, axis=ax) - - def ch_fun(h: field_t) -> field_t: - output = numpy.empty_like(h) - output[0] = dh(h[2], 1) - output[1] = dh(h[0], 2) - output[2] = dh(h[1], 0) - output[0] -= dh(h[1], 2) - output[1] -= dh(h[2], 0) - output[2] -= dh(h[0], 1) - return output - - return ch_fun - - -def curl_e(dxes: dx_lists_t = None) -> field_updater: - """ - Curl operator for use with the E field. - - :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types - :return: Function for taking the discretized curl of the E-field, F(E) -> curlE - """ - if dxes is not None: - dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij') - - def de(f, ax): - return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax] - else: - def de(f, ax): - return numpy.roll(f, -1, axis=ax) - f - - def ce_fun(e: field_t) -> field_t: - output = numpy.empty_like(e) - output[0] = de(e[2], 1) - output[1] = de(e[0], 2) - output[2] = de(e[1], 0) - output[0] -= de(e[1], 2) - output[1] -= de(e[2], 0) - output[2] -= de(e[0], 1) - return output - - return ce_fun - - def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater: - curl_h_fun = curl_h(dxes) + curl_h_fun = curl_back(dxes[1]) def me_fun(e: field_t, h: field_t, epsilon: field_t): e += dt * curl_h_fun(h) / epsilon @@ -78,7 +22,7 @@ def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater: def maxwell_h(dt: float, dxes: dx_lists_t = None) -> field_updater: - curl_e_fun = curl_e(dxes) + curl_e_fun = curl_forward(dxes[0]) def mh_fun(e: field_t, h: field_t): h -= dt * curl_e_fun(e) diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py index 8644646..d5ea1dc 100644 --- a/meanas/fdtd/energy.py +++ b/meanas/fdtd/energy.py @@ -35,6 +35,7 @@ def poynting_divergence(s: field_t = None, if s is None: s = poynting(e, h, dxes=dxes) + #TODO use deriv operators ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) + (s[1] - numpy.roll(s[1], 1, axis=1)) + (s[2] - numpy.roll(s[2], 1, axis=2))) diff --git a/meanas/test/__init__.py b/meanas/test/__init__.py index e69de29..e02b636 100644 --- a/meanas/test/__init__.py +++ b/meanas/test/__init__.py @@ -0,0 +1,3 @@ +""" +Tests (run with `python3 -m pytest -rxPXs | tee results.txt`) +""" diff --git a/meanas/test/test_fdfd_pml.py b/meanas/test/test_fdfd_pml.py index 112193d..eafa14f 100644 --- a/meanas/test/test_fdfd_pml.py +++ b/meanas/test/test_fdfd_pml.py @@ -101,8 +101,8 @@ def j_distribution(request, shape, epsilon, dxes, omega, src_polarity): slices[dim] = slice(shape[dim + 1] // 2, shape[dim + 1] // 2 + 1) - j = fdfd.waveguide_mode.compute_source(E=e, wavenumber=wavenumber_corrected, omega=omega, dxes=dxes, - axis=dim, polarity=src_polarity, slices=slices, epsilon=epsilon) + j = fdfd.waveguide_3d.compute_source(E=e, wavenumber=wavenumber_corrected, omega=omega, dxes=dxes, + axis=dim, polarity=src_polarity, slices=slices, epsilon=epsilon) yield j @@ -145,4 +145,3 @@ def sim(request, shape, epsilon, dxes, j_distribution, omega, pec, pmc): ) return sim - diff --git a/meanas/types.py b/meanas/types.py index 7dc5c1c..667a286 100644 --- a/meanas/types.py +++ b/meanas/types.py @@ -6,17 +6,20 @@ 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) +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)""" + +dx_lists_t = List[List[numpy.ndarray]] ''' '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_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/pdoc_templates/config.mako b/pdoc_templates/config.mako new file mode 100644 index 0000000..0642566 --- /dev/null +++ b/pdoc_templates/config.mako @@ -0,0 +1,46 @@ +<%! + # Template configuration. Copy over in your template directory + # (used with --template-dir) and adapt as required. + html_lang = 'en' + show_inherited_members = False + extract_module_toc_into_sidebar = True + list_class_variables_in_index = True + sort_identifiers = True + show_type_annotations = True + + # Show collapsed source code block next to each item. + # Disabling this can improve rendering speed of large modules. + show_source_code = True + + # If set, format links to objects in online source code repository + # according to this template. Supported keywords for interpolation + # are: commit, path, start_line, end_line. + #git_link_template = 'https://github.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}' + #git_link_template = 'https://gitlab.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}' + #git_link_template = 'https://bitbucket.org/USER/PROJECT/src/{commit}/{path}#lines-{start_line}:{end_line}' + #git_link_template = 'https://CGIT_HOSTNAME/PROJECT/tree/{path}?id={commit}#n{start-line}' + git_link_template = None + + # A prefix to use for every HTML hyperlink in the generated documentation. + # No prefix results in all links being relative. + link_prefix = '' + + # Enable syntax highlighting for code/source blocks by including Highlight.js + syntax_highlighting = True + + # Set the style keyword such as 'atom-one-light' or 'github-gist' + # Options: https://github.com/highlightjs/highlight.js/tree/master/src/styles + # Demo: https://highlightjs.org/static/demo/ + hljs_style = 'github' + + # If set, insert Google Analytics tracking code. Value is GA + # tracking id (UA-XXXXXX-Y). + google_analytics = '' + + # If set, render LaTeX math syntax within \(...\) (inline equations), + # or within \[...\] or $$...$$ or `.. math::` (block equations) + # as nicely-formatted math formulas using MathJax. + # Note: in Python docstrings, either all backslashes need to be escaped (\\) + # or you need to use raw r-strings. + latex_math = True +%> diff --git a/pdoc_templates/css.mako b/pdoc_templates/css.mako new file mode 100644 index 0000000..39a77ed --- /dev/null +++ b/pdoc_templates/css.mako @@ -0,0 +1,389 @@ +<%! + from pdoc.html_helpers import minify_css +%> + +<%def name="mobile()" filter="minify_css"> + .flex { + display: flex !important; + } + + body { + line-height: 1.5em; + background: black; + color: #DDD; + } + + #content { + padding: 20px; + } + + #sidebar { + padding: 30px; + overflow: hidden; + } + + .http-server-breadcrumbs { + font-size: 130%; + margin: 0 0 15px 0; + } + + #footer { + font-size: .75em; + padding: 5px 30px; + border-top: 1px solid #ddd; + text-align: right; + } + #footer p { + margin: 0 0 0 1em; + display: inline-block; + } + #footer p:last-child { + margin-right: 30px; + } + + h1, h2, h3, h4, h5 { + font-weight: 300; + } + h1 { + font-size: 2.5em; + line-height: 1.1em; + } + h2 { + font-size: 1.75em; + margin: 1em 0 .50em 0; + } + h3 { + font-size: 1.4em; + margin: 25px 0 10px 0; + } + h4 { + margin: 0; + font-size: 105%; + } + + a { + color: #999; + text-decoration: none; + transition: color .3s ease-in-out; + } + a:hover { + color: #18d; + } + + .title code { + font-weight: bold; + } + h2[id^="header-"] { + margin-top: 2em; + } + .ident { + color: #7ff; + } + + pre code { + background: transparent; + font-size: .8em; + line-height: 1.4em; + } + code { + background: #0d0d0e; + padding: 1px 4px; + overflow-wrap: break-word; + } + h1 code { background: transparent } + + pre { + background: #111; + border: 0; + border-top: 1px solid #ccc; + border-bottom: 1px solid #ccc; + margin: 1em 0; + padding: 1ex; + } + + #http-server-module-list { + display: flex; + flex-flow: column; + } + #http-server-module-list div { + display: flex; + } + #http-server-module-list dt { + min-width: 10%; + } + #http-server-module-list p { + margin-top: 0; + } + + .toc ul, + #index { + list-style-type: none; + margin: 0; + padding: 0; + } + #index code { + background: transparent; + } + #index h3 { + border-bottom: 1px solid #ddd; + } + #index ul { + padding: 0; + } + #index h4 { + font-weight: bold; + } + #index h4 + ul { + margin-bottom:.6em; + } + /* Make TOC lists have 2+ columns when viewport is wide enough. + Assuming ~20-character identifiers and ~30% wide sidebar. */ + @media (min-width: 200ex) { #index .two-column { column-count: 2 } } + @media (min-width: 300ex) { #index .two-column { column-count: 3 } } + + dl { + margin-bottom: 2em; + } + dl dl:last-child { + margin-bottom: 4em; + } + dd { + margin: 0 0 1em 3em; + } + #header-classes + dl > dd { + margin-bottom: 3em; + } + dd dd { + margin-left: 2em; + } + dd p { + margin: 10px 0; + } + .name { + background: #111; + font-weight: bold; + font-size: .85em; + padding: 5px 10px; + display: inline-block; + min-width: 40%; + } + .name:hover { + background: #101010; + } + .name > span:first-child { + white-space: nowrap; + } + .name.class > span:nth-child(2) { + margin-left: .4em; + } + .inherited { + color: #777; + border-left: 5px solid #eee; + padding-left: 1em; + } + .inheritance em { + font-style: normal; + font-weight: bold; + } + + /* Docstrings titles, e.g. in numpydoc format */ + .desc h2 { + font-weight: 400; + font-size: 1.25em; + } + .desc h3 { + font-size: 1em; + } + .desc dt code { + background: inherit; /* Don't grey-back parameters */ + } + + .source summary, + .git-link-div { + color: #aaa; + text-align: right; + font-weight: 400; + font-size: .8em; + text-transform: uppercase; + } + .source summary > * { + white-space: nowrap; + cursor: pointer; + } + .git-link { + color: inherit; + margin-left: 1em; + } + .source pre { + max-height: 500px; + overflow: auto; + margin: 0; + } + .source pre code { + font-size: 12px; + overflow: visible; + } + .hlist { + list-style: none; + } + .hlist li { + display: inline; + } + .hlist li:after { + content: ',\2002'; + } + .hlist li:last-child:after { + content: none; + } + .hlist .hlist { + display: inline; + padding-left: 1em; + } + + img { + max-width: 100%; + } + + .admonition { + padding: .1em .5em; + margin-bottom: 1em; + } + .admonition-title { + font-weight: bold; + } + .admonition.note, + .admonition.info, + .admonition.important { + background: #610; + } + .admonition.todo, + .admonition.versionadded, + .admonition.tip, + .admonition.hint { + background: #202; + } + .admonition.warning, + .admonition.versionchanged, + .admonition.deprecated { + background: #02b; + } + .admonition.error, + .admonition.danger, + .admonition.caution { + background: darkpink; + } + + +<%def name="desktop()" filter="minify_css"> + @media screen and (min-width: 700px) { + #sidebar { + width: 30%; + } + #content { + width: 70%; + max-width: 100ch; + padding: 3em 4em; + border-left: 1px solid #ddd; + } + pre code { + font-size: 1em; + } + .item .name { + font-size: 1em; + } + main { + display: flex; + flex-direction: row-reverse; + justify-content: flex-end; + } + .toc ul ul, + #index ul { + padding-left: 1.5em; + } + .toc > ul > li { + margin-top: .5em; + } + } + + +<%def name="print()" filter="minify_css"> +@media print { + #sidebar h1 { + page-break-before: always; + } + .source { + display: none; + } +} +@media print { + * { + background: transparent !important; + color: #000 !important; /* Black prints faster: h5bp.com/s */ + box-shadow: none !important; + text-shadow: none !important; + } + + a[href]:after { + content: " (" attr(href) ")"; + font-size: 90%; + } + /* Internal, documentation links, recognized by having a title, + don't need the URL explicity stated. */ + a[href][title]:after { + content: none; + } + + abbr[title]:after { + content: " (" attr(title) ")"; + } + + /* + * Don't show links for images, or javascript/internal links + */ + + .ir a:after, + a[href^="javascript:"]:after, + a[href^="#"]:after { + content: ""; + } + + pre, + blockquote { + border: 1px solid #999; + page-break-inside: avoid; + } + + thead { + display: table-header-group; /* h5bp.com/t */ + } + + tr, + img { + page-break-inside: avoid; + } + + img { + max-width: 100% !important; + } + + @page { + margin: 0.5cm; + } + + p, + h2, + h3 { + orphans: 3; + widows: 3; + } + + h1, + h2, + h3, + h4, + h5, + h6 { + page-break-after: avoid; + } +} + diff --git a/pdoc_templates/html.mako b/pdoc_templates/html.mako new file mode 100644 index 0000000..9cf1137 --- /dev/null +++ b/pdoc_templates/html.mako @@ -0,0 +1,421 @@ +<% + import os + + import pdoc + from pdoc.html_helpers import extract_toc, glimpse, to_html as _to_html, format_git_link + + + def link(d, name=None, fmt='{}'): + name = fmt.format(name or d.qualname + ('()' if isinstance(d, pdoc.Function) else '')) + if not isinstance(d, pdoc.Doc) or isinstance(d, pdoc.External) and not external_links: + return name + url = d.url(relative_to=module, link_prefix=link_prefix, + top_ancestor=not show_inherited_members) + return '{}'.format(d.refname, url, name) + + + def to_html(text): + return _to_html(text, module=module, link=link, latex_math=latex_math) +%> + +<%def name="ident(name)">${name} + +<%def name="show_source(d)"> + % if (show_source_code or git_link_template) and d.source and d.obj is not getattr(d.inherits, 'obj', None): + <% git_link = format_git_link(git_link_template, d) %> + % if show_source_code: +
+ + Expand source code + % if git_link: + Browse git + %endif + + 
${d.source | h}
+
+ % elif git_link: +
Browse git
+ %endif + %endif + + +<%def name="show_desc(d, short=False)"> + <% + inherits = ' inherited' if d.inherits else '' + docstring = glimpse(d.docstring) if short or inherits else d.docstring + %> + % if d.inherits: +

+ Inherited from: + % if hasattr(d.inherits, 'cls'): + ${link(d.inherits.cls)}.${link(d.inherits, d.name)} + % else: + ${link(d.inherits)} + % endif +

+ % endif +
${docstring | to_html}
+ % if not isinstance(d, pdoc.Module): + ${show_source(d)} + % endif + + +<%def name="show_module_list(modules)"> +

Python module list

+ +% if not modules: +

No modules found.

+% else: +
+ % for name, desc in modules: +
+
${name}
+
${desc | glimpse, to_html}
+
+ % endfor +
+% endif + + +<%def name="show_column_list(items)"> + <% + two_column = len(items) >= 6 and all(len(i.name) < 20 for i in items) + %> + + + +<%def name="show_module(module)"> + <% + variables = module.variables(sort=sort_identifiers) + classes = module.classes(sort=sort_identifiers) + functions = module.functions(sort=sort_identifiers) + submodules = module.submodules() + %> + + <%def name="show_func(f)"> +
+ <% + params = ', '.join(f.params(annotate=show_type_annotations, link=link)) + returns = show_type_annotations and f.return_annotation(link=link) or '' + if returns: + returns = ' ->\N{NBSP}' + returns + %> + ${f.funcdef()} ${ident(f.name)}(${params})${returns} +
+
${show_desc(f)}
+ + +
+ % if http_server: + + % endif +

${'Namespace' if module.is_namespace else 'Module'} ${module.name}

+
+ +
+ ${module.docstring | to_html} + ${show_source(module)} +
+ +
+ % if submodules: +

Sub-modules

+
+ % for m in submodules: +
${link(m)}
+
${show_desc(m, short=True)}
+ % endfor +
+ % endif +
+ +
+ % if variables: +

Global variables

+
+ % for v in variables: +
var ${ident(v.name)}
+
${show_desc(v)}
+ % endfor +
+ % endif +
+ +
+ % if functions: +

Functions

+
+ % for f in functions: + ${show_func(f)} + % endfor +
+ % endif +
+ +
+ % if classes: +

Classes

+
+ % for c in classes: + <% + class_vars = c.class_variables(show_inherited_members, sort=sort_identifiers) + smethods = c.functions(show_inherited_members, sort=sort_identifiers) + inst_vars = c.instance_variables(show_inherited_members, sort=sort_identifiers) + methods = c.methods(show_inherited_members, sort=sort_identifiers) + mro = c.mro() + subclasses = c.subclasses() + params = ', '.join(c.params(annotate=show_type_annotations, link=link)) + %> +
+ class ${ident(c.name)} + % if params: + (${params}) + % endif +
+ +
${show_desc(c)} + + % if mro: +

Ancestors

+
    + % for cls in mro: +
  • ${link(cls)}
    • + % endfor +
+ %endif + + % if subclasses: +

Subclasses

+
    + % for sub in subclasses: +
  • ${link(sub)}
    • + % endfor +
+ % endif + % if class_vars: +

Class variables

+
+ % for v in class_vars: +
var ${ident(v.name)}
+
${show_desc(v)}
+ % endfor +
+ % endif + % if smethods: +

Static methods

+
+ % for f in smethods: + ${show_func(f)} + % endfor +
+ % endif + % if inst_vars: +

Instance variables

+
+ % for v in inst_vars: +
var ${ident(v.name)}
+
${show_desc(v)}
+ % endfor +
+ % endif + % if methods: +

Methods

+
+ % for f in methods: + ${show_func(f)} + % endfor +
+ % endif + + % if not show_inherited_members: + <% + members = c.inherited_members() + %> + % if members: +

Inherited members

+
    + % for cls, mems in members: +
  • ${link(cls)}: +
      + % for m in mems: +
    • ${link(m, name=m.name)}
      • + % endfor +
    + +
    • + % endfor +
+ % endif + % endif + +
+ % endfor +
+ % endif +
+ + +<%def name="module_index(module)"> + <% + variables = module.variables(sort=sort_identifiers) + classes = module.classes(sort=sort_identifiers) + functions = module.functions(sort=sort_identifiers) + submodules = module.submodules() + supermodule = module.supermodule + %> + + + + + + + + + + +<% + module_list = 'modules' in context.keys() # Whether we're showing module list in server mode +%> + + % if module_list: + Python module list + + % else: + ${module.name} API documentation + + % endif + + + + % if syntax_highlighting: + + %endif + + <%namespace name="css" file="css.mako" /> + + + + + % if google_analytics: + + % endif + + % if latex_math: + + % endif + + <%include file="head.mako"/> + + +
+ % if module_list: +
+ ${show_module_list(modules)} +
+ % else: +
+ ${show_module(module)} +
+ ${module_index(module)} + % endif +
+ + + +% if syntax_highlighting: + + +% endif + +% if http_server and module: ## Auto-reload on file change in dev mode + +% endif + +