diff --git a/tests/test_components/test_mode.py b/tests/test_components/test_mode.py index df08f470a3..30f81a82de 100644 --- a/tests/test_components/test_mode.py +++ b/tests/test_components/test_mode.py @@ -244,3 +244,10 @@ def get_mode_sim_data(): def test_mode_sim_data(): sim_data = get_mode_sim_data() _ = sim_data.plot_field("Ey", ax=AX, mode_index=0, f=FS[0]) + + +def test_mode_boundary(): + _ = td.ModeSpec(boundary="PEC") + _ = td.ModeSpec(boundary="PMC") + with pytest.raises(pydantic.ValidationError): + _ = td.ModeSpec(boundary="PBC") diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index 911916d4b3..44d44709fc 100644 --- a/tests/test_plugins/test_mode_solver.py +++ b/tests/test_plugins/test_mode_solver.py @@ -167,11 +167,12 @@ def mock_download(resource_id, remote_filename, to_file, *args, **kwargs): ) -def test_compute_modes(): +@pytest.mark.parametrize("boundary", ("PEC", "PMC")) +def test_compute_modes(boundary): """Test direct call to `compute_modes`.""" eps_cross = np.random.rand(10, 10) coords = np.arange(11) - mode_spec = td.ModeSpec(num_modes=3, target_neff=2.0) + mode_spec = td.ModeSpec(num_modes=3, target_neff=2.0, boundary=boundary) _ = compute_modes( eps_cross=[eps_cross] * 9, coords=[coords, coords], @@ -736,14 +737,16 @@ def test_mode_bend_radius(): assert np.allclose(data1.n_complex, data2.n_complex * 5.5 / 5.0) -def test_mode_solver_2D(): +@pytest.mark.parametrize("boundary", ("PEC", "PMC")) +def test_mode_solver_2D(boundary): """Run mode solver in 2D simulations.""" mode_spec = td.ModeSpec( num_modes=3, - filter_pol="te", + filter_pol="te" if boundary == "PEC" else "tm", precision="double", num_pml=(0, 10), track_freq="central", + boundary=boundary, ) simulation = td.Simulation( size=(0, SIM_SIZE[1], SIM_SIZE[2]), @@ -764,9 +767,10 @@ def test_mode_solver_2D(): mode_spec = td.ModeSpec( num_modes=3, - filter_pol="te", + filter_pol="te" if boundary == "PEC" else "tm", precision="double", num_pml=(10, 0), + boundary=boundary, ) simulation = td.Simulation( size=(SIM_SIZE[0], SIM_SIZE[1], 0), diff --git a/tidy3d/components/mode/derivatives.py b/tidy3d/components/mode/derivatives.py index 4c13ca7835..6c3f24451c 100644 --- a/tidy3d/components/mode/derivatives.py +++ b/tidy3d/components/mode/derivatives.py @@ -6,13 +6,15 @@ from ...constants import EPSILON_0, ETA_0 -def make_dxf(dls, shape, pmc): +def make_dxf(dls, shape, pmc_neg, pmc_pos): """Forward derivative in x.""" Nx, Ny = shape if Nx == 1: return sp.csr_matrix((Ny, Ny)) dxf = sp.csr_matrix(sp.diags([-1, 1], [0, 1], shape=(Nx, Nx))) - if not pmc: + if pmc_pos: + dxf[-1, -1] = 0.0 + if not pmc_neg: dxf[0, 0] = 0.0 dxf = sp.diags(1 / dls).dot(dxf) dxf = sp.kron(dxf, sp.eye(Ny)) @@ -34,13 +36,15 @@ def make_dxb(dls, shape, pmc): return dxb -def make_dyf(dls, shape, pmc): +def make_dyf(dls, shape, pmc_neg, pmc_pos): """Forward derivative in y.""" Nx, Ny = shape if Ny == 1: return sp.csr_matrix((Nx, Nx)) dyf = sp.csr_matrix(sp.diags([-1, 1], [0, 1], shape=(Ny, Ny))) - if not pmc: + if pmc_pos: + dyf[-1, -1] = 0.0 + if not pmc_neg: dyf[0, 0] = 0.0 dyf = sp.diags(1 / dls).dot(dyf) dyf = sp.kron(sp.eye(Nx), dyf) @@ -62,15 +66,15 @@ def make_dyb(dls, shape, pmc): return dyb -def create_d_matrices(shape, dls, dmin_pmc=(False, False)): +def create_d_matrices(shape, dls, dmin_pmc=(False, False), dmax_pmc=False): """Make the derivative matrices without PML. If dmin_pmc is True, the 'backward' derivative in that dimension will be set to implement PMC boundary, otherwise it will be set to PEC.""" dlf, dlb = dls - dxf = make_dxf(dlf[0], shape, dmin_pmc[0]) + dxf = make_dxf(dlf[0], shape, dmin_pmc[0], dmax_pmc) dxb = make_dxb(dlb[0], shape, dmin_pmc[0]) - dyf = make_dyf(dlf[1], shape, dmin_pmc[1]) + dyf = make_dyf(dlf[1], shape, dmin_pmc[1], dmax_pmc) dyb = make_dyb(dlb[1], shape, dmin_pmc[1]) return (dxf, dxb, dyf, dyb) diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index 2cd12ef4de..37c328c244 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -203,7 +203,8 @@ def compute_modes( always impose PEC boundary at the xmax and ymax interfaces, and on the xmin and ymin interfaces unless PMC symmetry is present. If so, the PMC boundary is imposed through the backward derivative matrices.""" - dmin_pmc = [sym == 1 for sym in symmetry] + dmax_pmc = mode_spec.boundary == "PMC" + dmin_pmc = [sym == 1 or (sym == 0 and dmax_pmc) for sym in symmetry] # Primal grid steps for E-field derivatives dl_f = [new_cs[1:] - new_cs[:-1] for new_cs in new_coords] @@ -213,7 +214,7 @@ def compute_modes( dls = (dl_f, dl_b) # Derivative matrices with PEC boundaries by default and optional PMC at the near end - der_mats_tmp = d_mats(Nxy, dls, dmin_pmc) + der_mats_tmp = d_mats(Nxy, dls, dmin_pmc, dmax_pmc) # PML matrices; do not impose PML on the bottom when symmetry present dmin_pml = np.array(symmetry) == 0 diff --git a/tidy3d/components/mode_spec.py b/tidy3d/components/mode_spec.py index 4f63465bce..15a2ac2fef 100644 --- a/tidy3d/components/mode_spec.py +++ b/tidy3d/components/mode_spec.py @@ -160,6 +160,12 @@ class ModeSpec(Tidy3dBaseModel): f"default of {GROUP_INDEX_STEP} is used.", ) + boundary: Literal["PEC", "PMC"] = pd.Field( + "PEC", + title="Boundary Conditions", + description="Boundary conditions to impose at the edges of the mode plane.", + ) + @pd.validator("bend_axis", always=True) @skip_if_fields_missing(["bend_radius"]) def bend_axis_given(cls, val, values):