From f04c0daf287722eb31b00dd3c9c4d09d8103fa04 Mon Sep 17 00:00:00 2001 From: jan Date: Thu, 5 Sep 2019 22:38:29 +0200 Subject: [PATCH] solve_waveguide_mode_2d -> vsolve_* - return (e_xy. wavenumber) - vectorized inputs since we returned a vectorized output - exy -> e_xy --- meanas/fdfd/waveguide_mode.py | 66 +++++++++++++++++------------------ 1 file changed, 32 insertions(+), 34 deletions(-) diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py index cfa888a..6bdbc73 100644 --- a/meanas/fdfd/waveguide_mode.py +++ b/meanas/fdfd/waveguide_mode.py @@ -1,4 +1,4 @@ -from typing import Dict, List +from typing import Dict, List, Tuple import numpy import scipy.sparse as sparse @@ -7,13 +7,13 @@ from . import operators, waveguide, functional from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration -def solve_waveguide_mode_2d(mode_number: int, - omega: complex, - dxes: dx_lists_t, - epsilon: field_t, - mu: field_t = None, - mode_margin: int = 2, - ) -> Dict[str, complex or field_t]: +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. @@ -25,39 +25,31 @@ def solve_waveguide_mode_2d(mode_number: int, :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': numpy.ndarray, 'H': numpy.ndarray, 'wavenumber': complex} + :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, vec(numpy.real(epsilon)), vec(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 + mode_margin) - exy = eigvecs[:, -(mode_number + 1)] + 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, vec(epsilon), vec(mu)) - eigval, exy = rayleigh_quotient_iteration(A, exy) + 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)) - e, h = waveguide.normalized_fields_e(exy, wavenumber, omega, dxes, vec(epsilon), vec(mu)) + return e_xy, wavenumber - shape = [d.size for d in dxes[0]] - fields = { - 'wavenumber': wavenumber, - 'E': unvec(e, shape), - 'H': unvec(h, shape), - } - - return fields def solve_waveguide_mode(mode_number: int, @@ -102,36 +94,42 @@ def solve_waveguide_mode(mode_number: int, # 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': [epsilon[i][slices].transpose(order) for i in order], - 'mu': [mu[i][slices].transpose(order) for i in order], + 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]), + 'mu': vec([mu[i][slices].transpose(order) for i in order]), } - fields_2d = solve_waveguide_mode_2d(mode_number, omega=omega, **args_2d) + 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. - print(fields_2d['wavenumber'] / (2/dx_prop * numpy.arcsin(fields_2d['wavenumber'] * dx_prop/2))) - print(fields_2d['wavenumber'].real / (2/dx_prop * numpy.arcsin(fields_2d['wavenumber'].real * dx_prop/2))) - fields_2d['wavenumber'] = 2/dx_prop * numpy.arcsin(fields_2d['wavenumber'] * dx_prop/2) + 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 - fields_2d['H'] *= polarity + h *= polarity # Apply phase shift to H-field - fields_2d['H'][:2] *= numpy.exp(-1j * polarity * 0.5 * fields_2d['wavenumber'] * dx_prop) - fields_2d['E'][2] *= numpy.exp(-1j * polarity * 0.5 * fields_2d['wavenumber'] * dx_prop) + 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)] = fields_2d['E'][o][:, :, None].transpose(reverse_order) - H[(a, *slices)] = fields_2d['H'][o][:, :, None].transpose(reverse_order) + E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order) + H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order) results = { - 'wavenumber': fields_2d['wavenumber'], + 'wavenumber': wavenumber, + 'wavenumber_2d': wavenumber_2d, 'H': H, 'E': E, }