- 
                Notifications
    You must be signed in to change notification settings 
- Fork 65
feat: surface monitors (FXC-1635) #2891
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: develop
Are you sure you want to change the base?
feat: surface monitors (FXC-1635) #2891
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Additional Comments (3)
- 
tidy3d/components/tcad/data/monitor_data/abstract.py, line 63 (link)logic: Method _symmetry_expanded_copy_basewas removed but is still being called here. This will cause an AttributeError at runtime.
- 
tidy3d/components/data/unstructured/base.py, line 235-237 (link)logic: Return type should be UnstructuredDatasetto match the new class hierarchy
- 
tidy3d/components/data/unstructured/base.py, line 641 (link)syntax: Return type annotation should be UnstructuredDatasetto match the method implementation
28 files reviewed, 16 comments
| Diff CoverageDiff: origin/develop...HEAD, staged and unstaged changes
 Summary
 tidy3d/components/data/dataset.pyLines 516-524   516     @property
  517     def intensity(self) -> TriangularSurfaceDataset:
  518         """Return the sum of the squared absolute electric field components."""
  519         if self.E is None:
! 520             raise ValueError(
  521                 "Could not calculate intensity: the dataset does not contain E field information."
  522             )
  523         intensity = self.E.norm(dim="axis") ** 2
  524         return intensityLines 553-561   553 
  554     @property
  555     def grid_locations(self) -> dict[str, str]:
  556         """Maps field components to the string key of their grid locations on the yee lattice."""
! 557         raise RuntimeError("Function 'grid_location' does not apply to surface monitors.")
  558 
  559     @property
  560     def symmetry_eigenvalues(self) -> dict[str, Callable[[Axis], float]]:
  561         """Maps field components to their (positive) symmetry eigenvalues."""tidy3d/components/data/monitor_data.pyLines 1565-1576   1565         :class:`ElectromagneticSurfaceFieldData`
  1566             A data object with the symmetry expanded fields.
  1567         """
  1568 
! 1569         if all(sym == 0 for sym in self.symmetry):
! 1570             return self
  1571 
! 1572         return self._updated(self._symmetry_update_dict)
  1573 
  1574     @property
  1575     def symmetry_expanded_copy(self) -> AbstractFieldData:
  1576         """Create a copy of the :class:`.ElectromagneticSurfaceFieldData` with fields expanded based on symmetry.tidy3d/components/data/unstructured/base.pyLines 375-383   375 
  376         # Defer to the implementation of the ufunc on unwrapped values.
  377         inputs = tuple(x.values if isinstance(x, type(self)) else x for x in inputs)
  378         if out:
! 379             kwargs["out"] = tuple(x.values if isinstance(x, type(self)) else x for x in out)
  380         result = getattr(ufunc, method)(*inputs, **kwargs)
  381 
  382         if type(result) is tuple:
  383             # multiple return valuesLines 921-929   921             boundary_edges_vtk = feature_edges.GetOutput()
  922 
  923         if boundary_edges_vtk.GetNumberOfCells() == 0:
  924             # Mesh is watertight, no boundary
! 925             return np.array([], dtype=int)
  926 
  927         # Get the original point indices
  928         original_ids_array = vtk["vtk_to_numpy"](
  929             boundary_edges_vtk.GetPointData().GetArray("vtkOriginalPointIds")Lines 964-972   964         if isinstance(symmetry, XrDataArray):
  965             value_dims = set(self.values.dims) - {"index"}
  966             sym_dims = set(symmetry.dims)
  967             if not sym_dims.issubset(value_dims):
! 968                 raise DataError(
  969                     f"Symmetry xarray dimensions {sym_dims} must be a subset of values dimensions {value_dims}"
  970                 )
  971             # Check that coordinates match along shared dimensions
  972             for dim in sym_dims:Lines 970-978   970                 )
  971             # Check that coordinates match along shared dimensions
  972             for dim in sym_dims:
  973                 if not np.array_equal(symmetry.coords[dim], self.values.coords[dim]):
! 974                     raise DataError(
  975                         f"Coordinate values for dimension '{dim}' must match between symmetry and values"
  976                     )
  977 
  978         if reflection_only:Lines 1176-1184   1176         spatial_dims_given = any(comp is not None for comp in [x, y, z])
  1177 
  1178         if spatial_dims_given:
  1179             if any(comp is None for comp in [x, y, z]):
