diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1cb8c9a13..e01809f16 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,10 +12,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Integration of the `documentation` alongside the main codebase repository.
 - Integration of the `tidy3d-notebooks` repository.
 - `tidy3d develop` CLI and development guide on the main documentation.
+- Added a convenience method `Simulation.subsection()` to a create a new simulation based on a subregion of another one.
 
 ### Changed
 - `poetry` based installation. Removal of `setup.py` and `requirements.txt`.
 - Upgrade to sphinx 6 for the documentation build, and change of theme.
+- Remote mode solver web api automatically reduces the associated `Simulation` object to the mode solver plane before uploading it to server.
 
 ### Fixed
 
diff --git a/tests/test_components/test_custom.py b/tests/test_components/test_custom.py
index e95731da0..b95802e48 100644
--- a/tests/test_components/test_custom.py
+++ b/tests/test_components/test_custom.py
@@ -351,7 +351,7 @@ def test_n_cfl():
     assert med.n_cfl >= 2
 
 
-def verify_custom_medium_methods(mat):
+def verify_custom_medium_methods(mat, reduced_fields=[]):
     """Verify that the methods in custom medium is producing expected results."""
     freq = 1.0
     assert isinstance(mat, AbstractCustomMedium)
@@ -362,6 +362,41 @@ def verify_custom_medium_methods(mat):
     for i in range(3):
         assert np.allclose(eps_grid[i].shape, [len(f) for f in coord_interp.to_list])
 
+    # check reducing data
+    subsection = td.Box(size=(0.3, 0.4, 0.35), center=(0.4, 0.4, 0.4))
+
+    mat_reduced = mat.sel_inside(subsection.bounds)
+
+    for field in reduced_fields:
+        original = getattr(mat, field)
+        reduced = getattr(mat_reduced, field)
+
+        if original is None:
+            assert reduced is None
+            continue
+
+        # data fields in medium classes could be SpatialArrays or 2d tuples of spatial arrays
+        # lets convert everything into 2d tuples of spatial arrays for uniform handling
+        if isinstance(original, td.SpatialDataArray):
+            original = [
+                [
+                    original,
+                ],
+            ]
+            reduced = [
+                [
+                    reduced,
+                ],
+            ]
+
+        for or_set, re_set in zip(original, reduced):
+            assert len(or_set) == len(re_set)
+
+            for ind in range(len(or_set)):
+                diff = (or_set[ind] - re_set[ind]).abs
+                assert diff.does_cover(subsection.bounds)
+                assert np.allclose(diff, 0)
+
     # construct sim
     struct = td.Structure(
         geometry=td.Box(size=(0.5, 0.5, 0.5)),
@@ -375,6 +410,8 @@ def verify_custom_medium_methods(mat):
         structures=(struct,),
     )
     _ = sim.grid
+    sim_reduced = sim.subsection(subsection, remove_outside_custom_mediums=False)
+    sim_reduced = sim.subsection(subsection, remove_outside_custom_mediums=True)
 
     # bkg
     sim = td.Simulation(
@@ -384,6 +421,8 @@ def verify_custom_medium_methods(mat):
         medium=mat,
     )
     _ = sim.grid
+    sim_reduced = sim.subsection(subsection, remove_outside_custom_mediums=False)
+    sim_reduced = sim.subsection(subsection, remove_outside_custom_mediums=True)
 
 
 def test_anisotropic_custom_medium():
@@ -438,7 +477,7 @@ def test_custom_isotropic_medium():
     with pytest.raises(pydantic.ValidationError):
         mat = CustomMedium(permittivity=permittivity, conductivity=sigmatmp)
     mat = CustomMedium(permittivity=permittivity, conductivity=sigmatmp, allow_gain=True)
-    verify_custom_medium_methods(mat)
+    verify_custom_medium_methods(mat, ["permittivity", "conductivity"])
 
     # inconsistent coords
     with pytest.raises(pydantic.ValidationError):
@@ -448,15 +487,15 @@ def test_custom_isotropic_medium():
         mat = CustomMedium(permittivity=permittivity, conductivity=sigmatmp)
 
     mat = CustomMedium(permittivity=permittivity, conductivity=conductivity)
-    verify_custom_medium_methods(mat)
+    verify_custom_medium_methods(mat, ["permittivity", "conductivity"])
 
     mat = CustomMedium(permittivity=permittivity)
-    verify_custom_medium_methods(mat)
+    verify_custom_medium_methods(mat, ["permittivity", "conductivity"])
 
 
-def verify_custom_dispersive_medium_methods(mat):
+def verify_custom_dispersive_medium_methods(mat, reduced_fields=[]):
     """Verify that the methods in custom dispersive medium is producing expected results."""
-    verify_custom_medium_methods(mat)
+    verify_custom_medium_methods(mat, reduced_fields)
     freq = 1.0
     for i in range(3):
         assert mat.eps_dataarray_freq(freq)[i].shape == (Nx, Ny, Nz)
@@ -515,7 +554,7 @@ def test_custom_pole_residue():
 
     eps_inf = td.SpatialDataArray(np.random.random((Nx, Ny, Nz)) + 1, coords=dict(x=X, y=Y, z=Z))
     mat = CustomPoleResidue(eps_inf=eps_inf, poles=((a, c),))
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "poles"])
     assert mat.n_cfl > 1
 
     # to custom non-dispersive medium
@@ -529,12 +568,12 @@ def test_custom_pole_residue():
         mat_medium = mat.to_medium()
     mat = CustomPoleResidue(eps_inf=eps_inf, poles=((a, c - 0.1),), allow_gain=True)
     mat_medium = mat.to_medium()
-    verify_custom_medium_methods(mat_medium)
+    verify_custom_medium_methods(mat_medium, ["permittivity", "conductivity"])
     assert mat_medium.n_cfl > 1
 
     # custom medium to pole residue
     mat = CustomPoleResidue.from_medium(mat_medium)
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "poles"])
     assert mat.n_cfl > 1
 
 
@@ -578,14 +617,14 @@ def test_custom_sellmeier():
         mat = CustomSellmeier(coeffs=((b1, c2), (btmp, c2)))
 
     mat = CustomSellmeier(coeffs=((b1, c1), (b2, c2)))
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["coeffs"])
     assert mat.n_cfl == 1
 
     # from dispersion
     n = td.SpatialDataArray(2 + np.random.random((Nx, Ny, Nz)), coords=dict(x=X, y=Y, z=Z))
     dn_dwvl = td.SpatialDataArray(-np.random.random((Nx, Ny, Nz)), coords=dict(x=X, y=Y, z=Z))
     mat = CustomSellmeier.from_dispersion(n=n, dn_dwvl=dn_dwvl, freq=2, interp_method="linear")
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["coeffs"])
     assert mat.n_cfl == 1
 
 
@@ -638,13 +677,13 @@ def test_custom_lorentz():
     mat = CustomLorentz(
         eps_inf=eps_inf, coeffs=((de1, f1, delta1), (detmp, f2, delta2)), allow_gain=True
     )
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "coeffs"])
     assert mat.n_cfl > 1
 
     mat = CustomLorentz(
         eps_inf=eps_inf, coeffs=((de1, f1, delta1), (de2, f2, delta2)), subpixel=True
     )
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "coeffs"])
     assert mat.n_cfl > 1
     assert mat.pole_residue.subpixel
 
@@ -681,7 +720,7 @@ def test_custom_drude():
         mat = CustomDrude(eps_inf=eps_inf, coeffs=((f1, delta1), (ftmp, delta2)))
 
     mat = CustomDrude(eps_inf=eps_inf, coeffs=((f1, delta1), (f2, delta2)))
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "coeffs"])
     assert mat.n_cfl > 1
 
 
@@ -728,11 +767,11 @@ def test_custom_debye():
         )
         mat = CustomDebye(eps_inf=eps_inf, coeffs=((eps1, tau1), (epstmp, tau2)))
     mat = CustomDebye(eps_inf=eps_inf, coeffs=((eps1, tau1), (epstmp, tau2)), allow_gain=True)
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "coeffs"])
     assert mat.n_cfl > 1
 
     mat = CustomDebye(eps_inf=eps_inf, coeffs=((eps1, tau1), (eps2, tau2)))
