Fix to mesher error in very specific cases

Warn when structures are too small compared to the generated mesh

Author: momchil-flex

CHANGELOG.md

@@ -8,10 +8,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - `ModeData.dispersion` and `ModeSolverData.dispersion` are calculated together with the group index.
 - String matching feature `contains_str` to `assert_log_level` testing utility.
+- Warning in automatic grid generation if a structure has a non-zero size along a given direction that is too small compared to a single mesh step.
 
 ### Changed
 - `jax` and `jaxlib` versions bumped to `0.4.*`.
 
+### Fixed
+- Error in automatic grid generation in specific cases with multiple thin structures.
+
 ## [2.5.0] - 2023-12-13
 
 ### Added

tests/test_components/test_meshgenerate.py

@@ -5,9 +5,10 @@
 
 import tidy3d as td
 from tidy3d.constants import fp_eps
-
 from tidy3d.components.grid.mesher import GradedMesher
 
+from ..utils import assert_log_level, log_capture
+
 np.random.seed(4)
 
 MESHER = GradedMesher()
@@ -650,6 +651,52 @@ def test_mesher_timeout():
     _ = sim.grid
 
 
+def test_small_structure_size(log_capture):
+    """Test that a warning is raised if a structure size is small during the auto meshing"""
+    box_size = 0.03
+    medium = td.Medium(permittivity=4)
+    box = td.Structure(geometry=td.Box(size=(box_size, td.inf, td.inf)), medium=medium)
+    src = td.UniformCurrentSource(
+        source_time=td.GaussianPulse(freq0=2e14, fwidth=1e13),
+        size=(0, 0, 0),
+        polarization="Ex",
+    )
+    sim = td.Simulation(
+        size=(10, 10, 10),
+        sources=[src],
+        structures=[box],
+        run_time=1e-12,
+        grid_spec=td.GridSpec.auto(wavelength=1),
+    )
+
+    # Warning raised as structure is too thin
+    assert_log_level(log_capture, "WARNING")
+
+    # Warning not raised if structure is higher index
+    log_capture.clear()
+    box2 = box.updated_copy(medium=td.Medium(permittivity=300))
+    sim2 = sim.updated_copy(structures=[box2])
+    assert len(log_capture) == 0
+
+    # Warning not raised if structure is covered by an override structure
+    log_capture.clear()
+    override = td.MeshOverrideStructure(geometry=box.geometry, dl=(box_size, td.inf, td.inf))
+    sim3 = sim.updated_copy(grid_spec=sim.grid_spec.updated_copy(override_structures=[override]))
+    assert len(log_capture) == 0
+    # Also check that the structure boundaries are in the grid
+    ind_mid_cell = int(sim3.grid.num_cells[0] // 2)
+    bounds = [-box_size / 2, box_size / 2]
+    assert np.allclose(bounds, sim3.grid.boundaries.x[ind_mid_cell : ind_mid_cell + 2])
+
+    # Test that the error coming from two thin slabs on top of each other is resolved
+    log_capture.clear()
+    box3 = td.Structure(
+        geometry=td.Box(center=(box_size, 0, 0), size=(box_size, td.inf, td.inf)), medium=medium
+    )
+    sim4 = sim.updated_copy(structures=[box3, box])
+    assert_log_level(log_capture, "WARNING")
+
+
 def test_shapely_strtree_warnings():
 
     with warnings.catch_warnings():

tests/test_components/test_simulation.py

@@ -769,7 +769,7 @@ def test_large_grid_size(log_capture, grid_size, log_level):
     assert_log_level(log_capture, log_level)
 
 
-@pytest.mark.parametrize("box_size,log_level", [(0.001, "INFO"), (9.9, "WARNING"), (20, "INFO")])
+@pytest.mark.parametrize("box_size,log_level", [(0.1, "INFO"), (9.9, "WARNING"), (20, "INFO")])
 def test_sim_structure_gap(log_capture, box_size, log_level):
     """Make sure the gap between a structure and PML is not too small compared to lambda0."""
     medium = td.Medium(permittivity=2)

tidy3d/components/grid/mesher.py

@@ -19,9 +19,15 @@
 from ..medium import AnisotropicMedium, Medium2D, PECMedium
 from ...exceptions import SetupError, ValidationError
 from ...constants import C_0, fp_eps
+from ...log import log
 
 _ROOTS_TOL = 1e-10
 
+# Shrink min_step a little so that if e.g. a structure has target dl = 0.1 and a width of 0.1,
+# a grid point will be added on both sides of the structure. Without this factor, the mesher
+# ``is_close`` check will deem that the second point is too close.
+MIN_STEP_SCALE = 0.9999
+
 
 class Mesher(Tidy3dBaseModel, ABC):
     """Abstract class for automatic meshing."""
@@ -122,7 +128,7 @@ def parse_structures(
             structures_ordered, wavelength, min_steps_per_wvl, dl_min, axis
         )
         # Smallest of the maximum steps
-        min_step = np.amin(structure_steps)
+        min_step = MIN_STEP_SCALE * np.amin(structure_steps)
 
         # If empty simulation, return
         if len(structures) == 1:
@@ -135,40 +141,74 @@ def parse_structures(
         tree = self.bounds_2d_tree(struct_bbox)
 
         intervals = {"coords": list(domain_bounds), "structs": [[]]}
-        # Iterate in reverse order as latter structures override earlier ones. To properly handle
-        # containment then we need to populate interval coordinates starting from the top.
-        # If a structure is found to be completely contained, the corresponding ``struct_bbox`` is
-        # set to ``None``.
-        for str_ind in range(len(structures_ordered) - 1, -1, -1):
-            # 3D and 2D bounding box of current structure
-            bbox = struct_bbox[str_ind]
-            if bbox is None:
-                # Structure has been removed because it is completely contained
-                continue
-            bbox_2d = shapely_box(bbox[0, 0], bbox[0, 1], bbox[1, 0], bbox[1, 1])
-
-            # List of structure indexes that may intersect the current structure in 2D
-            try:
-                query_inds = tree.query_items(bbox_2d)
-            except AttributeError:
-                query_inds = tree.query(bbox_2d)
-
-            # Remove all lower structures that the current structure completely contains
-            inds_lower = [
-                ind for ind in query_inds if ind < str_ind and struct_bbox[ind] is not None
-            ]
-            query_bbox = [struct_bbox[ind] for ind in inds_lower]
-            bbox_contains_inds = self.contains_3d(bbox, query_bbox)
-            for ind in bbox_contains_inds:
-                struct_bbox[inds_lower[ind]] = None
-
-            # List of structure bboxes that contain the current structure in 2D
-            inds_upper = [ind for ind in query_inds if ind > str_ind]
-            query_bbox = [struct_bbox[ind] for ind in inds_upper if struct_bbox[ind] is not None]
-            bbox_contained_2d = self.contained_2d(bbox, query_bbox)
-
-            # Handle insertion of the current structure bounds in the intervals
-            intervals = self.insert_bbox(intervals, str_ind, bbox, bbox_contained_2d, min_step)
+        """ Build the ``intervals`` dictionary. ``intervals["coords"]`` gets populated based on the
+        bounding boxes of all structures in the list (some filtering is done to exclude points that
+        will be too close together compared to the absolute lower desired step). At every point, the
+        ``"structs"`` list has length one lower than the ``"coords"`` list, and each element is
+        another list of all structure indexes that are found in the interval formed by ``coords[i]``
+        and ``coords[i + 1]``. This only includes structures that have a physical presence in the
+        interval, i.e. it excludes structures that are completely covered by higher-up ones.
+
+        To build this, we iterate in reverse order as latter structures override earlier ones.
+        We also handle containment in the following way (note - we work with bounding boxes only):
+         1. If a structure that is lower in the list than the current structure is found to be
+          completely contained in 3D, the corresponding ``struct_bbox`` is immediately set to
+          ``None`` and nothing more will be done using that structure.
+         2. If the current structure is covering an interval but there's a higher-up structure that
+          contains it in 2D and also covers the same interval, then it will not be added to the
+          ``intervals["structs"] list for that interval.
+         3. If the current structure is found to not cover any interval, its bounding box
+          is set to ``None``, so that it will not affect structures that lie below it as per point
+          2. A warning is also raised since the structure will have an unpredictable effect on the
+          material coefficients used in the simulation.
+        """
+
+        with log:
+            for str_ind in range(len(structures_ordered) - 1, -1, -1):
+                # 3D and 2D bounding box of current structure
+                bbox = struct_bbox[str_ind]
+                if bbox is None:
+                    # Structure has been removed because it is completely contained
+                    continue
+                bbox_2d = shapely_box(bbox[0, 0], bbox[0, 1], bbox[1, 0], bbox[1, 1])
+
+                # List of structure indexes that may intersect the current structure in 2D
+                try:
+                    query_inds = tree.query_items(bbox_2d)
+                except AttributeError:
+                    query_inds = tree.query(bbox_2d)
+
+                # Remove all lower structures that the current structure completely contains
+                inds_lower = [
+                    ind for ind in query_inds if ind < str_ind and struct_bbox[ind] is not None
+                ]
+                query_bbox = [struct_bbox[ind] for ind in inds_lower]
+                bbox_contains_inds = self.contains_3d(bbox, query_bbox)
+                for ind in bbox_contains_inds:
+                    struct_bbox[inds_lower[ind]] = None
+
+                # List of structure bboxes that contain the current structure in 2D
+                inds_upper = [ind for ind in query_inds if ind > str_ind]
+                query_bbox = [
+                    struct_bbox[ind] for ind in inds_upper if struct_bbox[ind] is not None
+                ]
+                bbox_contained_2d = self.contained_2d(bbox, query_bbox)
+
+                # Handle insertion of the current structure bounds in the intervals
+                # The intervals list is modified in-place
+                too_small = self.insert_bbox(intervals, str_ind, bbox, bbox_contained_2d, min_step)
+                if too_small and (bbox[1, 2] - bbox[0, 2]) > 0:
+                    # If the structure is too small (but not 0D), issue a warning
+                    log.warning(
+                        f"A structure has a nonzero dimension along axis {'xyz'[axis]}, which "
+                        "is however too small compared to the generated mesh step along that "
+                        "direction. This could produce unpredictable results. We recommend "
+                        "increasing the resolution, or adding a mesh override structure to ensure "
+                        "that all geometries are at least one pixel thick along all dimensions."
+                    )
+                    # Also remove this structure from the bbox list so that lower structures
+                    # will not be affected by it, as it was not added to any interval.
+                    struct_bbox[str_ind] = None
 
         # Truncate intervals to domain bounds
         coords = np.array(intervals["coords"])
@@ -226,17 +266,29 @@ def insert_bbox(
             List of 3D bounding boxes that contain the current structure in 2D.
         min_step : float
             Absolute minimum interval size to impose.
+
+        Returns
+        -------
+        structure_too_small : bool
+            True if the structure did not span any interval after coordinates were inserted. This
+            would happen if the structure size is too small compared to the minimum step.
+
+        Note
+        ----
+            This function modifies ``intervals`` in-place.
         """
 
         coords = intervals["coords"]
         structs = intervals["structs"]
 
+        min_step_check = MIN_STEP_SCALE * min_step
+
         # Left structure bound
         bound_coord = str_bbox[0, 2]
         indsmin = np.nonzero(bound_coord <= coords)[0]
-        indmin = int(indsmin[0])  # coordinate is in interval index ``indmin - 1````
-        is_close_l = self.is_close(bound_coord, coords, indmin - 1, min_step)
-        is_close_r = self.is_close(bound_coord, coords, indmin, min_step)
+        indmin = int(indsmin[0])  # coordinate is in interval index ``indmin - 1``
+        is_close_l = self.is_close(bound_coord, coords, indmin - 1, min_step_check)
+        is_close_r = self.is_close(bound_coord, coords, indmin, min_step_check)
         is_contained = self.is_contained(bound_coord, bbox_contained_2d)
 
         # Decide on whether coordinate should be inserted or indmin modified
@@ -254,8 +306,8 @@ def insert_bbox(
         bound_coord = str_bbox[1, 2]
         indsmax = np.nonzero(bound_coord >= coords)[0]
         indmax = int(indsmax[-1])  # coordinate is in interval index ``indmax``
-        is_close_l = self.is_close(bound_coord, coords, indmax, min_step)
-        is_close_r = self.is_close(bound_coord, coords, indmax + 1, min_step)
+        is_close_l = self.is_close(bound_coord, coords, indmax, min_step_check)
+        is_close_r = self.is_close(bound_coord, coords, indmax + 1, min_step_check)
         is_contained = self.is_contained(bound_coord, bbox_contained_2d)
 
         # Decide on whether coordinate should be inserted or indmax modified
@@ -278,7 +330,7 @@ def insert_bbox(
             if not self.is_contained(mid_coord, bbox_contained_2d):
                 structs[interval_ind].append(str_ind)
 
-        return {"coords": coords, "structs": structs}
+        return indmin >= indmax
 
     @staticmethod
     def reorder_structures_enforced_to_end(
@@ -501,7 +553,7 @@ def filter_min_step(
         coords_filter = [interval_coords[0]]
         steps_filter = []
         for coord_ind, coord in enumerate(interval_coords[1:]):
-            if coord - coords_filter[-1] > min_step:
+            if coord - coords_filter[-1] >= min_step:
                 coords_filter.append(coord)
                 steps_filter.append(max_steps[coord_ind])