diff --git a/CHANGELOG.md b/CHANGELOG.md index 9ca36109fa..fd95d6723e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added support for `tidy3d-extras`, an optional plugin that enables more accurate local mode solving via subpixel averaging. +- Added support for `symlog` and `log` scale plotting in `Scene.plot_eps()` and `Scene.plot_structures_property()` methods. The `symlog` scale provides linear behavior near zero and logarithmic behavior elsewhere, while 'log' is a base 10 logarithmic scale. ### Changed - Improved performance of antenna metrics calculation by utilizing cached wave amplitude calculations instead of recomputing wave amplitudes for each port excitation in the `TerminalComponentModelerData`. diff --git a/tests/test_components/test_scene.py b/tests/test_components/test_scene.py index 81958b8426..6f7e5741c3 100644 --- a/tests/test_components/test_scene.py +++ b/tests/test_components/test_scene.py @@ -9,6 +9,7 @@ import tidy3d as td from tidy3d.components.scene import MAX_GEOMETRY_COUNT, MAX_NUM_MEDIUMS +from tidy3d.exceptions import SetupError from ..utils import SIM_FULL, cartesian_to_unstructured @@ -370,7 +371,7 @@ def test_plot_property(): name="Si_MultiPhysics", ) - def try_plotting(mpm, display=False): + def try_plotting(mpm, display=False, scale=None): # Structure struct = td.Structure( geometry=td.Box(size=(2, 2, 2), center=(0, 0, 0)), @@ -384,10 +385,16 @@ def try_plotting(mpm, display=False): ) _, ax = plt.subplots(1, 4, figsize=(20, 4)) - scene.plot_structures_property(z=0, property="N_a", ax=ax[0]) - scene.plot_structures_property(z=0, property="N_d", ax=ax[1]) - scene.plot_structures_property(z=0, property="doping", ax=ax[2]) - scene.plot_structures_property(z=0, ax=ax[3]) # eps + if scale: + scene.plot_structures_property(z=0, property="N_a", ax=ax[0], scale=scale) + scene.plot_structures_property(z=0, property="N_d", ax=ax[1], scale=scale) + scene.plot_structures_property(z=0, property="doping", ax=ax[2], scale=scale) + scene.plot_structures_property(z=0, ax=ax[3], scale=scale) # eps + else: + scene.plot_structures_property(z=0, property="N_a", ax=ax[0]) + scene.plot_structures_property(z=0, property="N_d", ax=ax[1]) + scene.plot_structures_property(z=0, property="doping", ax=ax[2]) + scene.plot_structures_property(z=0, ax=ax[3]) # eps if display: plt.show() @@ -395,6 +402,7 @@ def try_plotting(mpm, display=False): const_doping = td.ConstantDoping(concentration=1e15) mpm = mpm.updated_copy(charge=semicon.updated_copy(N_a=[const_doping], N_d=[const_doping])) try_plotting(mpm, display=display_plots) + try_plotting(mpm, display=display_plots, scale="symlog") # add some Gaussian doping gaussian_box = td.GaussianDoping( @@ -402,6 +410,7 @@ def try_plotting(mpm, display=False): ) mpm = mpm.updated_copy(charge=semicon.updated_copy(N_a=[gaussian_box], N_d=[gaussian_box])) try_plotting(mpm, display=display_plots) + try_plotting(mpm, display=display_plots, scale="symlog") # now try with a custom doping x = np.linspace(-1, 1, 30) @@ -419,3 +428,60 @@ def try_plotting(mpm, display=False): custom_box1 = td.CustomDoping(center=(0, 0, 0), size=(2, 2, 2), concentration=concentration) mpm = mpm.updated_copy(charge=semicon.updated_copy(N_a=[custom_box1], N_d=[custom_box1])) try_plotting(mpm, display=display_plots) + try_plotting(mpm, display=display_plots, scale="symlog") + + +def test_log_scale_with_custom_limits(): + """Test log scale with custom limits.""" + + # Create a scene with different permittivity values + scene = td.Scene( + structures=[ + td.Structure( + geometry=td.Box(size=(1, 1, 1), center=(-1, 0, 0)), + medium=td.MultiPhysicsMedium( + optical=td.Medium(permittivity=1.0), + ), + ), + td.Structure( + geometry=td.Box(size=(1, 1, 1), center=(0, 0, 0)), + medium=td.MultiPhysicsMedium( + optical=td.Medium(permittivity=10.0), + ), + ), + td.Structure( + geometry=td.Box(size=(1, 1, 1), center=(1, 0, 0)), + medium=td.MultiPhysicsMedium( + optical=td.Medium(permittivity=100.0), + ), + ), + ], + medium=td.MultiPhysicsMedium( + optical=td.Medium(permittivity=1.0), + ), + ) + + # Test log scale with custom limits + _ = scene.plot_eps(x=0, scale="log", eps_lim=(1e-2, 100)) + plt.close() + + _ = scene.plot_eps(x=0, scale="log", eps_lim=(1e-5, 100)) + plt.close() + + with pytest.raises(SetupError, match="Log scale cannot be used with non-positive values."): + _ = scene.plot_eps(x=0, scale="log", eps_lim=(-1e-2, 100)) + plt.close() + + _ = scene.plot_structures_property(x=0, property="eps", scale="log", limits=(1e-2, 100)) + plt.close() + + with pytest.raises(SetupError, match="Log scale cannot be used with non-positive values."): + _ = scene.plot_structures_property(x=0, property="eps", scale="log", limits=(-2e-2, 100)) + plt.close() + + # Test that invalid scale raises error + with pytest.raises( + SetupError, match="The scale 'invalid' is not supported for plotting structures property." + ): + _ = scene.plot_structures_property(x=0, property="eps", scale="invalid") + plt.close() diff --git a/tests/test_data/test_sim_data.py b/tests/test_data/test_sim_data.py index 8d258f611e..dcbe5ab560 100644 --- a/tests/test_data/test_sim_data.py +++ b/tests/test_data/test_sim_data.py @@ -13,7 +13,7 @@ from tidy3d.components.data.sim_data import SimulationData from tidy3d.components.file_util import replace_values from tidy3d.components.monitor import FieldMonitor, FieldTimeMonitor, ModeMonitor -from tidy3d.exceptions import DataError, Tidy3dKeyError +from tidy3d.exceptions import DataError, SetupError, Tidy3dKeyError from ..utils import get_nested_shape from .test_data_arrays import FIELD_MONITOR, SIM, SIM_SYM @@ -527,3 +527,19 @@ def test_to_mat_file(tmp_path): sim_data = make_sim_data() path = str(tmp_path / "test.mat") sim_data.to_mat_file(path) + + +def test_plot_field_monitor_data_unsupported_scale(): + """Test plot_field_monitor_data with unsupported scale to trigger SetupError.""" + sim_data = make_sim_data() + + # Test with unsupported scale + with pytest.raises( + SetupError, match="The scale 'invalid' is not supported for plotting field data" + ): + sim_data.plot_field_monitor_data( + field_monitor_data=sim_data.monitor_data["field"], + field_name="Ex", + val="real", + scale="invalid", + ) diff --git a/tidy3d/components/data/sim_data.py b/tidy3d/components/data/sim_data.py index 639f80c7b3..e4b4b3b13a 100644 --- a/tidy3d/components/data/sim_data.py +++ b/tidy3d/components/data/sim_data.py @@ -27,7 +27,7 @@ from tidy3d.components.types import Ax, Axis, ColormapType, FieldVal, PlotScale, annotate_type from tidy3d.components.viz import add_ax_if_none, equal_aspect from tidy3d.constants import C_0, inf -from tidy3d.exceptions import DataError, FileError, Tidy3dKeyError +from tidy3d.exceptions import DataError, FileError, SetupError, Tidy3dKeyError from tidy3d.log import log from .data_array import FreqDataArray, TimeDataArray @@ -540,7 +540,7 @@ def plot_field_monitor_data( field_data = db_factor * np.log10(np.abs(field_data)) field_data.name += " (dB)" cmap_type = "sequential" - else: + elif scale == "lin": cmap_type = ( "cyclic" if val == "phase" @@ -550,6 +550,8 @@ def plot_field_monitor_data( else "sequential" ) ) + else: + raise SetupError(f"The scale '{scale}' is not supported for plotting field data.") # interp out any monitor.size==0 dimensions monitor = field_monitor_data.monitor diff --git a/tidy3d/components/scene.py b/tidy3d/components/scene.py index 162d40dad9..491b41b08a 100644 --- a/tidy3d/components/scene.py +++ b/tidy3d/components/scene.py @@ -59,6 +59,7 @@ InterpMethod, LengthUnit, PermittivityComponent, + PlotScale, PriorityMode, Shapely, Size, @@ -632,11 +633,19 @@ def _get_structure_plot_params( return plot_params @staticmethod - def _add_cbar(vmin: float, vmax: float, label: str, cmap: str, ax: Ax = None) -> None: + def _add_cbar( + vmin: float, + vmax: float, + label: str, + cmap: str, + ax: Ax = None, + norm: mpl.colors.Normalize = None, + ) -> None: """Add a colorbar to plot.""" - norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.15) + if norm is None: + norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) mappable = mpl.cm.ScalarMappable(norm=norm, cmap=cmap) plt.colorbar(mappable, cax=cax, label=label) @@ -804,6 +813,8 @@ def plot_eps( ax: Ax = None, hlim: Optional[tuple[float, float]] = None, vlim: Optional[tuple[float, float]] = None, + eps_lim: tuple[Union[float, None], Union[float, None]] = (None, None), + scale: PlotScale = "lin", ) -> Ax: """Plot each of scene's components on a plane defined by one nonzero x,y,z coordinate. The permittivity is plotted in grayscale based on its value at the specified frequency. @@ -828,6 +839,10 @@ def plot_eps( The x range if plotting on xy or xz planes, y range if plotting on yz plane. vlim : Tuple[float, float] = None The z range if plotting on xz or yz planes, y plane if plotting on xy plane. + eps_lim : Tuple[float, float] = None + Custom limits for eps coloring. + scale : PlotScale = "lin" + Scale for the plot. Either 'lin' for linear, 'log' for log10, 'symlog' for symmetric logarithmic (linear near zero, logarithmic elsewhere), or 'dB' for decibel scale. Returns ------- @@ -838,7 +853,17 @@ def plot_eps( hlim, vlim = Scene._get_plot_lims(bounds=self.bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim) ax = self.plot_structures_eps( - freq=freq, cbar=True, alpha=alpha, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim + freq=freq, + cbar=True, + alpha=alpha, + ax=ax, + x=x, + y=y, + z=z, + hlim=hlim, + vlim=vlim, + eps_lim=eps_lim, + scale=scale, ) ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim) return ax @@ -855,6 +880,7 @@ def plot_structures_eps( cbar: bool = True, reverse: bool = False, eps_lim: tuple[Union[float, None], Union[float, None]] = (None, None), + scale: PlotScale = "lin", ax: Ax = None, hlim: Optional[tuple[float, float]] = None, vlim: Optional[tuple[float, float]] = None, @@ -885,6 +911,8 @@ def plot_structures_eps( Defaults to the structure default alpha. eps_lim : Tuple[float, float] = None Custom limits for eps coloring. + scale : PlotScale = "lin" + Scale for the plot. Either 'lin' for linear, 'log' for log10, 'symlog' for symmetric logarithmic (linear near zero, logarithmic elsewhere), or 'dB' for decibel scale. ax : matplotlib.axes._subplots.Axes = None Matplotlib axes to plot on, if not specified, one is created. hlim : Tuple[float, float] = None @@ -911,6 +939,7 @@ def plot_structures_eps( cbar=cbar, reverse=reverse, limits=eps_lim, + scale=scale, ax=ax, hlim=hlim, vlim=vlim, @@ -931,6 +960,7 @@ def plot_structures_property( cbar: bool = True, reverse: bool = False, limits: tuple[Union[float, None], Union[float, None]] = (None, None), + scale: PlotScale = "lin", ax: Ax = None, hlim: Optional[tuple[float, float]] = None, vlim: Optional[tuple[float, float]] = None, @@ -962,6 +992,9 @@ def plot_structures_property( Defaults to the structure default alpha. limits : Tuple[float, float] = None Custom coloring limits for the property to plot. + scale : PlotScale = "lin" + Scale for the plot. Either 'lin' for linear, 'log' for log10, or 'dB' for decibel scale. + For log scale with negative values, the absolute value is taken before log transformation. ax : matplotlib.axes._subplots.Axes = None Matplotlib axes to plot on, if not specified, one is created. hlim : Tuple[float, float] = None @@ -1018,23 +1051,52 @@ def plot_structures_property( if property_min is None or property_max is None: if property == "eps": eps_min_sim, eps_max_sim = self.eps_bounds(freq=freq, eps_component=eps_component) - if property_min is None: - property_min = eps_min_sim - - if property_max is None: - property_max = eps_max_sim + property_min = property_min if property_min is not None else eps_min_sim + property_max = property_max if property_max is not None else eps_max_sim + linthresh = 1e-2 * property_min if property in ["N_d", "N_a", "doping"]: acceptor_limits, donor_limits = self.doping_bounds() + acceptor_abs_min, donor_abs_min = self.doping_absolute_minimum() if property == "N_d": - property_min = donor_limits[0] - property_max = donor_limits[1] + property_min = property_min if property_min is not None else donor_limits[0] + property_max = property_max if property_max is not None else donor_limits[1] + linthresh = donor_abs_min elif property == "N_a": - property_min = acceptor_limits[0] - property_max = acceptor_limits[1] + property_min = property_min if property_min is not None else acceptor_limits[0] + property_max = property_max if property_max is not None else acceptor_limits[1] + linthresh = acceptor_abs_min elif property == "doping": - property_min = -donor_limits[1] - property_max = acceptor_limits[1] + property_min = property_min if property_min is not None else -donor_limits[1] + property_max = property_max if property_max is not None else acceptor_limits[1] + linthresh = min(acceptor_abs_min, donor_abs_min) + + if np.isclose(linthresh, 0.0): + # fallback to default linthresh of 1e-3 + linthresh = 1e-3 + else: + if np.isclose(property_min, 0.0) or property_min < 0.0: + linthresh = 1e-3 + else: + linthresh = 1e-2 * np.abs(property_min) + + if scale == "lin": + norm = mpl.colors.Normalize(vmin=property_min, vmax=property_max) + elif scale == "log": + # LogNorm doesn't work with negative values, so we need to handle this case + if property_min <= 0 or property_max <= 0: + raise SetupError( + f"Log scale cannot be used with non-positive values. " + f"Property range: [{property_min}, {property_max}]. " + f"Consider using 'symlog' scale instead." + ) + norm = mpl.colors.LogNorm(vmin=property_min, vmax=property_max) + elif scale == "symlog": + norm = mpl.colors.SymLogNorm(linthresh=linthresh, vmin=property_min, vmax=property_max) + else: + raise SetupError( + f"The scale '{scale}' is not supported for plotting structures property." + ) for medium, shape in medium_shapes: if property in ["doping", "N_a", "N_d"]: @@ -1051,7 +1113,17 @@ def plot_structures_property( ) else: self._pcolormesh_shape_doping_box( - x, y, z, alpha, medium, property_min, property_max, shape, ax, property + x, + y, + z, + alpha, + medium, + property_min, + property_max, + shape, + ax, + property, + norm, ) else: # if the background medium is custom medium, it needs to be rendered separately @@ -1069,6 +1141,7 @@ def plot_structures_property( shape=shape, ax=ax, eps_component=eps_component, + norm=norm, ) else: # For custom medium, apply pcolormesh clipped by the shape. @@ -1086,6 +1159,7 @@ def plot_structures_property( ax, grid, eps_component=eps_component, + norm=norm, ) if cbar: @@ -1096,10 +1170,11 @@ def plot_structures_property( label=r"$\rm{Doping} \#/cm^3$", cmap=HEAT_SOURCE_CMAP, ax=ax, + norm=norm, ) else: self._add_cbar_eps( - eps_min=property_min, eps_max=property_max, ax=ax, reverse=reverse + eps_min=property_min, eps_max=property_max, ax=ax, reverse=reverse, norm=norm ) # clean up the axis display @@ -1113,7 +1188,13 @@ def plot_structures_property( return ax @staticmethod - def _add_cbar_eps(eps_min: float, eps_max: float, ax: Ax = None, reverse: bool = False) -> None: + def _add_cbar_eps( + eps_min: float, + eps_max: float, + ax: Ax = None, + reverse: bool = False, + norm: Optional[mpl.colors.Normalize] = None, + ) -> None: """Add a permittivity colorbar to plot.""" Scene._add_cbar( vmin=eps_min, @@ -1121,6 +1202,7 @@ def _add_cbar_eps(eps_min: float, eps_max: float, ax: Ax = None, reverse: bool = label=r"$\epsilon_r$", cmap=STRUCTURE_EPS_CMAP if not reverse else STRUCTURE_EPS_CMAP_R, ax=ax, + norm=norm, ) @staticmethod @@ -1182,6 +1264,7 @@ def _pcolormesh_shape_custom_medium_structure_eps( ax: Ax, grid: Grid, eps_component: Optional[PermittivityComponent] = None, + norm: mpl.colors.Normalize = None, ): """ Plot shape made of custom medium with ``pcolormesh``. @@ -1324,10 +1407,9 @@ def _pcolormesh_shape_custom_medium_structure_eps( eps_shape, clip_path=(polygon_path(shape), ax.transData), cmap=STRUCTURE_EPS_CMAP, - vmin=eps_min, - vmax=eps_max, alpha=alpha, clip_box=ax.bbox, + norm=norm, ) @staticmethod @@ -1339,6 +1421,7 @@ def _get_structure_eps_plot_params( reverse: bool = False, alpha: Optional[float] = None, eps_component: Optional[PermittivityComponent] = None, + norm: Optional[mpl.colors.Normalize] = None, ) -> PlotParams: """Constructs the plot parameters for a given medium in scene.plot_eps().""" @@ -1364,11 +1447,19 @@ def _get_structure_eps_plot_params( plot_params = plot_params.copy(update={"edgecolor": "k", "linewidth": 1}) else: eps_medium = medium._eps_plot(frequency=freq, eps_component=eps_component) - delta_eps = eps_medium - eps_min - delta_eps_max = eps_max - eps_min + 1e-5 - eps_fraction = delta_eps / delta_eps_max - color = eps_fraction if reverse else 1 - eps_fraction - color = min(1, max(color, 0)) # clip in case of custom eps limits + if norm is not None: + # Use the same normalization as the colorbar for consistency + color = norm(eps_medium) + if reverse: + color = 1 - color + color = min(1, max(color, 0)) # clip in case of custom eps limits + else: + # Fallback to linear mapping for backward compatibility + delta_eps = eps_medium - eps_min + delta_eps_max = eps_max - eps_min + 1e-5 + eps_fraction = delta_eps / delta_eps_max + color = eps_fraction if reverse else 1 - eps_fraction + color = min(1, max(color, 0)) # clip in case of custom eps limits plot_params = plot_params.copy(update={"facecolor": str(color)}) return plot_params @@ -1384,6 +1475,7 @@ def _plot_shape_structure_eps( reverse: bool = False, alpha: Optional[float] = None, eps_component: Optional[PermittivityComponent] = None, + norm: Optional[mpl.colors.Normalize] = None, ) -> Ax: """Plot a structure's cross section shape for a given medium, grayscale for permittivity.""" plot_params = self._get_structure_eps_plot_params( @@ -1394,6 +1486,7 @@ def _plot_shape_structure_eps( alpha=alpha, reverse=reverse, eps_component=eps_component, + norm=norm, ) ax = self.box.plot_shape(shape=shape, plot_params=plot_params, ax=ax) return ax @@ -1935,6 +2028,67 @@ def doping_bounds(self): limits[1] = max_value return acceptors_lims, donors_lims + def doping_absolute_minimum(self): + """Get the absolute minimum values of the doping concentrations. + + Returns + ------- + Tuple[float, float] + Absolute minimum values for acceptors and donors respectively. + """ + # Use more reasonable initial values + acceptors_abs_min = np.inf + donors_abs_min = np.inf + + for struct in self.all_structures: + if isinstance(struct.medium.charge, SemiconductorMedium): + electric_spec = struct.medium.charge + + # Process acceptors + acceptors_min = self._get_absolute_minimum_from_doping(electric_spec.N_a) + if acceptors_min < acceptors_abs_min: + acceptors_abs_min = acceptors_min + + # Process donors + donors_min = self._get_absolute_minimum_from_doping(electric_spec.N_d) + if donors_min < donors_abs_min: + donors_abs_min = donors_min + + return acceptors_abs_min, donors_abs_min + + def _get_absolute_minimum_from_doping(self, doping): + """Helper method to get absolute minimum from a single doping specification. + + Parameters + ---------- + doping : Union[float, SpatialDataArray, tuple] + Doping specification to analyze. + + Returns + ------- + float + Absolute minimum value found in the doping specification. + """ + if isinstance(doping, float): + return np.abs(doping) + + # NOTE: This will be deprecated. + if isinstance(doping, SpatialDataArray): + return np.min(np.abs(doping.data.flatten())) + + if isinstance(doping, tuple): + min_values = [] + for doping_box in doping: + if isinstance(doping_box, ConstantDoping): + min_values.append(np.abs(doping_box.concentration)) + elif isinstance(doping_box, GaussianDoping): + min_values.append(np.abs(doping_box.ref_con)) + elif isinstance(doping_box, CustomDoping): + min_values.append(np.min(np.abs(doping_box.concentration.data.flatten()))) + return min(min_values) if min_values else np.inf + + return np.inf + def _pcolormesh_shape_doping_box( self, x: float, @@ -1947,6 +2101,7 @@ def _pcolormesh_shape_doping_box( shape: Shapely, ax: Ax, plt_type: str = "doping", + norm: mpl.colors.Normalize = None, ): """ Plot shape made of structure defined with doping. @@ -2015,10 +2170,9 @@ def _pcolormesh_shape_doping_box( struct_doping_to_plot, clip_path=(polygon_path(shape), ax.transData), cmap=HEAT_SOURCE_CMAP, - vmin=doping_min, - vmax=doping_max, alpha=alpha, clip_box=ax.bbox, + norm=norm, ) def plot_3d(self, width=800, height=800) -> None: diff --git a/tidy3d/components/types/base.py b/tidy3d/components/types/base.py index 56d5c037cd..60c8c8ea03 100644 --- a/tidy3d/components/types/base.py +++ b/tidy3d/components/types/base.py @@ -238,7 +238,7 @@ def __modify_schema__(cls, field_schema): PlotVal = Literal["real", "imag", "abs"] FieldVal = Literal["real", "imag", "abs", "abs^2", "phase"] RealFieldVal = Literal["real", "abs", "abs^2"] -PlotScale = Literal["lin", "dB"] +PlotScale = Literal["lin", "dB", "log", "symlog"] ColormapType = Literal["divergent", "sequential", "cyclic"] """ mode solver """