-    verify_custom_dispersive_medium_methods(mat)
+    verify_custom_dispersive_medium_methods(mat, ["eps_inf", "coeffs"])
     assert mat.n_cfl > 1
 
 
diff --git a/tests/test_components/test_simulation.py b/tests/test_components/test_simulation.py
index 3b738637a..549534fda 100644
--- a/tests/test_components/test_simulation.py
+++ b/tests/test_components/test_simulation.py
@@ -2248,3 +2248,76 @@ def test_to_gds(tmp_path):
     assert np.allclose(areas[(2, 1)], 0.5)
     assert np.allclose(areas[(1, 0)], 0.25 * np.pi * 1.4**2, atol=1e-2)
     assert np.allclose(areas[(0, 0)], 0.25 * np.pi * 1.4**2, atol=1e-2)
+
+
+def test_sim_subsection():
+    region = td.Box(size=(0.3, 0.5, 0.7), center=(0.1, 0.05, 0.02))
+
+    sim_red = SIM_FULL.subsection(region=region)
+    assert sim_red.structures != SIM_FULL.structures
+    sim_red = SIM_FULL.subsection(
+        region=region,
+        symmetry=(1, 0, -1),
+        monitors=[
+            mnt
+            for mnt in SIM_FULL.monitors
+            if not isinstance(mnt, (td.ModeMonitor, td.ModeSolverMonitor))
+        ],
+    )
+    assert sim_red.symmetry == (1, 0, -1)
+    sim_red = SIM_FULL.subsection(
+        region=region, boundary_spec=td.BoundarySpec.all_sides(td.Periodic())
+    )
+    sim_red = SIM_FULL.subsection(region=region, sources=[], grid_spec=td.GridSpec.uniform(dl=20))
+    assert len(sim_red.sources) == 0
+    sim_red = SIM_FULL.subsection(region=region, monitors=[])
+    assert len(sim_red.monitors) == 0
+    sim_red = SIM_FULL.subsection(region=region, remove_outside_structures=False)
+    assert sim_red.structures == SIM_FULL.structures
+    sim_red = SIM_FULL.subsection(region=region, remove_outside_custom_mediums=True)
+
+    fine_custom_medium = td.CustomMedium(
+        permittivity=td.SpatialDataArray(
+            1 + np.random.random((11, 12, 13)),
+            coords=dict(
+                x=np.linspace(-0.51, 0.52, 11),
+                y=np.linspace(-1.02, 1.04, 12),
+                z=np.linspace(-1.51, 1.51, 13),
+            ),
+        )
+    )
+
+    sim = SIM_FULL.updated_copy(
+        structures=[
+            td.Structure(
+                geometry=td.Box(size=(1, 2, 3)),
+                medium=fine_custom_medium,
+            )
+        ],
+        medium=fine_custom_medium,
+    )
+    sim_red = sim.subsection(region=region, remove_outside_custom_mediums=True)
+
+    # check automatic symmetry expansion
+    sim_sym = SIM_FULL.updated_copy(
+        symmetry=(-1, 0, 1),
+        sources=[src for src in SIM_FULL.sources if not isinstance(src, td.TFSF)],
+    )
+    sim_red = sim_sym.subsection(region=region)
+    assert np.allclose(sim_red.center, (0, 0.05, 0.0))
+
+    # check grid is preserved when requested
+    sim_red = SIM_FULL.subsection(
+        region=region, grid_spec="identical", boundary_spec=td.BoundarySpec.all_sides(td.Periodic())
+    )
+    grids_1d = SIM_FULL.grid.boundaries
+    grids_1d_red = sim_red.grid.boundaries
+    tol = 1e-8
+    for full_grid, red_grid in zip(
+        [grids_1d.x, grids_1d.y, grids_1d.z], [grids_1d_red.x, grids_1d_red.y, grids_1d_red.z]
+    ):
+        # find index into full grid at which reduced grid is starting
+        start = red_grid[0]
+        ind = np.argmax(np.logical_and(full_grid >= start - tol, full_grid <= start + tol))
+        # compare
+        assert np.allclose(red_grid, full_grid[ind : ind + len(red_grid)])
diff --git a/tests/test_components/test_time_modulation.py b/tests/test_components/test_time_modulation.py
index 43b88b940..6fc7aedeb 100644
--- a/tests/test_components/test_time_modulation.py
+++ b/tests/test_components/test_time_modulation.py
@@ -32,6 +32,74 @@
 # medium modulation spec
 MODULATION_SPEC = td.ModulationSpec()
 
+SUBSECTION = td.Box(size=(0.3, 0.4, 0.35), center=(0.4, 0.4, 0.4))
+
+
+def reduce(obj):
+
+    return obj.sel_inside(SUBSECTION.bounds)
+
+
+def check_reduction(obj, obj_reduced):
+
+    for field in ["amplitude", "phase"]:
+        original = getattr(obj, field)
+        reduced = getattr(obj_reduced, field)
+
+        if isinstance(original, float):
+            assert reduced == original
+            continue
+
+        diff = (original - reduced).abs
+        assert diff.does_cover(SUBSECTION.bounds)
+        assert np.allclose(diff, 0)
+
+
+def check_sp_reduction(sp):
+
+    check_reduction(sp, reduce(sp))
+
+
+def check_st_reduction(st):
+
+    check_reduction(st.space_modulation, reduce(st).space_modulation)
+
+
+def check_med_reduction(med):
+
+    med_red = reduce(med)
+
+    for field in ["permittivity", "conductivity"]:
+
+        field_mod = getattr(med.modulation_spec, field)
+        field_mod_red = getattr(med_red.modulation_spec, field)
+        if field_mod is None:
+            assert field_mod_red is None
+        else:
+            check_reduction(field_mod.space_modulation, field_mod_red.space_modulation)
+
+
+def check_ani_med_reduction(med):
+
+    reduced_med = reduce(med)
+
+    for comp, comp_red in zip(
+        [med.xx, med.yy, med.zz], [reduced_med.xx, reduced_med.yy, reduced_med.zz]
+    ):
+
+        if comp.modulation_spec is None:
+            assert comp_red.modulation_spec is None
+        else:
+
+            for field in ["permittivity", "conductivity"]:
+
+                field_mod = getattr(comp.modulation_spec, field)
+                field_mod_red = getattr(comp_red.modulation_spec, field)
+                if field_mod is None:
+                    assert field_mod_red is None
+                else:
+                    check_reduction(field_mod.space_modulation, field_mod_red.space_modulation)
+
 
 def test_time_modulation():
     """time modulation: only supporting CW for now."""
@@ -47,6 +115,7 @@ def test_space_modulation():
     """uniform or custom space modulation"""
     # uniform in both amplitude and phase
     assert isclose(SP_UNIFORM.max_modulation, 1)
+    check_sp_reduction(SP_UNIFORM)
 
     # uniform in phase, but custom in amplitude
     with pytest.raises(pydantic.ValidationError):
@@ -54,35 +123,42 @@ def test_space_modulation():
 
     sp = SP_UNIFORM.updated_copy(amplitude=ARRAY)
     assert isclose(sp.max_modulation, np.max(ARRAY))
+    check_sp_reduction(sp)
 
     # uniform in amplitude, but custom in phase
     with pytest.raises(pydantic.ValidationError):
         sp = SP_UNIFORM.updated_copy(phase=ARRAY_CMP)
     sp = SP_UNIFORM.updated_copy(phase=ARRAY)
     assert isclose(sp.max_modulation, 1)
+    check_sp_reduction(sp)
 
     # custom in both
     with pytest.raises(pydantic.ValidationError):
         sp = SP_UNIFORM.updated_copy(phase=ARRAY_CMP, amplitude=ARRAY_CMP)
     sp = SP_UNIFORM.updated_copy(phase=ARRAY, amplitude=ARRAY)
+    check_sp_reduction(sp)
 
 
 def test_space_time_modulation():
     # cw modulation, uniform in space
     assert isclose(ST.max_modulation, AMP_TIME)
     assert not ST.negligible_modulation
+    check_st_reduction(ST)
 
     # cw modulation, but 0 amplitude
     st = ST.updated_copy(time_modulation=CW.updated_copy(amplitude=0))
     assert st.negligible_modulation
