moving to backend

Author: momchil-flex

tests/test_plugins/test_mode_solver.py

@@ -1237,87 +1237,6 @@ 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=(-1, 0, 0),
-        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)
-
-    # Modes with and without an angle are close
-    assert np.allclose(modes_zero.n_complex.values, modes_angle.n_complex.values, atol=2e-2)
-    # Modes with angle (tensorial mode solver) have a small k_eff
-    assert np.all(np.abs(modes_angle.k_eff) < 1e-10)
-
-
 def test_gauge_robustness():
     array = np.zeros((5, 5), dtype=float)
     ij = np.arange(5)
@@ -1367,7 +1286,3 @@ def test_translated_dot():
 
     assert np.allclose(data2.outer_dot(data_translated), data2.outer_dot(data2), atol=atol)
     assert np.allclose(data_translated.outer_dot(data2), data2.outer_dot(data2), atol=atol)
-
-
-if __name__ == "__main__":
-    test_mode_solver_angle_symmetry()

tidy3d/components/mode/derivatives.py

@@ -240,119 +240,3 @@ 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."""
-
-    # return make_cxf(dls, shape, pmc).T
-
-    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."""
-
-    # return make_cyf(dls, shape, pmc).T
-
-    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)