! 1180                 raise DataError("Must provide either all or none of 'x', 'y', and 'z'")
  1181 
  1182             # For spatial interpolation, use legacy default of 0 if not specified
  1183             if fill_value is None:
  1184                 log.warning(Lines 1242-1250   1242         interp_kwargs = {"method": method}
  1243         if fill_value != "extrapolate":
  1244             interp_kwargs["kwargs"] = {"fill_value": fill_value}
  1245         else:
! 1246             interp_kwargs["kwargs"] = {"fill_value": "extrapolate"}
  1247 
  1248         return self.updated_copy(
  1249             values=self.values.interp(**coords_kwargs_only_lists, **interp_kwargs)
  1250         )tidy3d/components/monitor.pyLines 1649-1672   1649         In general, this is severely overestimated for surface monitors.
  1650         """
  1651 
  1652         # estimation based on triangulated surface when it crosses cells in xy plane
! 1653         num_tris = num_cells * 6
! 1654         num_points = num_cells * 4
  1655 
  1656         # storing 3 coordinate components per point
! 1657         storage = 3 * BYTES_REAL * num_points
  1658 
  1659         # storing 3 indices per triangle
! 1660         storage += 3 * BYTES_REAL * num_tris
  1661 
  1662         # EH field values + normal field
! 1663         storage += (
  1664             BYTES_COMPLEX * num_points * len(self.freqs) * len(self.fields) * 3
  1665             + 3 * num_points * BYTES_REAL
  1666         )
  1667 
! 1668         return storage
  1669 
  1670     def _storage_size_solver(self, num_cells: int, tmesh: ArrayFloat1D) -> int:
  1671         """Size of intermediate data recorded by the monitor during a solver run."""Lines 1670-1686   1670     def _storage_size_solver(self, num_cells: int, tmesh: ArrayFloat1D) -> int:
  1671         """Size of intermediate data recorded by the monitor during a solver run."""
  1672 
  1673         # fields
! 1674         storage = BYTES_COMPLEX * num_cells * len(self.freqs) * len(self.fields) * 3
  1675 
  1676         # fields valid map
! 1677         storage += BYTES_REAL * num_cells * len(self.freqs) * len(self.fields) * 3
  1678 
  1679         # auxiliary variables (normals and locations)
! 1680         storage += BYTES_REAL * num_cells * 7 * 4
  1681 
! 1682         return storage
  1683 
  1684 
  1685 class SurfaceFieldTimeMonitor(AbstractSurfaceMonitor, TimeMonitor):
  1686     """:class:`Monitor` that records electromagnetic fields in the time domain on PEC surfaces.Lines 1716-1758   1716     def storage_size(self, num_cells: int, tmesh: ArrayFloat1D) -> int:
  1717         """Size of monitor storage given the number of points after discretization.
  1718         In general, this is severely overestimated for surface monitors.
  1719         """
! 1720         num_steps = self.num_steps(tmesh)
  1721 
  1722         # estimation based on triangulated surface when it crosses cells in xy plane
! 1723         num_tris = num_cells * 6
! 1724         num_points = num_cells * 4
  1725 
  1726         # storing 3 coordinate components per point
! 1727         storage = 3 * BYTES_REAL * num_points
  1728 
  1729         # storing 3 indices per triangle
! 1730         storage += 3 * BYTES_REAL * num_tris
  1731 
  1732         # EH field values + normal field
! 1733         storage += (
  1734             BYTES_COMPLEX * num_points * num_steps * len(self.fields) * 3
  1735             + 3 * num_points * BYTES_REAL
  1736         )
  1737 
! 1738         return storage
  1739 
  1740     def _storage_size_solver(self, num_cells: int, tmesh: ArrayFloat1D) -> int:
  1741         """Size of intermediate data recorded by the monitor during a solver run."""
  1742 
! 1743         num_steps = self.num_steps(tmesh)
  1744 
  1745         # fields
! 1746         storage = BYTES_COMPLEX * num_cells * num_steps * len(self.fields) * 3
  1747 
  1748         # fields valid map
! 1749         storage += BYTES_REAL * num_cells * num_steps * len(self.fields) * 3
  1750 
  1751         # auxiliary variables (normals and locations)
! 1752         storage += BYTES_REAL * num_cells * 7 * 4
  1753 
! 1754         return storage
  1755 
  1756 
  1757 SurfaceMonitorType = Union[
  1758     SurfaceFieldMonitor,tidy3d/components/simulation.pyLines 3967-3988   3967 
  3968     def _get_surface_monitor_bounds_self(self, monitor: SurfaceMonitorType) -> list[Bound]:
  3969         """Intersect a surface monitor with the bounding box of each PEC structure."""
  3970 
! 3971         sim_box = Box(center=self.center, size=self.size)
! 3972         mnt_bounds = Box.bounds_intersection(monitor.bounds, sim_box.bounds)
  3973 
! 3974         if _is_pec_like(self.medium):
! 3975             return [mnt_bounds]
  3976 
! 3977         bounds = []
! 3978         for structure in self.structures:
! 3979             if _is_pec_like(structure.medium):
! 3980                 intersection_bounds = Box.bounds_intersection(mnt_bounds, structure.geometry.bounds)
! 3981                 if all(bmin <= bmax for bmin, bmax in zip(*intersection_bounds)):
! 3982                     bounds.append(intersection_bounds)
  3983 
! 3984         return bounds
  3985 
  3986     @pydantic.validator("monitors", always=True)
  3987     @skip_if_fields_missing(["medium", "structures", "size"])
  3988     def error_empty_surface_monitor(cls, val, values):tidy3d/components/viz/axes_utils.pyLines 151-161   151 
  152         ipython = get_ipython()
  153         if ipython is None:
  154             return False
! 155         return "IPKernelApp" in ipython.config
! 156     except (ImportError, AttributeError):
! 157         return False
  158 
  159 
  160 def add_plotter_if_none(plot):
  161     """Decorates ``plot(*args, **kwargs, plotter=None)`` function for PyVista.Lines 206-214   206         plotter = plot(*args, **kwargs)
  207 
  208         # Show if we created the plotter and show=True
  209         if plotter_created and show:
! 210             return plotter.show()
  211 
  212         return plotter
  213 
  214     return _plottidy3d/packaging.pyLines 171-180   171                 import pyvista as pv
  172 
  173                 pyvista["mod"] = pv
  174 
! 175             except ImportError as exc:
! 176                 raise Tidy3dImportError(
  177                     "The package 'pyvista' is required for this operation, but it was not found. "
  178                     "Please install the 'pyvista' dependencies using, for example, "
  179                     "'pip install pyvista' or 'pip install tidy3d[pyvista]'."
  180                 ) from exc | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks @dbochkov-flexcompute! Great refactor of the unstructured grid and the data array tests!
Left some minor comments but overall looks good to me
        
          
                tests/test_data/test_data_arrays.py
              
                Outdated
          
        
      | # SURFACE_FIELD_MONITOR, | ||
| # SURFACE_FIELD_TIME_MONITOR, | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We can remove this, right?
        
          
                tests/test_data/test_monitor_data.py
              
                Outdated
          
        
      | assert np.all(np.isreal(intensity.values.values)) | ||
|  | ||
|  | ||
| def test_surface_field_time_data(): | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This test and the one above are very similar. I wonder if they could be combined using @pytest.mark.parametrize().
I guess some of the duplicated asserts would be averted though maybe is less clear/explicit? What do you think?
The following two I think could benefit from combining them and have 2 instead of 4 tests?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Big effort! Looks great overall just a couple of comments.
| grid_spec=td.GridSpec.auto(wavelength=1), | ||
| ) | ||
|  | ||
| # monitor must be volumetric | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Wondering about 2d simulations - does the surface monitor work for those?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
currently, not, which reminds me to add validator for that. I should be able to add this option later either in this PR or in a follow-up PR. To keep it simple, I would still probably return a 3d surface (ribbon) with a finite width in the zero size dimension instead of a 2d curve
| @classmethod | ||
| def _point_dims(cls) -> pd.PositiveInt: | ||
| """Dimensionality of stored surface point coordinates.""" | ||
| return 3 | ||
|  | ||
| @classmethod | ||
| def _cell_num_vertices(cls) -> pd.PositiveInt: | ||
| """Number of vertices in a cell.""" | ||
| return 3 | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Wondering whether these should rather be properties?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
we need access to these in validators and in class method from_vtu, and it seems like @classmethod + @Property is depreciated in python 11 python/cpython#119628
        
          
                tidy3d/components/viz/axes_utils.py
              
                Outdated
          
        
      | def _plot(*args, **kwargs): | ||
| """New plot function using a generated plotter if None.""" | ||
| # Extract decorator-specific parameters | ||
| plotter = kwargs.get("plotter") | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
plotter is not required to be a keyword argument as per the definition of plot(), so if you were to pass a plotter positionally, then the decorator here would still inject a keyword copy and you'd end up with an exception. This would be more robustly handled via inspect.signature().
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed
| sel_kwargs_only_lists = { | ||
| key: value if isinstance(value, list) else [value] for key, value in sel_kwargs.items() | ||
| } | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since this converts all kwargs to lists, wouldn't this end up breaking "non-indexing" arguments? .sel has a couple of those: https://docs.xarray.dev/en/latest/generated/xarray.DataArray.sel.html#xarray.DataArray.sel
Or do we not care since this is a private method and we know we're not going to use those arguments?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed
| if all(sym == 0 for sym in self.symmetry): | ||
| return data | ||
|  | ||
| new_data = copy.copy(data) | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess this is safe but it doesn't look to me like the copy is necessary here? new_data is never mutated.
| def _is_notebook() -> bool: | ||
| """Detect if running in Jupyter notebook. | ||
| Returns | ||
| ------- | ||
| bool | ||
| True if running in Jupyter notebook, False otherwise. | ||
| """ | ||
| try: | ||
| from IPython import get_ipython | ||
|  | ||
| ipython = get_ipython() | ||
| if ipython is None: | ||
| return False | ||
| return "IPKernelApp" in ipython.config | ||
| except (ImportError, AttributeError): | ||
| return False | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I haven't played around with this much but I'm wondering whether the approach is robust against different notebook environments like JupyterLab, VSCode, Colab? It might be more robust to check the class name of the shell directly, i.e., based on get_ipython().__class__.__name__. Would need to test though.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
expanded this function using recommendations from the web and tested in jupyter-lab, vscode, and colab - all seem to work alright
| colocate: Literal[False] = pydantic.Field( | ||
| False, | ||
| title="Colocate Fields", | ||
| description="For surface monitors fields are always colocated on surface.", | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
if they are always colocated on the surface, would colocate always be True instead of False?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yeah, it's a bit confusing here for some internal reasons. But I will change it to True
        
          
                tidy3d/components/medium.py
              
                Outdated
          
        
      | return Lorentz.from_nk(n, k, freq, **kwargs) | ||
|  | ||
|  | ||
| def _is_pec_like(medium: MediumType3D) -> bool: | 
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is there a reason to have this be a separate function instead of a property of the mediums? I.e. - medium.is_pec_like?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
no reason, implemented now
a few examples (unfortunately, need to re-ran to see visualization):
Greptile Overview
Updated On: 2025-10-14 22:39:42 UTC
Greptile Summary
This PR introduces comprehensive surface monitoring capabilities to Tidy3D, allowing users to monitor electromagnetic fields directly on Perfect Electric Conductor (PEC) and lossy metal surfaces. The implementation adds two new monitor types:
SurfaceFieldMonitorfor frequency-domain monitoring andSurfaceFieldTimeMonitorfor time-domain monitoring. These monitors capture electric and magnetic fields on triangular surface meshes rather than regular Cartesian grids.The core functionality includes new type definitions (
EMSurfaceField), unstructured dataset support (TriangularSurfaceDataset), and corresponding data classes (SurfaceFieldData,SurfaceFieldTimeData) that handle field storage and manipulation on surface meshes. The monitors integrate with the existing simulation validation framework to ensure they only operate where PEC or lossy metal structures exist. PyVista support has been added as an optional dependency for 3D surface visualization, following the established pattern of optional visualization libraries.The implementation extends the existing monitor infrastructure while maintaining backward compatibility, adding specialized handling for surface current density calculations through cross products of magnetic field discontinuities across conductor boundaries. All new functionality includes comprehensive test coverage and follows the established code patterns for monitor validation, data handling, and visualization.
Important Files Changed
Changed Files
Confidence score: 3/5
AbstractTcadFieldData.symmetry_expanded_copy(line 63) calling removed_symmetry_expanded_copy_basemethod, type annotation inconsistencies in base classes, and complex class hierarchy changes that could affect existing functionalitytidy3d/components/tcad/data/monitor_data/abstract.pywhich has a critical broken method call, and the unstructured data base class refactoring which may impact inheritance behaviorContext used:
dashboard- When modifying a piece of logic, ensure the change is propagated to all independent functions or cal... (source)