+    check_st_reduction(st)
 
     st = ST.updated_copy(space_modulation=td.SpaceModulation(amplitude=0))
     assert st.negligible_modulation
+    check_st_reduction(st)
 
     # cw modulation, nonuniform in space
     st = ST.updated_copy(space_modulation=td.SpaceModulation(amplitude=ARRAY, phase=ARRAY))
     assert not st.negligible_modulation
     assert isclose(st.max_modulation, AMP_TIME * np.max(ARRAY))
+    check_st_reduction(st)
 
 
 def test_modulated_medium():
@@ -91,10 +167,12 @@ def test_modulated_medium():
     medium = td.Medium()
     assert medium.modulation_spec is None
     assert medium.time_modulated == False
+    reduce(medium)
 
     assert MODULATION_SPEC.applied_modulation == False
     medium = medium.updated_copy(modulation_spec=MODULATION_SPEC)
     assert medium.time_modulated == False
+    reduce(medium)
 
     # permittivity modulated
     modulation_spec = MODULATION_SPEC.updated_copy(permittivity=ST)
@@ -103,6 +181,7 @@ def test_modulated_medium():
         medium = td.Medium(modulation_spec=modulation_spec)
     medium = td.Medium(permittivity=2, modulation_spec=modulation_spec)
     assert isclose(medium.n_cfl, np.sqrt(2 - AMP_TIME))
+    check_med_reduction(medium)
 
     # conductivity modulated
     modulation_spec = MODULATION_SPEC.updated_copy(conductivity=ST)
@@ -111,6 +190,8 @@ def test_modulated_medium():
         medium = td.Medium(modulation_spec=modulation_spec)
     medium_sometimes_active = td.Medium(modulation_spec=modulation_spec, allow_gain=True)
     medium = td.Medium(conductivity=2, modulation_spec=modulation_spec)
+    check_med_reduction(medium)
+    check_med_reduction(medium_sometimes_active)
 
     # both modulated, but different time modulation: error
     st_freq2 = ST.updated_copy(
@@ -126,6 +207,7 @@ def test_modulated_medium():
         conductivity=1,
         modulation_spec=modulation_spec,
     )
+    check_med_reduction(medium)
 
 
 def test_unsupported_modulated_medium_types():
@@ -178,6 +260,8 @@ def test_supported_modulated_medium_types():
     with pytest.raises(pydantic.ValidationError):
         mat = mat_p.updated_copy(modulation_spec=modulation_both_spec)
     mat = mat_p.updated_copy(modulation_spec=modulation_both_spec, allow_gain=True)
+    check_med_reduction(mat)
+    check_med_reduction(mat_p)
 
     # custom
     permittivity = td.SpatialDataArray(np.ones((1, 1, 1)) * 2, coords=dict(x=[1], y=[1], z=[1]))
@@ -191,14 +275,18 @@ def test_supported_modulated_medium_types():
     with pytest.raises(pydantic.ValidationError):
         mat = mat_c.updated_copy(modulation_spec=modulation_both_spec)
     mat = mat_c.updated_copy(modulation_spec=modulation_both_spec, allow_gain=True)
+    check_med_reduction(mat_c)
+    check_med_reduction(mat)
 
     # anisotropic medium component
     mat = td.AnisotropicMedium(xx=td.Medium(), yy=mat_p, zz=td.Medium())
     assert mat.time_modulated
     assert isclose(mat.n_cfl, np.sqrt(2 - AMP_TIME))
+    check_ani_med_reduction(mat)
 
     # custom anistropic medium component
     mat_uc = td.CustomMedium(permittivity=permittivity)
     mat = td.CustomAnisotropicMedium(xx=mat_uc, yy=mat_c, zz=mat_uc)
     assert mat.time_modulated
     assert isclose(mat.n_cfl, np.sqrt(2 - AMP_TIME))
+    check_ani_med_reduction(mat)
diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py
index 8f9e7f846..bfda830fe 100644
--- a/tests/test_plugins/test_mode_solver.py
+++ b/tests/test_plugins/test_mode_solver.py
@@ -214,6 +214,18 @@ def verify_dtype(ms):
         assert dtype == field.dtype
 
 
+def check_ms_reduction(ms):
+
+    ms_red = ms.reduced_simulation_copy
+    grids_1d = ms._solver_grid.boundaries
+    grids_1d_red = ms_red._solver_grid.boundaries
+    assert np.allclose(grids_1d.x, grids_1d_red.x)
+    assert np.allclose(grids_1d.y, grids_1d_red.y)
+    assert np.allclose(grids_1d.z, grids_1d_red.z)
+    modes_red = ms.solve()
+    assert np.allclose(ms.data.n_eff.values, modes_red.n_eff.values)
+
+
 def test_mode_solver_validation():
     """Test invalidate mode solver setups."""
 
@@ -319,6 +331,7 @@ def test_mode_solver_simple(mock_remote_api, local):
         verify_pol_fraction(ms)
         verify_dtype(ms)
         dataframe = ms.data.to_dataframe()
+        check_ms_reduction(ms)
 
     else:
         _ = msweb.run(ms)
@@ -379,6 +392,9 @@ def test_mode_solver_custom_medium(mock_remote_api, local, tmp_path):
         modes = ms.solve() if local else msweb.run(ms)
         n_eff.append(modes.n_eff.values)
 
+        if local:
+            check_ms_reduction(ms)
+
         fname = str(tmp_path / "ms_custom_medium.hdf5")
         ms.to_file(fname)
         m2 = ModeSolver.from_file(fname)
@@ -445,6 +461,9 @@ def test_mode_solver_straight_vs_angled():
         direction="-",
     )
 
+    check_ms_reduction(ms)
+    check_ms_reduction(ms_angled)
+
     for key, val in ms.data.modes_info.items():
         tol = 1e-2
         if key == "mode area":
@@ -483,6 +502,7 @@ def test_mode_solver_angle_bend():
     verify_pol_fraction(ms)
     verify_dtype(ms)
     dataframe = ms.data.to_dataframe()
+    check_ms_reduction(ms)
 
     # Plot field
     _, ax = plt.subplots(1)
@@ -519,6 +539,7 @@ def test_mode_solver_2D():
     verify_pol_fraction(ms)
     verify_dtype(ms)
     dataframe = ms.data.to_dataframe()
