Introducing colocation matrices in tensorial mode solver

Author: momchil-flex

docs/notebooks

@@ -1237,6 +1237,84 @@ def test_high_order_mode_normalization():
     assert (values[: values.size // 3] > 0).all()
 
 
+# Test that a mode solver with an angle and PEC symmetry does not produce unphysical modes
+def test_mode_solver_angle_symmetry():
+    # Angle to test
+    theta = np.pi / 180
+    theta = np.pi / 6
+
+    # Other parameters
+    plane = td.Box(center=(0, 0, 0), size=(5, 0, 8))
+    sim_size = (10, 10, 10)
+    # sim_size = (4, 3, 3)
+    freq0 = td.C_0 / 1.0
+
+    # Make the default waveguide lossless to also test numerical loss artifact in the angled case
+    wg = WAVEGUIDE.updated_copy(medium=td.Medium(permittivity=4))
+
+    # ModeSource at an angle for visualization aid
+    src = td.ModeSource(
+        source_time=td.GaussianPulse(freq0=2e14, fwidth=1e13),
+        center=plane.center,
+        size=plane.size,
+        mode_spec=td.ModeSpec(num_modes=1, angle_theta=theta),
+        mode_index=0,
+        direction="+",
+    )
+
+    # Mode solver without angle
+    simulation = td.Simulation(
+        size=sim_size,
+        grid_spec=td.GridSpec.auto(wavelength=1.0, min_steps_per_wvl=10),
+        # grid_spec=td.GridSpec.uniform(dl=0.04),
+        structures=[wg],
+        run_time=1e-12,
+        symmetry=(0, 0, -1),
+        sources=[src],
+        # medium=td.Medium(permittivity=4.5),
+    )
+
+    # simulation = td.Simulation.from_file("/home/momchil/flexcompute/test_fdtd/lucas/FDTD_Setup_P1@0_v1.json")
+    # print(simulation.grid.num_cells)
+    # freq0 = simulation.sources[0].source_time.freq0
+    # plane = td.Box(center=simulation.monitors[1].center, size=simulation.monitors[1].size)
+
+    mode_spec = td.ModeSpec(num_modes=5)
+    ms = ModeSolver(
+        simulation=simulation,
+        plane=plane,
+        mode_spec=mode_spec,
+        freqs=[freq0],
+    )
+    modes_zero = ms.solve()
+
+    # Compute modes with angle
+    wg = wg.updated_copy(geometry=wg.geometry.rotated(angle=-theta, axis=[0, 0, 1]))
+    simulation = simulation.updated_copy(structures=[wg])
+    mode_spec = mode_spec.updated_copy(angle_theta=theta)
+    ms_angle = ms.updated_copy(mode_spec=mode_spec, simulation=simulation)
+    modes_angle = ms_angle.solve()
+
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(1, 3)
+    simulation.plot(x=0, ax=ax[0])
+    simulation.plot(y=0, ax=ax[1])
+    simulation.plot(z=0, ax=ax[2])
+    plt.show()
+
+    # Plot the mode_index = 0 mode for both
+    _, ax = plt.subplots(1, 2)
+    ms.plot_field("E", mode_index=0, ax=ax[0], robust=False)
+    ms_angle.plot_field("E", mode_index=0, ax=ax[1], robust=False)
+    plt.show()
+
+    print(modes_zero.n_complex.values)
+    print(modes_angle.n_complex.values)
+
+    assert np.allclose(modes_zero.n_complex.values, modes_angle.n_complex.values, atol=3e-2)
+
+
 def test_gauge_robustness():
     array = np.zeros((5, 5), dtype=float)
     ij = np.arange(5)

tests/test_plugins/test_mode_solver.py

@@ -240,3 +240,113 @@ def s_value(
     kappa = kappa_min + (kappa_max - kappa_min) * step**order
     sigma = sigma_max * avg_speed / (ETA_0 * dl) * step**order
     return kappa + 1j * sigma / (omega * EPSILON_0)
+
+
+def make_cxf(dls, shape, pmc):
+    """Forward colocation in x."""
+    import scipy.sparse as sp
+
+    Nx, Ny = shape
+    if Nx == 1:
+        return sp.csr_matrix((Ny, Ny))
+    # Create matrix with [1,1,0,...,0], [0,1,1,...,0], etc pattern
+    cxf = 0.5 * sp.csr_matrix(sp.diags([1, 1], [0, 1], shape=(Nx, Nx)))
+
+    # Apply boundary conditions before Kronecker product
+    if not pmc:
+        cxf[0, 0] = 0
+
+    return sp.kron(cxf, sp.eye(Ny))
+
+
+def make_cxb(dls, shape, pmc):
+    """Backward colocation in x."""
+    import scipy.sparse as sp
+
+    Nx, Ny = shape
+    if Nx == 1:
+        return sp.csr_matrix((Ny, Ny))
+
+    # Calculate weights for each position
+    p1 = dls[:-1]  # First grid spacing
+    p2 = dls[1:]  # Second grid spacing
+
+    weights_curr = np.ones(Nx)
+    weights_curr[1:] = p1 / (p1 + p2)
+
+    weights_prev = np.zeros(Nx)
+    weights_prev[:-1] = p2 / (p1 + p2)
+
+    # Create the matrix using diagonals
+    cxb = sp.diags([weights_curr, weights_prev], [0, -1], shape=(Nx, Nx))
+
+    # Apply boundary conditions before Kronecker product
+    # The matrix is already set up for PEC (cxb[0, 0] = 1)
+    if pmc:
+        cxb[0, 0] = 0
+
+    return sp.kron(cxb, sp.eye(Ny))
+
+
+def make_cyf(dls, shape, pmc):
+    """Forward colocation in y."""
+    import scipy.sparse as sp
+
+    Nx, Ny = shape
+    if Ny == 1:
+        return sp.csr_matrix((Nx, Nx))
+    # Create matrix with [1,1,0,...,0], [0,1,1,...,0], etc pattern
+    cyf = 0.5 * sp.csr_matrix(sp.diags([1, 1], [0, 1], shape=(Ny, Ny)))
+
+    # Apply boundary conditions before Kronecker product
+    if not pmc:
+        cyf[0, 0] = 0
+
+    return sp.kron(sp.eye(Nx), cyf)
+
+
+def make_cyb(dls, shape, pmc):
+    """Backward colocation in y."""
+    import scipy.sparse as sp
+
+    Nx, Ny = shape
+    if Ny == 1:
+        return sp.csr_matrix((Nx, Nx))
+
+    # Calculate weights for each position
+    p1 = dls[:-1]  # First grid spacing
+    p2 = dls[1:]  # Second grid spacing
+
+    weights_curr = np.ones(Ny)
+    weights_curr[1:] = p1 / (p1 + p2)
+
+    weights_prev = np.zeros(Ny)
+    weights_prev[:-1] = p2 / (p1 + p2)
+
+    # Create the matrix using diagonals
+    cyb = sp.diags([weights_curr, weights_prev], [0, -1], shape=(Ny, Ny))
+
+    # Apply boundary conditions before Kronecker product
+    # The matrix is already set up for PEC (cyb[0, 0] = 1)
+    if pmc:
+        cyb[0, 0] = 0
+
+    return sp.kron(sp.eye(Nx), cyb)
+
+
+def create_c_matrices(shape, dls, dmin_pmc=(False, False)):
+    """Make the colocation matrices. If dmin_pmc is True, the matrices will be modified
+    to implement PMC boundary conditions, otherwise they will implement PEC."""
+
+    dlf, _ = dls
+    cxf = make_cxf(dlf[0], shape, dmin_pmc[0])
+    cxb = make_cxb(dlf[0], shape, dmin_pmc[0])  # Note: uses primal grid spacings
+    cyf = make_cyf(dlf[1], shape, dmin_pmc[1])
+    cyb = make_cyb(dlf[1], shape, dmin_pmc[1])  # Note: uses primal grid spacings
+
+    # cxf = sp.eye(shape[0] * shape[1])
+    # cxb = sp.eye(shape[0] * shape[1])
+    # cyf = sp.eye(shape[0] * shape[1])
+    # cyb = sp.eye(shape[0] * shape[1])
+
+    return (cxf, cxb, cyf, cyb)