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

@@ -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,
+def normalized_fields_e(e_xy: numpy.ndarray,
-                      wavenumber: complex,
+                        wavenumber: complex,
-                      omega: complex,
+                        omega: complex,
-                      dxes: dx_lists_t,
+                        dxes: dx_lists_t,
-                      epsilon: vfield_t,
+                        epsilon: vfield_t,
-                      mu: vfield_t = None,
+                        mu: vfield_t = None,
-                      dx_prop: float = 0,
+                        dx_prop: float = 0,
-                      ) -> Tuple[vfield_t, vfield_t]:
+                        ) -> 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 e_xy: Vector containing E_x and E_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
@@ -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)
+    e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy
-    h = v2h(v, wavenumber, dxes, mu=mu)
+    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,
+def exy2h(wavenumber: complex,
-        wavenumber: complex,
+          omega: complex,
-        dxes: dx_lists_t,
+          dxes: dx_lists_t,
-        mu: vfield_t = None
+          epsilon: vfield_t,
-        ) -> vfield_t:
+          mu: vfield_t = None
+          ) -> sparse.spmatrix:
     """
-    Given a vector v containing the vectorized H_x and H_y fields,
+    Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields,
-     returns a vectorized H including all three H components.
+     into a vectorized H containing all three H components
 
-    :param v: Vector containing H_x and H_y fields
+    :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
-    :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: 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)
+    n_pts = dxes[1][0].size * dxes[1][1].size
-    return numpy.hstack((v, w)).flatten()
+    op = sparse.vstack((sparse.eye(2 * n_pts),
+                        hxy2hz))
+    return op
 
 
-def v2e(v: numpy.ndarray,
+def exy2e(wavenumber: complex,
-        wavenumber: complex,
+          dxes: dx_lists_t,
-        omega: complex,
+          epsilon: vfield_t,
-        dxes: dx_lists_t,
+          ) -> sparse.spmatrix:
-        epsilon: vfield_t,
-        mu: vfield_t = None
-        ) -> vfield_t:
     """
-    Given a vector v containing the vectorized H_x and H_y fields,
+    Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields,
-     returns a vectorized E containing all three E components
+     into a vectorized E containing all three E components
 
-    :param v: Vector containing H_x and H_y fields
+    :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
-    :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: Sparse matrix representing the operator
-    :return: Vectorized E field with all vector components.
     """
-    h2eop = h2e(wavenumber, omega, dxes, epsilon)
+    Dbx, Dby = operators.deriv_back(dxes[1])
-    return h2eop @ v2h(v, wavenumber, dxes, mu)
+    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,

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)
+    eigvals, eigvecs = signed_eigensolve(A_r, mode_number + mode_margin)
-    v = eigvecs[:, -(mode_number + 1)]
+    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)
+    A = waveguide.operator_e(omega, dxes, epsilon, mu)
-    eigval, v = rayleigh_quotient_iteration(A, v)
+    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 = {