+    check_ms_reduction(ms)
 
     mode_spec = td.ModeSpec(
         num_modes=3,
@@ -540,6 +561,7 @@ def test_mode_solver_2D():
     compare_colocation(ms)
     # verify_pol_fraction(ms)
     dataframe = ms.data.to_dataframe()
+    check_ms_reduction(ms)
 
     # The simulation and the mode plane are both 0D along the same dimension
     simulation = td.Simulation(
@@ -552,6 +574,7 @@ def test_mode_solver_2D():
     ms = ModeSolver(simulation=simulation, plane=PLANE, mode_spec=mode_spec, freqs=[td.C_0 / 1.0])
     compare_colocation(ms)
     verify_pol_fraction(ms)
+    check_ms_reduction(ms)
 
 
 @pytest.mark.parametrize("local", [True, False])
@@ -605,6 +628,7 @@ def test_group_index(mock_remote_api, log_capture, local):
         assert len(log_capture) == 2
         _ = modes.dispersion
         assert len(log_capture) == 2
+        check_ms_reduction(ms)
 
     # Group index calculated
     ms = ModeSolver(
@@ -623,6 +647,7 @@ def test_group_index(mock_remote_api, log_capture, local):
         assert (modes.dispersion.sel(mode_index=0).values < 1500).all()
         assert (modes.dispersion.sel(mode_index=1).values > -16500).all()
         assert (modes.dispersion.sel(mode_index=1).values < -15000).all()
+        check_ms_reduction(ms)
 
 
 def test_pml_params():
@@ -685,6 +710,7 @@ def test_mode_solver_nan_pol_fraction():
     )
 
     md = ms.solve()
+    check_ms_reduction(ms)
 
     assert list(np.where(np.isnan(md.pol_fraction.te))[1]) == [8, 9]
 
diff --git a/tidy3d/components/data/data_array.py b/tidy3d/components/data/data_array.py
index 683396aeb..18fab26f4 100644
--- a/tidy3d/components/data/data_array.py
+++ b/tidy3d/components/data/data_array.py
@@ -1,6 +1,7 @@
 """Storing tidy3d data at it's most fundamental level as xr.DataArray objects"""
 from __future__ import annotations
 from typing import Dict, List
+from abc import ABC
 
 import xarray as xr
 import numpy as np
@@ -240,17 +241,8 @@ class MixedModeDataArray(DataArray):
     _dims = ("f", "mode_index_0", "mode_index_1")
 
 
-class SpatialDataArray(DataArray):
-    """Spatial distribution.
-
-    Example
-    -------
-    >>> x = [1,2]
-    >>> y = [2,3,4]
-    >>> z = [3,4,5,6]
-    >>> coords = dict(x=x, y=y, z=z)
-    >>> fd = SpatialDataArray((1+1j) * np.random.random((2,3,4)), coords=coords)
-    """
+class AbstractSpatialDataArray(DataArray, ABC):
+    """Spatial distribution."""
 
     __slots__ = ()
     _dims = ("x", "y", "z")
@@ -274,33 +266,40 @@ def sel_inside(self, bounds: Bound) -> SpatialDataArray:
 
         inds_list = []
 
-        for coord, smin, smax in zip(self.coords.values(), bounds[0], bounds[1]):
+        coords = (self.x, self.y, self.z)
+
+        for coord, smin, smax in zip(coords, bounds[0], bounds[1]):
 
             length = len(coord)
 
-            # if data does not cover structure at all take the closest index
-            if smax < coord[0]:  # structure is completely on the left side
+            # one point along direction, assume invariance
+            if length == 1:
+                comp_inds = [0]
+            else:
+
+                # if data does not cover structure at all take the closest index
+                if smax < coord[0]:  # structure is completely on the left side
 
-                # take 2 if possible, so that linear interpolation is possible
-                comp_inds = np.arange(0, max(2, length))
+                    # take 2 if possible, so that linear iterpolation is possible
+                    comp_inds = np.arange(0, max(2, length))
 
-            elif smin > coord[-1]:  # structure is completely on the right side
+                elif smin > coord[-1]:  # structure is completely on the right side
 
-                # take 2 if possible, so that linear interpolation is possible
-                comp_inds = np.arange(min(0, length - 2), length)
+                    # take 2 if possible, so that linear iterpolation is possible
+                    comp_inds = np.arange(min(0, length - 2), length)
 
-            else:
-                if smin < coord[0]:
-                    ind_min = 0
                 else:
-                    ind_min = max(0, (coord >= smin).argmax().data - 1)
+                    if smin < coord[0]:
+                        ind_min = 0
+                    else:
+                        ind_min = max(0, (coord >= smin).argmax().data - 1)
 
-                if smax > coord[-1]:
-                    ind_max = length - 1
-                else:
-                    ind_max = (coord >= smax).argmax().data
+                    if smax > coord[-1]:
+                        ind_max = length - 1
+                    else:
+                        ind_max = (coord >= smax).argmax().data
 
-                comp_inds = np.arange(ind_min, ind_max + 1)
+                    comp_inds = np.arange(ind_min, ind_max + 1)
 
             inds_list.append(comp_inds)
 
@@ -323,11 +322,27 @@ def does_cover(self, bounds: Bound) -> bool:
             Full cover check outcome.
         """
 
+        coords = (self.x, self.y, self.z)
         return all(
             (coord[0] <= smin and coord[-1] >= smax) or len(coord) == 1
-            for coord, smin, smax in zip(self.coords.values(), bounds[0], bounds[1])
+            for coord, smin, smax in zip(coords, bounds[0], bounds[1])
         )
 
+
+class SpatialDataArray(AbstractSpatialDataArray):
+    """Spatial distribution.
+
+    Example
+    -------
+    >>> x = [1,2]
+    >>> y = [2,3,4]
+    >>> z = [3,4,5,6]
+    >>> coords = dict(x=x, y=y, z=z)
+    >>> fd = SpatialDataArray((1+1j) * np.random.random((2,3,4)), coords=coords)
+    """
+
+    __slots__ = ()
+
     def reflect(self, axis: Axis, center: float) -> SpatialDataArray:
         """Reflect data across the plane define by parameters ``axis`` and ``center`` from right to
         left.
@@ -380,7 +395,7 @@ def reflect(self, axis: Axis, center: float) -> SpatialDataArray:
         return SpatialDataArray(new_data, coords=coords_dict)
 
 
-class ScalarFieldDataArray(DataArray):
+class ScalarFieldDataArray(AbstractSpatialDataArray):
     """Spatial distribution in the frequency-domain.
 
     Example
@@ -395,10 +410,9 @@ class ScalarFieldDataArray(DataArray):
 
     __slots__ = ()
     _dims = ("x", "y", "z", "f")
-    _data_attrs = {"long_name": "field value"}
 
 
-class ScalarFieldTimeDataArray(DataArray):
+class ScalarFieldTimeDataArray(AbstractSpatialDataArray):
     """Spatial distribution in the time-domain.
 
     Example
@@ -413,10 +427,9 @@ class ScalarFieldTimeDataArray(DataArray):
 
     __slots__ = ()
     _dims = ("x", "y", "z", "t")
-    _data_attrs = {"long_name": "field value"}
 
 
-class ScalarModeFieldDataArray(DataArray):
+class ScalarModeFieldDataArray(AbstractSpatialDataArray):
     """Spatial distribution of a mode in frequency-domain as a function of mode index.
 
     Example
@@ -432,7 +445,6 @@ class ScalarModeFieldDataArray(DataArray):
 
     __slots__ = ()
     _dims = ("x", "y", "z", "f", "mode_index")
-    _data_attrs = {"long_name": "field value"}
 
 
 class FluxDataArray(DataArray):
diff --git a/tidy3d/components/medium.py b/tidy3d/components/medium.py
index 327d2fbf5..7e22941d4 100644
--- a/tidy3d/components/medium.py
+++ b/tidy3d/components/medium.py
@@ -816,6 +816,28 @@ def is_pec(self):
         """Whether the medium is a PEC."""
         return False
 
+    def sel_inside(self, bounds: Bound) -> AbstractMedium:
+        """Return a new medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        AbstractMedium
+            Medium with reduced data.
+        """
+
+        if self.modulation_spec is not None:
+            modulation_reduced = self.modulation_spec.sel_inside(bounds)
+            return self.updated_copy(modulation_spec=modulation_reduced)
+
+        return self
+
 
 class AbstractCustomMedium(AbstractMedium, ABC):
     """A spatially varying medium."""
@@ -960,6 +982,31 @@ def _validate_isreal_dataarray_tuple(dataarray_tuple: Tuple[SpatialDataArray, ..
         """Validate that the dataarray is real"""
         return np.all([AbstractCustomMedium._validate_isreal_dataarray(f) for f in dataarray_tuple])
 
+    @abstractmethod
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new medium that contains the minimal amount custom data necessary to cover
+        a spatial region defined by ``bounds``."""
+
+    def sel_inside(self, bounds: Bound) -> AbstractCustomMedium:
+        """Return a new medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        AbstractMedium
+            Medium with reduced data.
+        """
+
+        self_mod_data_reduced = super().sel_inside(bounds)
+
+        return self_mod_data_reduced._sel_custom_data_inside(bounds)
+
 
 """ Dispersionless Medium """
 
@@ -1261,6 +1308,39 @@ def eps_dataarray_freq(
         eps = self.eps_sigma_to_eps_complex(self.permittivity, conductivity, frequency)
         return (eps, eps, eps)
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomMedium
+            CustomMedium with reduced data.
+        """
+        if not self.permittivity.does_cover(bounds=bounds):
+            log.warning(
+                "Permittivity spatial data array does not fully cover the requested region."
+            )
+        perm_reduced = self.permittivity.sel_inside(bounds=bounds)
+        cond_reduced = None
+        if self.conductivity is not None:
+            if not self.conductivity.does_cover(bounds=bounds):
+                log.warning(
+                    "Conductivity spatial data array does not fully cover the requested region."
+                )
+            cond_reduced = self.conductivity.sel_inside(bounds=bounds)
+
+        return self.updated_copy(
+            permittivity=perm_reduced,
+            conductivity=cond_reduced,
+        )
+
 
 class CustomMedium(AbstractCustomMedium):
     """:class:`.Medium` with user-supplied permittivity distribution.
@@ -1844,6 +1924,55 @@ def make_bound_coords(coords: np.ndarray, pt_min: float, pt_max: float) -> List[
 
         return grids
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomMedium
+            CustomMedium with reduced data.
+        """
+
+        perm_reduced = None
+        if self.permittivity is not None:
+            if not self.permittivity.does_cover(bounds=bounds):
+                log.warning(
+                    "Permittivity spatial data array does not fully cover the requested region."
+                )
+            perm_reduced = self.permittivity.sel_inside(bounds=bounds)
+
+        cond_reduced = None
+        if self.conductivity is not None:
+            if not self.conductivity.does_cover(bounds=bounds):
+                log.warning(
+                    "Conductivity spatial data array does not fully cover the requested region."
+                )
+            cond_reduced = self.conductivity.sel_inside(bounds=bounds)
+
+        eps_reduced = None
+        if self.eps_dataset is not None:
+            eps_reduced_dict = {}
+            for key, comp in self.eps_dataset.field_components.items():
+                if not comp.does_cover(bounds=bounds):
+                    log.warning(
+                        f"{key} spatial data array does not fully cover the requested region."
+                    )
+                eps_reduced_dict[key] = comp.sel_inside(bounds=bounds)
+            eps_reduced = PermittivityDataset(**eps_reduced_dict)
+
+        return self.updated_copy(
+            permittivity=perm_reduced,
+            conductivity=cond_reduced,
+            eps_dataset=eps_reduced,
+        )
+
 
 """ Dispersive Media """
 
@@ -2603,6 +2732,37 @@ def loss_upper_bound(self) -> float:
         """Not implemented yet."""
         raise SetupError("To be implemented.")
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomPoleResidue
+            CustomPoleResidue with reduced data.
+        """
+        if not self.eps_inf.does_cover(bounds=bounds):
+            log.warning("eps_inf spatial data array does not fully cover the requested region.")
+        eps_inf_reduced = self.eps_inf.sel_inside(bounds=bounds)
+        poles_reduced = []
+        for pole, residue in self.poles:
+
+            if not pole.does_cover(bounds=bounds):
+                log.warning("Pole spatial data array does not fully cover the requested region.")
+
+            if not residue.does_cover(bounds=bounds):
+                log.warning("Residue spatial data array does not fully cover the requested region.")
+
+            poles_reduced.append((pole.sel_inside(bounds), residue.sel_inside(bounds)))
+
+        return self.updated_copy(eps_inf=eps_inf_reduced, poles=poles_reduced)
+
 
 class Sellmeier(DispersiveMedium):
     """A dispersive medium described by the Sellmeier model.
@@ -2891,6 +3051,38 @@ def from_dispersion(
             **kwargs,
         )
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomSellmeier
+            CustomSellmeier with reduced data.
+        """
+        coeffs_reduced = []
+        for b_coeff, c_coeff in self.coeffs:
+
+            if not b_coeff.does_cover(bounds=bounds):
+                log.warning(
+                    "Sellmeier B coeff spatial data array does not fully cover the requested region."
+                )
+
+            if not c_coeff.does_cover(bounds=bounds):
+                log.warning(
+                    "Sellmeier C coeff spatial data array does not fully cover the requested region."
+                )
+
+            coeffs_reduced.append((b_coeff.sel_inside(bounds), c_coeff.sel_inside(bounds)))
+
+        return self.updated_copy(coeffs=coeffs_reduced)
+
 
 class Lorentz(DispersiveMedium):
     """A dispersive medium described by the Lorentz model.
@@ -3198,6 +3390,48 @@ def eps_dataarray_freq(
         eps = Lorentz.eps_model(self, frequency)
         return (eps, eps, eps)
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomLorentz
+            CustomLorentz with reduced data.
+        """
+        if not self.eps_inf.does_cover(bounds=bounds):
+            log.warning("Eps inf spatial data array does not fully cover the requested region.")
+        eps_inf_reduced = self.eps_inf.sel_inside(bounds=bounds)
+        coeffs_reduced = []
+        for de, f, delta in self.coeffs:
+
+            if not de.does_cover(bounds=bounds):
+                log.warning(
+                    "Lorentz 'de' spatial data array does not fully cover the requested region."
+                )
+
+            if not f.does_cover(bounds=bounds):
+                log.warning(
+                    "Lorentz 'f' spatial data array does not fully cover the requested region."
+                )
+
+            if not delta.does_cover(bounds=bounds):
+                log.warning(
+                    "Lorentz 'delta' spatial data array does not fully cover the requested region."
+                )
+
+            coeffs_reduced.append(
+                (de.sel_inside(bounds), f.sel_inside(bounds), delta.sel_inside(bounds))
+            )
+
+        return self.updated_copy(eps_inf=eps_inf_reduced, coeffs=coeffs_reduced)
+
 
 class Drude(DispersiveMedium):
     """A dispersive medium described by the Drude model.
@@ -3385,6 +3619,41 @@ def eps_dataarray_freq(
         eps = Drude.eps_model(self, frequency)
         return (eps, eps, eps)
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomDrude
+            CustomDrude with reduced data.
+        """
+        if not self.eps_inf.does_cover(bounds=bounds):
+            log.warning("Eps inf spatial data array does not fully cover the requested region.")
+        eps_inf_reduced = self.eps_inf.sel_inside(bounds=bounds)
+        coeffs_reduced = []
+        for f, delta in self.coeffs:
+
+            if not f.does_cover(bounds=bounds):
+                log.warning(
+                    "Drude 'f' spatial data array does not fully cover the requested region."
+                )
+
+            if not delta.does_cover(bounds=bounds):
+                log.warning(
+                    "Drude 'delta' spatial data array does not fully cover the requested region."
+                )
+
+            coeffs_reduced.append((f.sel_inside(bounds), delta.sel_inside(bounds)))
+
+        return self.updated_copy(eps_inf=eps_inf_reduced, coeffs=coeffs_reduced)
+
 
 class Debye(DispersiveMedium):
     """A dispersive medium described by the Debye model.
@@ -3590,6 +3859,41 @@ def eps_dataarray_freq(
         eps = Debye.eps_model(self, frequency)
         return (eps, eps, eps)
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        """Return a new custom medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        CustomDebye
+            CustomDebye with reduced data.
+        """
+        if not self.eps_inf.does_cover(bounds=bounds):
+            log.warning("Eps inf spatial data array does not fully cover the requested region.")
+        eps_inf_reduced = self.eps_inf.sel_inside(bounds=bounds)
+        coeffs_reduced = []
+        for de, tau in self.coeffs:
+
+            if not de.does_cover(bounds=bounds):
+                log.warning(
+                    "Debye 'f' spatial data array does not fully cover the requested region."
+                )
+
+            if not tau.does_cover(bounds=bounds):
+                log.warning(
+                    "Debye 'tau' spatial data array does not fully cover the requested region."
+                )
+
+            coeffs_reduced.append((de.sel_inside(bounds), tau.sel_inside(bounds)))
+
+        return self.updated_copy(eps_inf=eps_inf_reduced, coeffs=coeffs_reduced)
+
 
 IsotropicUniformMediumType = Union[Medium, PoleResidue, Sellmeier, Lorentz, Debye, Drude, PECMedium]
 IsotropicCustomMediumType = Union[
@@ -3771,6 +4075,26 @@ def is_comp_pec(self, comp: Axis):
         """Whether the medium is a PEC."""
         return isinstance(self.components[["xx", "yy", "zz"][comp]], PECMedium)
 
+    def sel_inside(self, bounds: Bound):
+        """Return a new medium that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        AnisotropicMedium
+            AnisotropicMedium with reduced data.
+        """
+
+        new_comps = [comp.sel_inside(bounds) for comp in [self.xx, self.yy, self.zz]]
+
+        return self.updated_copy(**dict(zip(["xx", "yy", "zz"], new_comps)))
+
 
 class FullyAnisotropicMedium(AbstractMedium):
     """Fully anisotropic medium including all 9 components of the permittivity and conductivity
@@ -4176,6 +4500,9 @@ def eps_dataarray_freq(
             for ind, mat_component in enumerate(self.components.values())
         )
 
+    def _sel_custom_data_inside(self, bounds: Bound):
+        return self
+
 
 class CustomAnisotropicMediumInternal(CustomAnisotropicMedium):
     """Diagonally anisotropic medium with spatially varying permittivity in each component.
diff --git a/tidy3d/components/monitor.py b/tidy3d/components/monitor.py
index 47fe32785..b9a651b07 100644
--- a/tidy3d/components/monitor.py
+++ b/tidy3d/components/monitor.py
@@ -983,7 +983,7 @@ class FieldProjectionAngleMonitor(AbstractFieldProjectionMonitor):
     ...     name='n2f_monitor',
     ...     custom_origin=(1,2,3),
     ...     phi=[0, np.pi/2],
-    ...     theta=np.linspace(-np.pi/2, np.pi/2, 100)
+    ...     theta=np.linspace(-np.pi/2, np.pi/2, 100),
     ...     far_field_approx=True,
     ...     )
 
diff --git a/tidy3d/components/scene.py b/tidy3d/components/scene.py
index d3ac88156..69b53031c 100644
--- a/tidy3d/components/scene.py
+++ b/tidy3d/components/scene.py
@@ -1274,7 +1274,6 @@ def perturbed_mediums_copy(
 
         scene_dict = self.dict()
         structures = self.structures
-        scene_bounds = self.bounds
         array_dict = {
             "temperature": temperature,
             "electron_density": electron_density,
@@ -1289,12 +1288,7 @@ def perturbed_mediums_copy(
             med = structure.medium
             if isinstance(med, AbstractPerturbationMedium):
                 # get structure's bounding box
-                s_bounds = structure.geometry.bounds
-
-                bounds = [
-                    np.max([scene_bounds[0], s_bounds[0]], axis=0),
-                    np.min([scene_bounds[1], s_bounds[1]], axis=0),
-                ]
+                bounds = structure.geometry.bounds
 
                 # for each structure select a minimal subset of data that covers it
                 restricted_arrays = {}
@@ -1321,22 +1315,6 @@ def perturbed_mediums_copy(
         med = self.medium
         if isinstance(med, AbstractPerturbationMedium):
 
-            # get scene's bounding box
-            bounds = scene_bounds
-
-            # for each structure select a minimal subset of data that covers it
-            restricted_arrays = {}
-
-            for name, array in array_dict.items():
-                if array is not None:
-                    restricted_arrays[name] = array.sel_inside(bounds)
-
-                    # check provided data fully cover scene
-                    if not array.does_cover(bounds):
-                        log.warning(f"Provided '{name}' does not fully cover scene domain.")
-
-            scene_dict["medium"] = med.perturbed_copy(
-                **restricted_arrays, interp_method=interp_method
-            )
+            scene_dict["medium"] = med.perturbed_copy(**array_dict, interp_method=interp_method)
 
         return Scene.parse_obj(scene_dict)
diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py
index 9d0501c1a..2dd404b83 100644
--- a/tidy3d/components/simulation.py
+++ b/tidy3d/components/simulation.py
@@ -7,6 +7,7 @@
 import numpy as np
 import xarray as xr
 import matplotlib as mpl
+import math
 
 from .base import cached_property
 from .validators import assert_objects_in_sim_bounds
@@ -16,9 +17,10 @@
 from .geometry.mesh import TriangleMesh
 from .geometry.polyslab import PolySlab
 from .geometry.utils import flatten_groups, traverse_geometries
-from .types import Ax, FreqBound, Axis, annotate_type, InterpMethod
+from .types import Ax, FreqBound, Axis, annotate_type, InterpMethod, Symmetry
+from .types import Literal
 from .grid.grid import Coords1D, Grid, Coords
-from .grid.grid_spec import GridSpec, UniformGrid, AutoGrid
+from .grid.grid_spec import GridSpec, UniformGrid, AutoGrid, CustomGrid
 from .medium import Medium, MediumType, AbstractMedium
 from .medium import AbstractCustomMedium, Medium2D, MediumType3D
 from .medium import AnisotropicMedium, FullyAnisotropicMedium, AbstractPerturbationMedium
@@ -1233,7 +1235,10 @@ def diffraction_monitor_medium(cls, val, values):
             if isinstance(monitor, DiffractionMonitor):
                 medium_set = Scene.intersecting_media(monitor, structures)
                 medium = medium_set.pop() if medium_set else medium
-                _, index_k = medium.nk_model(frequency=np.array(monitor.freqs))
+                freqs = np.array(monitor.freqs)
+                if isinstance(medium, AbstractCustomMedium) and len(freqs) > 1:
+                    freqs = 0.5 * (np.min(freqs) + np.max(freqs))
+                _, index_k = medium.nk_model(frequency=freqs)
                 if not np.all(index_k == 0):
                     raise SetupError("Diffraction monitors must not lie in a lossy medium.")
         return val
@@ -3419,7 +3424,7 @@ def from_scene(cls, scene: Scene, **kwargs) -> Simulation:
         scene : :class:.`Scene`
             Size of object in x, y, and z directions.
         **kwargs
-            Other arguments
+            Other arguments passed to new simulation instance.
 
         Example
         -------
@@ -3445,3 +3450,201 @@ def from_scene(cls, scene: Scene, **kwargs) -> Simulation:
             medium=scene.medium,
             **kwargs,
         )
+
+    def subsection(
+        self,
+        region: Box,
+        boundary_spec: BoundarySpec = None,
+        grid_spec: Union[GridSpec, Literal["identical"]] = None,
+        symmetry: Tuple[Symmetry, Symmetry, Symmetry] = None,
+        sources: Tuple[SourceType, ...] = None,
+        monitors: Tuple[MonitorType, ...] = None,
+        remove_outside_structures: bool = True,
+        remove_outside_custom_mediums: bool = False,
+        **kwargs,
+    ) -> Simulation:
+        """Generate a simulation instance containing only the ``region``.
+
+        Parameters
+        ----------
+        region : :class:.`Box`
+            New simulation domain.
+        boundary_spec : :class:.`BoundarySpec` = None
+            New boundary specification. If ``None``, then it is inherited from the original
+            simulation.
+        grid_spec : :class:.`GridSpec` = None
+            New grid specification. If ``None``, then it is inherited from the original
+            simulation. If ``identical``, then the original grid is transferred directly as a
+            :class:.`CustomGrid`. Note that in the latter case the region of the new simulation is
+            snapped to the original grid lines.
+        symmetry : Tuple[Literal[0, -1, 1], Literal[0, -1, 1], Literal[0, -1, 1]] = None
+            New simulation symmetry. If ``None``, then it is inherited from the original
+            simulation. Note that in this case the size and placement of new simulation domain
+            must be commensurate with the original symmetry.
+        sources : Tuple[SourceType, ...] = None
+            New list of sources. If ``None``, then the sources intersecting the new simulation
+            domain are inherited from the original simulation.
+        monitors : Tuple[MonitorType, ...] = None
+            New list of monitors. If ``None``, then the monitors intersecting the new simulation
+            domain are inherited from the original simulation.
+        remove_outside_structures : bool = True
+            Remove structures outside of the new simulation domain.
+        remove_outside_custom_mediums : bool = True
+            Remove custom medium data outside of the new simulation domain.
+        **kwargs
+            Other arguments passed to new simulation instance.
+        """
+
+        # must intersect the original domain
+        if not self.intersects(region):
+            raise SetupError("Requested region does not intersect simulation domain")
+
+        # restrict to the original simulation domain
+        new_bounds = Box.bounds_intersection(self.bounds, region.bounds)
+        new_bounds = [list(new_bounds[0]), list(new_bounds[1])]
+
+        # grid spec inheritace
+        if grid_spec is None:
+            grid_spec = self.grid_spec
+        elif isinstance(grid_spec, str) and grid_spec == "identical":
+            # create a custom grid from existing one
+            grids_1d = self.grid.boundaries
+            grid_spec = GridSpec(
+                grid_x=CustomGrid(dl=tuple(np.diff(grids_1d.x)), custom_offset=grids_1d.x[0]),
+                grid_y=CustomGrid(dl=tuple(np.diff(grids_1d.y)), custom_offset=grids_1d.y[0]),
+                grid_z=CustomGrid(dl=tuple(np.diff(grids_1d.z)), custom_offset=grids_1d.z[0]),
+            )
+
+            # adjust region bounds to perfectly coincide with the grid
+            # note, sometimes (when a box already seems to perfrecty align with the grid)
+            # this causes the new region to expand one more pixel because of numerical roundoffs
+            aux_box = Box.from_bounds(*new_bounds)
+            grid_inds = self.grid.discretize_inds(box=aux_box)
+
+            new_bounds = [
+                [
+                    grids_1d.x[grid_inds[0][0]],
+                    grids_1d.y[grid_inds[1][0]],
+                    grids_1d.z[grid_inds[2][0]],
+                ],
+                [
+                    grids_1d.x[grid_inds[0][1]],
+                    grids_1d.y[grid_inds[1][1]],
+                    grids_1d.z[grid_inds[2][1]],
+                ],
+            ]
+
+        # if symmetry is not overriden we inherit it from the original simulation where is needed
+        if symmetry is None:
+
+            # start with no symmetry
+            symmetry = [0, 0, 0]
+
+            # now check in each dimension whether we cross symmetry plane
+            for dim in range(3):
+                if self.symmetry[dim] != 0:
+                    crosses_symmetry = (
+                        new_bounds[0][dim] < self.center[dim]
+                        and new_bounds[1][dim] > self.center[dim]
+                    )
+
+                    # inherit symmetry only if we cross symmetry plane, otherwise we don't impose
+                    # symmetry even if the original simulation had symmetry
+                    if crosses_symmetry:
+                        symmetry[dim] = self.symmetry[dim]
+                        center = (new_bounds[0][dim] + new_bounds[1][dim]) / 2
+
+                        if not math.isclose(center, self.center[dim]):
+                            log.warning(
+                                f"The original simulation is symmetric along {'xyz'[dim]} direction. "
+                                "The requested new simulation region does cross the symmetry plane but is "
+                                "not symmetric with respect to it. To preserve correct symmetry, "
+                                "the requested simulation region is expanded symmetrically."
+                            )
+                            new_bounds[0][dim] = 2 * self.center[dim] - new_bounds[1][dim]
+
+        # symmetry and grid spec treatments could change new simulation bounds
+        # thus, recreate a box instance
+        new_box = Box.from_bounds(*new_bounds)
+
+        # inheritance of structures, sources, monitors, and boundary specs
+        if remove_outside_structures:
+            new_structures = [strc for strc in self.structures if new_box.intersects(strc.geometry)]
+        else:
+            new_structures = self.structures
+
+        if sources is None:
+            sources = [src for src in self.sources if new_box.intersects(src)]
+
+        if monitors is None:
+            monitors = [mnt for mnt in self.monitors if new_box.intersects(mnt)]
+
+        if boundary_spec is None:
+            boundary_spec = self.boundary_spec
+
+        # reduction of custom medium data
+        new_sim_medium = self.medium
+        if remove_outside_custom_mediums:
+
+            # check for special treatment in case of PML
+            if any(
+                any(isinstance(edge, (PML, StablePML, Absorber)) for edge in boundary)
+                for boundary in boundary_spec.to_list
+            ):
+                # if we need to cut out outside custom medium we have to be careful about PML/Absorber
+                # we should include data in PML so that there is no artificial reflection at PML boundaries
+
+                # to do this, we first create an auxiliary simulation
+                aux_sim = self.updated_copy(
+                    center=new_box.center,
+                    size=new_box.size,
+                    grid_spec=grid_spec,
+                    boundary_spec=boundary_spec,
+                    monitors=[],
+                    sources=sources,  # need wavelength in case of auto grid
+                    symmetry=symmetry,
+                    structures=new_structures,
+                )
+
+                # then use its bounds as region for data cut off
+                new_bounds = aux_sim.simulation_bounds
+
+                # Note that this is not a full proof strategy. For example, if grid_spec is AutoGrid
+                # then after outside custom medium data is removed the grid sizes and, thus,
+                # pml extents can change as well
+
+            # now cut out custom medium data
+            new_structures_reduced_data = []
+
+            for structure in new_structures:
+                medium = structure.medium
+                if isinstance(medium, AbstractCustomMedium):
+                    new_structure_bounds = Box.bounds_intersection(
+                        new_bounds, structure.geometry.bounds
+                    )
+                    new_medium = medium.sel_inside(bounds=new_structure_bounds)
+                    new_structure = structure.updated_copy(medium=new_medium)
+                    new_structures_reduced_data.append(new_structure)
+                else:
+                    new_structures_reduced_data.append(structure)
+
+            new_structures = new_structures_reduced_data
+
+            if isinstance(self.medium, AbstractCustomMedium):
+                new_sim_medium = self.medium.sel_inside(bounds=new_bounds)
+
+        # finally, create an updated copy with all modifications
+        new_sim = self.updated_copy(
+            center=new_box.center,
+            size=new_box.size,
+            medium=new_sim_medium,
+            grid_spec=grid_spec,
+            boundary_spec=boundary_spec,
+            monitors=monitors,
+            sources=sources,
+            symmetry=symmetry,
+            structures=new_structures,
+            **kwargs,
+        )
+
+        return new_sim
diff --git a/tidy3d/components/time_modulation.py b/tidy3d/components/time_modulation.py
index 3b2aa29fc..600c9a4fd 100644
--- a/tidy3d/components/time_modulation.py
+++ b/tidy3d/components/time_modulation.py
@@ -9,7 +9,7 @@
 import numpy as np
 
 from .base import Tidy3dBaseModel, cached_property
-from .types import InterpMethod
+from .types import InterpMethod, Bound
 from .time import AbstractTimeDependence
 from .data.data_array import SpatialDataArray
 from ..exceptions import ValidationError
@@ -164,6 +164,34 @@ def max_modulation(self) -> float:
         """Estimated maximum modulation amplitude."""
         return np.max(abs(np.array(self.amplitude)))
 
+    def sel_inside(self, bounds: Bound) -> SpaceModulation:
+        """Return a new space modulation that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        SpaceModulation
+            SpaceModulation with reduced data.
+        """
+
+        if isinstance(self.amplitude, SpatialDataArray):
+            amp_reduced = self.amplitude.sel_inside(bounds)
+        else:
+            amp_reduced = self.amplitude
+
+        if isinstance(self.phase, SpatialDataArray):
+            phase_reduced = self.phase.sel_inside(bounds)
+        else:
+            phase_reduced = self.phase
+
+        return self.updated_copy(amplitude=amp_reduced, phase=phase_reduced)
+
 
 SpaceModulationType = Union[SpaceModulation]
 
@@ -210,6 +238,24 @@ def negligible_modulation(self) -> bool:
             return True
         return False
 
+    def sel_inside(self, bounds: Bound) -> SpaceTimeModulation:
+        """Return a new space-time modulation that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        SpaceTimeModulation
+            SpaceTimeModulation with reduced data.
+        """
+
+        return self.updated_copy(space_modulation=self.space_modulation.sel_inside(bounds))
+
 
 class ModulationSpec(Tidy3dBaseModel):
     """Specification adding space-time modulation to the non-dispersive part of medium
@@ -245,3 +291,29 @@ def _same_modulation_frequency(cls, val, values):
     def applied_modulation(self) -> bool:
         """Check if any modulation has been applied to `permittivity` or `conductivity`."""
         return self.permittivity is not None or self.conductivity is not None
+
+    def sel_inside(self, bounds: Bound) -> ModulationSpec:
+        """Return a new modulation specficiation that contains the minimal amount data necessary to cover
+        a spatial region defined by ``bounds``.
+
+
+        Parameters
+        ----------
+        bounds : Tuple[float, float, float], Tuple[float, float float]
+            Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
+
+        Returns
+        -------
+        ModulationSpec
+            ModulationSpec with reduced data.
+        """
+
+        perm_reduced = None
+        if self.permittivity is not None:
+            perm_reduced = self.permittivity.sel_inside(bounds)
+
+        cond_reduced = None
+        if self.conductivity is not None:
+            cond_reduced = self.conductivity.sel_inside(bounds)
+
+        return self.updated_copy(permittivity=perm_reduced, conductivity=cond_reduced)
diff --git a/tidy3d/components/validators.py b/tidy3d/components/validators.py
index b7661d18f..2536f911d 100644
--- a/tidy3d/components/validators.py
+++ b/tidy3d/components/validators.py
@@ -171,7 +171,7 @@ def objects_in_sim_bounds(cls, val, values):
             if not sim_box.intersects(geometric_object.geometry):
 
                 message = (
-                    f"'{geometric_object}' (at `simulation.{field_name}[{position_index}]`) "
+                    f"'simulation.{field_name}[{position_index}]'"
                     "is completely outside of simulation domain."
                 )
                 custom_loc = [field_name, position_index]
diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py
index 68b634e1a..83a96901d 100644
--- a/tidy3d/plugins/mode/mode_solver.py
+++ b/tidy3d/plugins/mode/mode_solver.py
@@ -12,6 +12,7 @@
 
 from ...log import log
 from ...components.base import Tidy3dBaseModel, cached_property
+from ...components.boundary import PECBoundary, BoundarySpec, Boundary, PML, StablePML, Absorber
 from ...components.geometry.base import Box
 from ...components.simulation import Simulation
 from ...components.grid.grid import Grid
@@ -139,20 +140,33 @@ def solver_symmetry(self) -> Tuple[Symmetry, Symmetry]:
         _, solver_sym = self.plane.pop_axis(mode_symmetry, axis=self.normal_axis)
         return solver_sym
 
-    @cached_property
-    def _solver_grid(self) -> Grid:
-        """Grid for the mode solver, not snapped to plane or simulation zero dims, and also with
-        a small correction for symmetries. We don't do the snapping yet because 0-sized cells are
-        currently confusing to the subpixel averaging. The final data coordinates along the
-        plane normal dimension and dimensions where the simulation domain is 2D will be correctly
-        set after the solve."""
+    def _get_solver_grid(
+        self, preserve_layer_behind: bool = False, truncate_symmetry: bool = True
+    ) -> Grid:
+        """Grid for the mode solver, not snapped to plane or simulation zero dims, and optionally
+        corrected for symmetries.
+
+        Parameters
+        ----------
+        preserve_layer_behind : bool = False
+            Do not discard the layer of cells behind the main layer of cells. Together they
+            represent the region where custom medium data is needed for proper subpixel.
+        truncate_symmetry : bool = True
+            Truncate to symmetry quadrant if symmetry present.
+
+        Returns
+        -------
+        :class:.`Grid`
+            The resulting grid.
+        """
 
         monitor = self.to_mode_solver_monitor(name=MODE_MONITOR_NAME, colocate=False)
 
         span_inds = self.simulation._discretize_inds_monitor(monitor)
 
         # Remove extension along monitor normal
-        span_inds[self.normal_axis, 0] += 1
+        if not preserve_layer_behind:
+            span_inds[self.normal_axis, 0] += 1
         span_inds[self.normal_axis, 1] -= 1
 
         # Do not extend if simulation has a single pixel along a dimension
@@ -161,13 +175,24 @@ def _solver_grid(self) -> Grid:
                 span_inds[dim] = [0, 1]
 
         # Truncate to symmetry quadrant if symmetry present
-        _, plane_inds = monitor.pop_axis([0, 1, 2], self.normal_axis)
-        for dim, sym in enumerate(self.solver_symmetry):
-            if sym != 0:
-                span_inds[plane_inds[dim], 0] += np.diff(span_inds[plane_inds[dim]]) // 2
+        if truncate_symmetry:
+            _, plane_inds = Box.pop_axis([0, 1, 2], self.normal_axis)
+            for dim, sym in enumerate(self.solver_symmetry):
+                if sym != 0:
+                    span_inds[plane_inds[dim], 0] += np.diff(span_inds[plane_inds[dim]]) // 2
 
         return self.simulation._subgrid(span_inds=span_inds)
 
+    @cached_property
+    def _solver_grid(self) -> Grid:
+        """Grid for the mode solver, not snapped to plane or simulation zero dims, and also with
+        a small correction for symmetries. We don't do the snapping yet because 0-sized cells are
+        currently confusing to the subpixel averaging. The final data coordinates along the
+        plane normal dimension and dimensions where the simulation domain is 2D will be correctly
+        set after the solve."""
+
+        return self._get_solver_grid(preserve_layer_behind=False, truncate_symmetry=True)
+
     @cached_property
     def _num_cells_freqs_modes(self) -> Tuple[int, int, int]:
         """Get the number of spatial points, number of freqs, and number of modes requested."""
@@ -940,3 +965,48 @@ def _validate_modes_size(self):
 
     def validate_pre_upload(self, source_required: bool = True):
         self._validate_modes_size()
+
+    @cached_property
+    def reduced_simulation_copy(self):
+        """Strip objects not used by the mode solver from simulation object.
+        This might significantly reduce upload time in the presence of custom mediums.
+        """
+
+        # we preserve extracells along the normal direction to ensure
+        extended_grid = self._get_solver_grid(preserve_layer_behind=True, truncate_symmetry=False)
+        grids_1d = extended_grid.boundaries
+        new_sim_box = Box.from_bounds(
+            rmin=(grids_1d.x[0], grids_1d.y[0], grids_1d.z[0]),
+            rmax=(grids_1d.x[-1], grids_1d.y[-1], grids_1d.z[-1]),
+        )
+
+        # remove PML, Absorers, etc, to avoid unnecessary cells
+        bspec = self.simulation.boundary_spec
+
+        new_bspec_dict = {}
+        for axis in "xyz":
+            bcomp = bspec["x"]
+            for bside, sign in zip([bcomp.plus, bcomp.minus], "+-"):
+                if isinstance(bside, (PML, StablePML, Absorber)):
+                    new_bspec_dict[axis + sign] = PECBoundary()
+                else:
+                    new_bspec_dict[axis + sign] = bside
+
+        new_bspec = BoundarySpec(
+            x=Boundary(plus=new_bspec_dict["x+"], minus=new_bspec_dict["x-"]),
+            y=Boundary(plus=new_bspec_dict["y+"], minus=new_bspec_dict["y-"]),
+            z=Boundary(plus=new_bspec_dict["z+"], minus=new_bspec_dict["z-"]),
+        )
+
+        # extract sub-simulation removing everything irrelevant
+        new_sim = self.simulation.subsection(
+            region=new_sim_box,
+            monitors=[],
+            sources=[],
+            grid_spec="identical",
+            boundary_spec=new_bspec,
+            remove_outside_custom_mediums=True,
+            remove_outside_structures=True,
+        )
+
+        return self.updated_copy(simulation=new_sim)
diff --git a/tidy3d/web/api/mode.py b/tidy3d/web/api/mode.py
index d40013298..f82fa0685 100644
--- a/tidy3d/web/api/mode.py
+++ b/tidy3d/web/api/mode.py
@@ -46,6 +46,7 @@ def run(
     verbose: bool = True,
     progress_callback_upload: Callable[[float], None] = None,
     progress_callback_download: Callable[[float], None] = None,
+    reduce_simulation: bool = True,
 ) -> ModeSolverData:
     """Submits a :class:`.ModeSolver` to server, starts running, monitors progress, downloads,
     and loads results as a :class:`.ModeSolverData` object.
@@ -68,6 +69,9 @@ def run(
         Optional callback function called when uploading file with ``bytes_in_chunk`` as argument.
     progress_callback_download : Callable[[float], None] = None
         Optional callback function called when downloading file with ``bytes_in_chunk`` as argument.
+    reduce_simulation : bool = True
+        Restrict simulation to mode solver region. This can help reduce the amount of uploaded and
+        dowloaded data.
 
     Returns
     -------
@@ -79,6 +83,9 @@ def run(
     if verbose:
         console = get_logging_console()
 
+    if reduce_simulation:
+        mode_solver = mode_solver.reduced_simulation_copy
+
     task = ModeSolverTask.create(mode_solver, task_name, mode_solver_name, folder_name)
     if verbose:
         console.log(