From 5951f2bdb12bb57bc5ddaaa0b28de5d1f96acf34 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Mon, 5 Aug 2019 00:20:06 -0700 Subject: [PATCH] various fixes and improvements --- examples/fdfd.py | 17 +++++++++----- meanas/fdfd/farfield.py | 10 ++++---- meanas/fdfd/functional.py | 2 +- meanas/fdfd/operators.py | 2 +- meanas/fdfd/scpml.py | 2 ++ meanas/fdfd/waveguide.py | 17 +++++++------- meanas/fdfd/waveguide_mode.py | 43 ++++++++--------------------------- meanas/test/test_fdtd.py | 1 + 8 files changed, 40 insertions(+), 54 deletions(-) diff --git a/examples/fdfd.py b/examples/fdfd.py index 5ee3477..c20fc3e 100644 --- a/examples/fdfd.py +++ b/examples/fdfd.py @@ -4,7 +4,7 @@ from numpy.linalg import norm import meanas from meanas import vec, unvec -from meanas.fdfd import waveguide_mode, functional, scpml +from meanas.fdfd import waveguide_mode, functional, scpml, operators from meanas.fdfd.solvers import generic as generic_solver import gridlock @@ -56,18 +56,23 @@ def test0(solver=generic_solver): dxes = [grid.dxyz, grid.autoshifted_dxyz()] for a in (0, 1, 2): for p in (-1, 1): - dxes = meanas.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, - thickness=pml_thickness) + dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, + thickness=pml_thickness) J = [numpy.zeros_like(grid.grids[0], dtype=complex) for _ in range(3)] - J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1e5 + J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1 ''' Solve! ''' + sim_args = { + 'omega': omega, + 'dxes': dxes, + 'epsilon': vec(grid.grids), + } x = solver(J=vec(J), **sim_args) - A = functional.e_full(omega, dxes, vec(grid.grids)).tocsr() + A = operators.e_full(omega, dxes, vec(grid.grids)).tocsr() b = -1j * omega * vec(J) print('Norm of the residual is ', norm(A @ x - b)) @@ -208,7 +213,7 @@ def module_available(name): if __name__ == '__main__': - # test0() + test0() if module_available('opencl_fdfd'): from opencl_fdfd import cg_solver as opencl_solver diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py index 84a04ba..faa25b5 100644 --- a/meanas/fdfd/farfield.py +++ b/meanas/fdfd/farfield.py @@ -6,9 +6,11 @@ import numpy from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift from numpy import pi +from .. import field_t -def near_to_farfield(E_near: List[numpy.ndarray], - H_near: List[numpy.ndarray], + +def near_to_farfield(E_near: field_t, + H_near: field_t, dx: float, dy: float, padded_size: List[int] = None @@ -115,8 +117,8 @@ def near_to_farfield(E_near: List[numpy.ndarray], -def far_to_nearfield(E_far: List[numpy.ndarray], - H_far: List[numpy.ndarray], +def far_to_nearfield(E_far: field_t, + H_far: field_t, dkx: float, dky: float, padded_size: List[int] = None diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py index e57fe88..dd94421 100644 --- a/meanas/fdfd/functional.py +++ b/meanas/fdfd/functional.py @@ -8,7 +8,7 @@ e.g. E = [E_x, E_y, E_z] where each component has shape (X, Y, Z) from typing import List, Callable import numpy -from . import dx_lists_t, field_t +from .. import dx_lists_t, field_t __author__ = 'Jan Petykiewicz' diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 3b2de68..774c3d9 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -32,7 +32,7 @@ from typing import List, Tuple import numpy import scipy.sparse as sparse -from . import vec, dx_lists_t, vfield_t +from .. import vec, dx_lists_t, vfield_t __author__ = 'Jan Petykiewicz' diff --git a/meanas/fdfd/scpml.py b/meanas/fdfd/scpml.py index c4091a0..897d43a 100644 --- a/meanas/fdfd/scpml.py +++ b/meanas/fdfd/scpml.py @@ -5,6 +5,8 @@ Functions for creating stretched coordinate PMLs. from typing import List, Callable import numpy +from .. import dx_lists_t + __author__ = 'Jan Petykiewicz' diff --git a/meanas/fdfd/waveguide.py b/meanas/fdfd/waveguide.py index 48a1510..91b023c 100644 --- a/meanas/fdfd/waveguide.py +++ b/meanas/fdfd/waveguide.py @@ -23,7 +23,7 @@ import numpy from numpy.linalg import norm import scipy.sparse as sparse -from . import vec, unvec, dx_lists_t, field_t, vfield_t +from .. import vec, unvec, dx_lists_t, field_t, vfield_t from . import operators @@ -82,7 +82,8 @@ def normalized_fields(v: numpy.ndarray, omega: complex, dxes: dx_lists_t, epsilon: vfield_t, - mu: vfield_t = None + mu: vfield_t = None, + dx_prop: float = 0, ) -> Tuple[vfield_t, vfield_t]: """ Given a vector v containing the vectorized H_x and H_y fields, @@ -94,6 +95,7 @@ def normalized_fields(v: numpy.ndarray, :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 = v2e(v, wavenumber, omega, dxes, epsilon, mu=mu) @@ -105,11 +107,10 @@ def normalized_fields(v: numpy.ndarray, E = unvec(e, shape) H = unvec(h, shape) - S1 = E[0] * numpy.roll(numpy.conj(H[1]), 1, axis=0) * dxes_real[0][1] * dxes_real[1][0] - S2 = E[1] * numpy.roll(numpy.conj(H[0]), 1, axis=1) * dxes_real[0][0] * dxes_real[1][1] - S = 0.25 * ((S1 + numpy.roll(S1, 1, axis=0)) - - (S2 + numpy.roll(S2, 1, axis=1))) - P = 0.5 * numpy.real(S.sum()) + phase = numpy.exp(-1j * wavenumber * dx_prop / 2) + S1 = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0] + S2 = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1] + P = numpy.real(S1.sum() - S2.sum()) assert P > 0, 'Found a mode propagating in the wrong direction! P={}'.format(P) energy = epsilon * e.conj() * e @@ -120,8 +121,6 @@ def normalized_fields(v: numpy.ndarray, # Try to break symmetry to assign a consistent sign [experimental] E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) - logger.debug('norm_angle = {}'.format(norm_angle)) - logger.debug('norm_sign = {}'.format(sign) norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle) diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py index 6e2b192..7182377 100644 --- a/meanas/fdfd/waveguide_mode.py +++ b/meanas/fdfd/waveguide_mode.py @@ -2,9 +2,9 @@ from typing import Dict, List import numpy import scipy.sparse as sparse -from . import vec, unvec, dx_lists_t, vfield_t, field_t +from .. import vec, unvec, dx_lists_t, vfield_t, field_t from . import operators, waveguide, functional -from .eigensolvers import signed_eigensolve, rayleigh_quotient_iteration +from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration def solve_waveguide_mode_2d(mode_number: int, @@ -12,7 +12,7 @@ def solve_waveguide_mode_2d(mode_number: int, dxes: dx_lists_t, epsilon: vfield_t, mu: vfield_t = None, - wavenumber_correction: bool = True, + dx_prop: float = 0, ) -> Dict[str, complex or field_t]: """ Given a 2d region, attempts to solve for the eigenmode with the specified mode number. @@ -22,8 +22,8 @@ def solve_waveguide_mode_2d(mode_number: int, :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types :param epsilon: Dielectric constant :param mu: Magnetic permeability (default 1 everywhere) - :param wavenumber_correction: Whether to correct the wavenumber to - account for numerical dispersion (default True) + :param dx_prop: The cell width in the the propagation direction, used to apply a + correction to the wavenumber. Default 0 (i.e. continuous propagation direction) :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} """ @@ -51,15 +51,9 @@ def solve_waveguide_mode_2d(mode_number: int, ''' Perform correction on wavenumber to account for numerical dispersion. - - See Numerical Dispersion in Taflove's FDTD book. - This correction term reduces the error in emitted power, but additional - error is introduced into the E_err and H_err terms. This effect becomes - more pronounced as the wavenumber increases. ''' - if wavenumber_correction: - dx_mean = (numpy.hstack(dxes[0]) + numpy.hstack(dxes[1])).mean() / 2 #TODO figure out what dx to use here - wavenumber -= 2 * numpy.sin(numpy.real(wavenumber * dx_mean / 2)) / dx_mean - numpy.real(wavenumber) + if dx_prop != 0: + wavenumber = 2 / dx_prop * numpy.sin(wavenumber * dx_prop / 2) shape = [d.size for d in dxes[0]] fields = { @@ -79,7 +73,6 @@ def solve_waveguide_mode(mode_number: int, slices: List[slice], epsilon: field_t, mu: field_t = None, - wavenumber_correction: bool = True ) -> Dict[str, complex or numpy.ndarray]: """ Given a 3D grid, selects a slice from the grid and attempts to @@ -94,8 +87,6 @@ def solve_waveguide_mode(mode_number: int, as the waveguide cross-section. slices[axis] should select only one :param epsilon: Dielectric constant :param mu: Magnetic permeability (default 1 everywhere) - :param wavenumber_correction: Whether to correct the wavenumber to - account for numerical dispersion (default True) :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} """ if mu is None: @@ -115,7 +106,7 @@ def solve_waveguide_mode(mode_number: int, 'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes], 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]), 'mu': vec([mu[i][slices].transpose(order) for i in order]), - 'wavenumber_correction': wavenumber_correction, + 'dx_prop': dxes[0][order[2]][slices[order[2]]], } fields_2d = solve_waveguide_mode_2d(mode_number, omega=omega, **args_2d) @@ -175,9 +166,6 @@ def compute_source(E: field_t, :param mu: Magnetic permeability (default 1 everywhere) :return: J distribution for the unidirectional source """ - if mu is None: - mu = numpy.ones(3) - J = numpy.zeros_like(E, dtype=complex) M = numpy.zeros_like(E, dtype=complex) @@ -275,9 +263,9 @@ def solve_waveguide_mode_cylindrical(mode_number: int, dxes: dx_lists_t, epsilon: vfield_t, r0: float, - wavenumber_correction: bool = True, ) -> Dict[str, complex or field_t]: """ + TODO: fixup Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode of the bent waveguide with the specified mode number. @@ -288,8 +276,6 @@ def solve_waveguide_mode_cylindrical(mode_number: int, :param epsilon: Dielectric constant :param r0: Radius of curvature for the simulation. This should be the minimum value of r within the simulation domain. - :param wavenumber_correction: Whether to correct the wavenumber to - account for numerical dispersion (default True) :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} """ @@ -313,16 +299,7 @@ def solve_waveguide_mode_cylindrical(mode_number: int, wavenumber = numpy.sqrt(eigval) wavenumber *= numpy.sign(numpy.real(wavenumber)) - ''' - Perform correction on wavenumber to account for numerical dispersion. - - See Numerical Dispersion in Taflove's FDTD book. - This correction term reduces the error in emitted power, but additional - error is introduced into the E_err and H_err terms. This effect becomes - more pronounced as the wavenumber increases. - ''' - if wavenumber_correction: - wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - numpy.real(wavenumber) + # TODO: Perform correction on wavenumber to account for numerical dispersion. shape = [d.size for d in dxes[0]] v = numpy.hstack((v, numpy.zeros(shape[0] * shape[1]))) diff --git a/meanas/test/test_fdtd.py b/meanas/test/test_fdtd.py index 8443afc..08b678e 100644 --- a/meanas/test/test_fdtd.py +++ b/meanas/test/test_fdtd.py @@ -9,6 +9,7 @@ from meanas import fdtd prng = numpy.random.RandomState(12345) + def assert_fields_close(a, b, *args, **kwargs): numpy.testing.assert_allclose(a, b, verbose=False, err_msg='Fields did not match:\n{}\n{}'.format(numpy.rollaxis(a, -1), numpy.rollaxis(b, -1)), *args, **kwargs)