From bb4322ded9d6c7f66d54b32ff8d0426608c53756 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Thu, 25 Jun 2026 12:00:48 -0700 Subject: [PATCH] [waveguide_cyl] adjust operators to match waveguide_2d in the large-radius limit --- meanas/fdfd/waveguide_cyl.py | 100 +++++++++++++++-------------------- 1 file changed, 42 insertions(+), 58 deletions(-) diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index e1f8fd0..96a2ff6 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -43,39 +43,9 @@ T_b &= \operatorname{diag}(r_b / r_{\min}). $$ With the same forward/backward derivative notation used in `waveguide_2d`, the -coordinate-transformed discrete curl equations used here are - -$$ -\begin{aligned} --\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\ --\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r - - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\ --\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\ -\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\ -\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y - - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\ -\imath \omega E_z &= T_a \epsilon_{zz}^{-1} - \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right). -\end{aligned} -$$ - -The first three equations are the cylindrical analogue of the straight-guide -relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted -versions of the straight-guide identities for `\imath \beta H_y` and -`\imath \beta H_r`, and the last equation plays the same role as the -longitudinal `E_z` reconstruction in `waveguide_2d`. - -Following the same elimination steps as in `waveguide_2d`, apply -`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the -equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`, -and then eliminate `H_z` with - -$$ -H_z = \frac{1}{-\imath \omega \mu_{zz}} -\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right). -$$ - -This yields the transverse electric eigenproblem implemented by +implementation treats the transverse electric eigenproblem as the canonical +cylindrical discretization. It reduces to `waveguide_2d.operator_e(...)` in the +large-radius limit `T_a, T_b \to I`. The eigenproblem implemented by `cylindrical_operator(...)`: $$ @@ -111,6 +81,33 @@ T_a \epsilon_{zz}^{-1} \begin{bmatrix} E_r \\ E_y \end{bmatrix}. $$ +The full fields reconstructed by `exy2e(...)` and `e2h(...)` use the matching +large-radius-compatible identities + +$$ +E_z = +\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1} +\begin{bmatrix} +\hat{\partial}_r T_b \epsilon_{rr} & +\hat{\partial}_y T_a \epsilon_{yy} +\end{bmatrix} +\begin{bmatrix} E_r \\ E_y \end{bmatrix}, +$$ + +and + +$$ +\begin{bmatrix} H_r \\ H_y \\ H_z \end{bmatrix} += +\frac{1}{-\imath \omega}\mu^{-1} +\begin{bmatrix} +0 & \imath\beta T_a^{-1} & \tilde{\partial}_y \\ +-\imath\beta T_b^{-1} & 0 & -T_b^{-1}\tilde{\partial}_r T_a \\ +-\tilde{\partial}_y & \tilde{\partial}_r & 0 +\end{bmatrix} +\begin{bmatrix} E_r \\ E_y \\ E_z \end{bmatrix}. +$$ + Since `\beta = m / r_{\min}`, the solver implemented in this file returns the angular wavenumber `m`, while the operator itself is most naturally written in terms of the linear quantity `\beta`. The helpers below reconstruct the full @@ -367,7 +364,6 @@ def exy2h( def exy2e( angular_wavenumber: complex, - omega: float, dxes: dx_lists2_t, rmin: float, epsilon: vfdslice, @@ -383,7 +379,6 @@ def exy2e( angular_wavenumber: Wavenumber assuming fields have theta-dependence of `exp(-i * angular_wavenumber * theta)`. It should satisfy `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy` - omega: The angular frequency of the system dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D) rmin: Radius at the left edge of the simulation domain (at minimum 'x') epsilon: Vectorized dielectric constant grid @@ -391,30 +386,22 @@ def exy2e( Returns: Sparse matrix representing the operator. """ - Dfx, Dfy = deriv_forward(dxes[0]) Dbx, Dby = deriv_back(dxes[1]) wavenumber = angular_wavenumber / rmin Ta, Tb = dxes2T(dxes=dxes, rmin=rmin) - Tai = sparse.diags_array(1 / Ta.diagonal()) - #Tbi = sparse.diags_array(1 / Tb.diagonal()) epsilon_parts = numpy.split(epsilon, 3) epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2]) epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2]) n_pts = dxes[0][0].size * dxes[0][1].size - zeros = sparse.coo_array((n_pts, n_pts)) - - mu_z = numpy.ones(n_pts) - mu_z_inv = sparse.diags_array(1 / mu_z) - exy2hz = 1 / (-1j * omega) * mu_z_inv @ sparse.hstack((Dfy, -Dfx)) - hxy2ez = 1 / (1j * omega) * epsilon_z_inv @ sparse.hstack((Dby, -Dbx)) - - exy2hy = Tb / (1j * wavenumber) @ (-1j * omega * sparse.hstack((epsilon_x, zeros)) - Dby @ exy2hz) - exy2hx = Tb / (1j * wavenumber) @ ( 1j * omega * sparse.hstack((zeros, epsilon_y)) - Tai @ Dbx @ Tb @ exy2hz) - - exy2ez = hxy2ez @ sparse.vstack((exy2hx, exy2hy)) + exy2ez = ( + Ta @ epsilon_z_inv + @ sparse.hstack((Dbx @ Tb @ epsilon_x, + Dby @ Ta @ epsilon_y)) + / (1j * wavenumber) + ) op = sparse.vstack((sparse.eye_array(2 * n_pts), exy2ez)) @@ -460,9 +447,9 @@ def e2h( Tbi = sparse.diags_array(1 / Tb.diagonal()) jB = 1j * angular_wavenumber / rmin - op = sparse.block_array([[ None, -jB * Tai, -Dfy], - [jB * Tbi, None, Tbi @ Dfx @ Ta], - [ Dfy, -Dfx, None]]) / (-1j * omega) + op = sparse.block_array([[ None, jB * Tai, Dfy], + [-jB * Tbi, None, -Tbi @ Dfx @ Ta], + [ -Dfy, Dfx, None]]) / (-1j * omega) if mu is not None: op = sparse.diags_array(1 / mu) @ op return op @@ -565,19 +552,16 @@ def _normalized_fields( The normalization procedure is: - 1. Flip the magnetic field sign so the reconstructed `(e, h)` pair follows - the same forward-power convention as `waveguide_2d`. - 2. Compute the time-averaged forward power with + 1. Compute the time-averaged forward power with `waveguide_2d.inner_product(..., conj_h=True)`. - 3. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power. - 4. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic + 2. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power. + 3. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic so symmetric modes choose a stable representative. `prop_phase` has the same meaning as in `waveguide_2d`: it compensates for the half-cell longitudinal staggering between the E and H fields when the propagation direction is itself discretized. """ - h *= -1 shape = [s.size for s in dxes[0]] # Find time-averaged Sz and normalize to it