From af8efd00eb410542dbdbd2bf3cc58a7828864cb4 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Mon, 26 Aug 2019 00:25:36 -0700 Subject: [PATCH] Add E-field versions of waveguide mode operators, rename v->e_xy or h_xy, and add ability to specify mode margin in solve_waveguide_mode_2d --- meanas/fdfd/waveguide.py | 205 ++++++++++++++++++++++++++-------- meanas/fdfd/waveguide_mode.py | 16 ++- 2 files changed, 170 insertions(+), 51 deletions(-) diff --git a/meanas/fdfd/waveguide.py b/meanas/fdfd/waveguide.py index 611cb57..63ecb98 100644 --- a/meanas/fdfd/waveguide.py +++ b/meanas/fdfd/waveguide.py @@ -31,11 +31,36 @@ from . import operators __author__ = 'Jan Petykiewicz' -def operator(omega: complex, +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 @@ -71,27 +96,27 @@ def operator(omega: complex, mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1]))) mu_z_inv = sparse.diags(1 / mu_parts[2]) - op = omega ** 2 * eps_yx @ mu_xy + \ + 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(v: numpy.ndarray, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None, - dx_prop: float = 0, - ) -> Tuple[vfield_t, vfield_t]: +def normalized_fields_e(e_xy: numpy.ndarray, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + 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, + 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 v: Vector containing H_x and H_y fields - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v + :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 @@ -99,9 +124,51 @@ def normalized_fields(v: numpy.ndarray, :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) - h = v2h(v, wavenumber, dxes, mu=mu) + 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, wavenumber=wavenumber, omega=omega, + dxes=dxes, epsilon=epsilon, mu=mu, dx_prop=dx_prop) + 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, + dx_prop: 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, wavenumber=wavenumber, omega=omega, + dxes=dxes, epsilon=epsilon, mu=mu, dx_prop=dx_prop) + return e_norm, h_norm + + +def _normalized_fields(e: numpy.ndarray, + h: numpy.ndarray, + wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None, + dx_prop: 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)] @@ -131,56 +198,104 @@ def normalized_fields(v: numpy.ndarray, return e, h -def v2h(v: numpy.ndarray, - wavenumber: complex, - dxes: dx_lists_t, - mu: vfield_t = None - ) -> vfield_t: +def exy2h(wavenumber: complex, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + mu: vfield_t = None + ) -> sparse.spmatrix: """ - Given a vector v containing the vectorized H_x and H_y fields, - returns a vectorized H including all three H components. + 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 v: Vector containing H_x and H_y fields - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v + :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: Vectorized H field with all vector components + :return: Sparse matrix representing the operator """ Dfx, Dfy = operators.deriv_forward(dxes[0]) - op = sparse.hstack((Dfx, Dfy)) + 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]) - op = mu_z_inv @ op @ mu_xy + hxy2hz = mu_z_inv @ hxy2hz @ mu_xy - w = op @ v / (1j * wavenumber) - return numpy.hstack((v, w)).flatten() + n_pts = dxes[1][0].size * dxes[1][1].size + op = sparse.vstack((sparse.eye(2 * n_pts), + hxy2hz)) + return op -def v2e(v: numpy.ndarray, - wavenumber: complex, - omega: complex, - dxes: dx_lists_t, - epsilon: vfield_t, - mu: vfield_t = None - ) -> vfield_t: +def exy2e(wavenumber: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + ) -> sparse.spmatrix: """ - Given a vector v containing the vectorized H_x and H_y fields, - returns a vectorized E containing all three E components + 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 v: Vector containing H_x and H_y fields - :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v - :param omega: The angular frequency of the system + :param 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 - :param mu: Vectorized magnetic permeability grid (default 1 everywhere) - :return: Vectorized E field with all vector components. + :return: Sparse matrix representing the operator """ - h2eop = h2e(wavenumber, omega, dxes, epsilon) - return h2eop @ v2h(v, wavenumber, dxes, mu) + 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, diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py index 3dbdfd8..300d3e6 100644 --- a/meanas/fdfd/waveguide_mode.py +++ b/meanas/fdfd/waveguide_mode.py @@ -12,6 +12,7 @@ def solve_waveguide_mode_2d(mode_number: int, dxes: dx_lists_t, epsilon: vfield_t, mu: vfield_t = None, + mode_margin: int = 2, ) -> Dict[str, complex or field_t]: """ Given a 2d region, attempts to solve for the eigenmode with the specified mode number. @@ -21,6 +22,9 @@ 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 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': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} """ @@ -28,23 +32,23 @@ def solve_waveguide_mode_2d(mode_number: int, Solve for the largest-magnitude eigenvalue of the real operator ''' dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] - A_r = waveguide.operator(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu)) + A_r = waveguide.operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu)) - eigvals, eigvecs = signed_eigensolve(A_r, mode_number+3) - v = eigvecs[:, -(mode_number + 1)] + eigvals, eigvecs = signed_eigensolve(A_r, mode_number + mode_margin) + exy = 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(omega, dxes, epsilon, mu) - eigval, v = rayleigh_quotient_iteration(A, v) + A = waveguide.operator_e(omega, dxes, epsilon, mu) + eigval, exy = rayleigh_quotient_iteration(A, exy) # Calculate the wave-vector (force the real part to be positive) wavenumber = numpy.sqrt(eigval) wavenumber *= numpy.sign(numpy.real(wavenumber)) - e, h = waveguide.normalized_fields(v, wavenumber, omega, dxes, epsilon, mu) + e, h = waveguide.normalized_fields_e(exy, wavenumber, omega, dxes, epsilon, mu) shape = [d.size for d in dxes[0]] fields = {