From 3d99ab96cbc6b54da020d5b99443f69845a18f48 Mon Sep 17 00:00:00 2001 From: tylerflex Date: Mon, 4 Apr 2022 10:54:09 -0700 Subject: [PATCH 1/8] pylint no longer ignore material library and plugins --- .pylintrc | 1 - tidy3d/material_library.py | 5 ++--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/.pylintrc b/.pylintrc index 04bb03c9e1..3a29be50e9 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,5 @@ [MASTER] extension-pkg-allow-list=pydantic -ignore=material_library.py, plugins good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz [BASIC] diff --git a/tidy3d/material_library.py b/tidy3d/material_library.py index 9192ec3691..26df001e56 100644 --- a/tidy3d/material_library.py +++ b/tidy3d/material_library.py @@ -1,4 +1,5 @@ -""" Holds dispersive models for several commonly used optical materials.""" +"""Holds dispersive models for several commonly used optical materials.""" +import json from .components.medium import PoleResidue @@ -6,8 +7,6 @@ def export_matlib_to_file(fname: str = "matlib.json") -> None: """Write the material library to a .json file.""" - import json - mat_lib_dict = {} for mat_name, mat in material_library.items(): mat_lib_dict[mat_name] = {} From 1157aef30f10d4aba90c2a96b790a6a35247c7de Mon Sep 17 00:00:00 2001 From: tylerflex Date: Mon, 4 Apr 2022 11:00:49 -0700 Subject: [PATCH 2/8] cleaned up the plotly and smatrix stuff --- tidy3d/plugins/__init__.py | 2 +- tidy3d/plugins/plotly/__init__.py | 5 +- tidy3d/plugins/plotly/app.py | 60 ++- tidy3d/plugins/plotly/component.py | 19 +- tidy3d/plugins/plotly/data.py | 804 ++++++++++++++++++++++------ tidy3d/plugins/plotly/simulation.py | 78 ++- tidy3d/plugins/plotly/utils.py | 1 + tidy3d/plugins/smatrix/smatrix.py | 11 +- 8 files changed, 758 insertions(+), 222 deletions(-) diff --git a/tidy3d/plugins/__init__.py b/tidy3d/plugins/__init__.py index 626b4923f8..801e227578 100644 --- a/tidy3d/plugins/__init__.py +++ b/tidy3d/plugins/__init__.py @@ -1,4 +1,4 @@ -# import the specific classes / functions needed for the plugins +"""External plugins that have tidy3d as dependency and add functionality.""" from .dispersion.fit import DispersionFitter from .dispersion.fit_web import StableDispersionFitter, AdvancedFitterParam diff --git a/tidy3d/plugins/plotly/__init__.py b/tidy3d/plugins/plotly/__init__.py index 70096d559d..d4a6656b47 100644 --- a/tidy3d/plugins/plotly/__init__.py +++ b/tidy3d/plugins/plotly/__init__.py @@ -1,2 +1,3 @@ -# from .geo import SimulationPlotly, GeometryPlotly, StructurePlotly -# from .data import * +"""Plotly-based app and Simulation plotter.""" +from .app import SimulationDataApp +from .simulation import SimulationPlotly diff --git a/tidy3d/plugins/plotly/app.py b/tidy3d/plugins/plotly/app.py index af490f384d..8b38c3942b 100644 --- a/tidy3d/plugins/plotly/app.py +++ b/tidy3d/plugins/plotly/app.py @@ -1,21 +1,18 @@ -"""Makes an app.""" +"""Makes an app to visualize SimulationData objects.""" from abc import ABC, abstractmethod -from typing import List, Union, Any from typing_extensions import Literal from jupyter_dash import JupyterDash from dash import Dash, dcc, html import pydantic as pd -import sys - -sys.path.append("../../tidy3d") -import tidy3d as td -from tidy3d.components.base import Tidy3dBaseModel from .simulation import SimulationPlotly from .data import DataPlotly +from ...components.base import Tidy3dBaseModel +from ...components.simulation import Simulation +from ...components.data import SimulationData -APP_MODE = Literal["python", "jupyter", "jupyterlab"] +AppMode = Literal["python", "jupyter", "jupyterlab"] DEFAULT_MODE = "jupyterlab" DASH_APP = "Dash App" @@ -23,23 +20,22 @@ class App(Tidy3dBaseModel, ABC): """Basic dash app template: initializes, makes layout, and fires up a server.""" - mode: APP_MODE = pd.Field( + mode: AppMode = pd.Field( DEFAULT_MODE, title="App Mode", - description='Run app differently based on `mode` in `"python"`, `"jupyter"`, `"jupyterlab"`', + description='Run app in mode that is one of `"python"`, `"jupyter"`, `"jupyterlab"`.', ) def _initialize_app(self) -> DASH_APP: """Creates an app based on specs.""" if "jupyter" in self.mode.lower(): return JupyterDash(__name__) - elif "python" in self.mode.lower(): + if "python" in self.mode.lower(): return Dash(__name__) - else: - raise NotImplementedError(f"doesnt support mode={mode}") + raise NotImplementedError(f"App doesn't support mode='{self.mode}'.") @abstractmethod - def _make_layout(self, app: DASH_APP) -> "Dash Layout": + def _make_layout(self, app: DASH_APP) -> dcc.Tabs: """Creates the layout for the app.""" raise NotImplementedError("Must implement in subclass.") @@ -54,7 +50,7 @@ def run(self, debug: bool = False) -> None: app = self._make_app() - if self.mode == "jupyterlab": + if "jupyter" in self.mode.lower(): app.run_server( mode="jupyterlab", port=8090, @@ -66,17 +62,17 @@ def run(self, debug: bool = False) -> None: elif "python" in self.mode.lower(): app.run_server(debug=debug, port=8090) else: - raise NotImplementedError(f"doesnt support mode={mode}") + raise NotImplementedError(f"App doesn't support mode='{self.mode}'.") class SimulationDataApp(App): """App for viewing contents of a :class:`.SimulationData` instance.""" - sim_data: td.SimulationData = pd.Field( + sim_data: SimulationData = pd.Field( ..., title="Simulation data", description="A :class:`.SimulationData` instance to view." ) - def _make_layout(self, app: DASH_APP) -> "Dash Layout": + def _make_layout(self, app: DASH_APP) -> dcc.Tabs: """Creates the layout for the app.""" layout = dcc.Tabs([]) @@ -91,12 +87,17 @@ def _make_layout(self, app: DASH_APP) -> "Dash Layout": data_plotly = DataPlotly.from_monitor_data( monitor_data=monitor_data, monitor_name=monitor_name ) + if data_plotly is None: + continue component = data_plotly.make_component(app) layout.children += [component] # log component = dcc.Tab( - [html.Div([html.Code(self.sim_data.log, style={"whiteSpace": "pre-wrap"})])], + [ + html.Div([html.H1("Solver Log")]), + html.Div([html.Code(self.sim_data.log, style={"whiteSpace": "pre-wrap"})]), + ], label="log", ) layout.children += [component] @@ -104,25 +105,26 @@ def _make_layout(self, app: DASH_APP) -> "Dash Layout": return layout @classmethod - def from_file(cls, path: str, mode: APP_MODE = DEFAULT_MODE): - """Load the SimulationDataApp from a tidy3d data file in .hdf5 format.""" - sim_data = td.SimulationData.from_file(path) - return cls(sim_data=sim_data, mode=mode) + def from_file(cls, path: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ + """Load the :class:`.SimulationDataApp` from a tidy3d data file in .hdf5 format.""" + sim_data = SimulationData.from_file(path) + sim_data_normalized = sim_data.normalize() + return cls(sim_data=sim_data_normalized, mode=mode) class SimulationApp(App): """TODO: App for viewing and editing a :class:`.Simulation`.""" - simulation: td.Simulation = pd.Field( - ..., title="Simulation", description="A Simulation instance to view." + simulation: Simulation = pd.Field( + ..., title="Simulation", description="A :class:`.Simulation` instance to view." ) - def _make_layout(self, app: DASH_APP) -> "Dash Layout": + def _make_layout(self, app: DASH_APP) -> dcc.Tabs: """Creates the layout for the app.""" - return cc.Tabs([]) + return dcc.Tabs([]) @classmethod - def from_file(cls, path: str, mode: APP_MODE = DEFAULT_MODE): + def from_file(cls, path: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ """Load the SimulationApp from a tidy3d Simulation file in .json or .yaml format.""" - simulation = td.Simulation.from_file(path) + simulation = Simulation.from_file(path) return cls(simulation=simulation, mode=mode) diff --git a/tidy3d/plugins/plotly/component.py b/tidy3d/plugins/plotly/component.py index ac4d7fa806..cf912e3985 100644 --- a/tidy3d/plugins/plotly/component.py +++ b/tidy3d/plugins/plotly/component.py @@ -1,10 +1,9 @@ +"""Defines the structure of the components displayed in the app tabs.""" from abc import ABC, abstractmethod from dash import dcc -import sys -sys.path.append("../../../") - -from tidy3d.components.base import Tidy3dBaseModel +from .utils import PlotlyFig +from ...components.base import Tidy3dBaseModel """ How the components work. @@ -24,15 +23,13 @@ """ -class UIComponent(Tidy3dBaseModel): - @abstractmethod - def make_component(self, app) -> dcc.Tab: - """Creates the dash component for this montor data.""" +class UIComponent(Tidy3dBaseModel, ABC): + """Base class for a UI component. Individual data storage wrappers override this.""" @abstractmethod - def make_figure(self) -> "Plotly.Figure": + def make_component(self, app) -> dcc.Tab: """Creates the dash component for this montor data.""" @abstractmethod - def plotly(self, **kwargs) -> "Plotly.Figure": - """Make a plotly figure using all of the optional kwargs.""" + def make_figure(self) -> PlotlyFig: + """Creates a plotly figure for this component given its current state.""" diff --git a/tidy3d/plugins/plotly/data.py b/tidy3d/plugins/plotly/data.py index b3aa4188af..6de73db30c 100644 --- a/tidy3d/plugins/plotly/data.py +++ b/tidy3d/plugins/plotly/data.py @@ -1,163 +1,529 @@ -from abc import ABC, abstractmethod -from typing import Union +"""Defines how the specific data objects render as UI components.""" +from abc import ABC +from typing import Union, Tuple, List from typing_extensions import Literal import numpy as np import plotly.graph_objects as go -import plotly.express as px -from dash import dcc, html, Output, Input +from dash import dcc, html, Output, Input, Dash +import pydantic as pd from .component import UIComponent +from .utils import PlotlyFig -import sys +from ...components.data import FluxData, FluxTimeData, FieldData, FieldTimeData +from ...components.data import ModeFieldData, ModeData, AbstractScalarFieldData +from ...components.geometry import Geometry +from ...components.types import Axis, Direction +from ...log import Tidy3dKeyError, log -sys.path.append("../../../") +# divide an SI unit by these to put in the unit in the name +# eg. 1e-12 (sec) / PICOSECOND = 1 (ps) +# multiply a value in the specified unit by these to put in SI +# eg. 1 (THz) * TERAHERTZ = 1e12 (Hz) +PICOSECOND = 1e-12 +TERAHERTZ = 1e12 -from tidy3d.components.data import FluxData, FluxTimeData, FieldData, Tidy3dData, FieldTimeData -from tidy3d.components.geometry import Geometry -from tidy3d.components.types import Axis -from tidy3d.log import Tidy3dKeyError +# supported data types +Tidy3dDataType = Union[FluxData, FluxTimeData, FieldData, FieldTimeData, ModeData, ModeFieldData] -class DataPlotly(UIComponent): +class DataPlotly(UIComponent, ABC): + """Base class for anything that generates dash components from tidy3d data objects.""" - monitor_name: str - data: Union[FluxData, FluxTimeData, FieldData, FieldTimeData] + monitor_name: str = pd.Field(..., tite="monitor name", description="Name of the monitor.") + data: Tidy3dDataType = pd.Field( + ..., tite="data", description="The Tidy3dData object wrapped by this UI component." + ) @property def label(self) -> str: - """Get tab label for component.""" - return f"monitor: {self.monitor_name}" + """Get tab label for this component.""" + return f"monitor: '{self.monitor_name}'" - @property - def id(self) -> str: - """Get unique id for component.""" - return f"monitor_{self.monitor_name}" + # @property + # def unique_id(self) -> str: + # """Get unique id for this component.""" + # return f"monitor_{self.monitor_name}" @staticmethod - def sel_by_val(data, val): - if "re" in val.lower(): + def sel_by_val(data, val: str) -> "data": + """Select the correct data type based on the `val` selection.""" + if val.lower() == "real": return data.real - if "im" in val.lower(): + if val.lower() == "imag": return data.imag - if "abs" in val.lower(): + if val.lower() == "abs": return abs(data) + if val.lower() == "abs^2": + return abs(data) ** 2 + raise ValueError(f"Could not find the right function to apply with {val}.") - @property - def ft_label_coords(self): - """Get the `freq` or `time` label and coords.""" - - scalar_field_data = self.data.data_dict[self.field_val] - - if "f" in self.scalar_field_data.data.coords: - ft_label = "freq" - ft_coords = scalar_field_data.data.coords["f"].values - elif "t" in self.scalar_field_data.data.coords: - ft_label = "time" - ft_coords = scalar_field_data.data.coords["t"].values - else: - raise Tidy3dKeyError(f"neither frequency nor time data found in this data object.") - - return ft_label, ft_coords - - def append_monitor_name(self, value): - """makes the ids unique for this element.""" + def append_monitor_name(self, value: str) -> str: + """Adds the monitor name to a value, used to make the ids unique across all monitors.""" return f"{value}_{self.monitor_name}" @classmethod - def from_monitor_data(self, monitor_name: str, monitor_data): + def from_monitor_data(cls, monitor_name: str, monitor_data: Tidy3dDataType) -> "cls": """Load a PlotlyData UI component from the monitor name and its data.""" - DATA_PLOTLY_MAP = { + # maps the supplied ``monitor_data`` argument to the corresponding plotly wrapper. + data_plotly_map = { FluxData: FluxDataPlotly, FluxTimeData: FluxTimeDataPlotly, FieldData: FieldDataPlotly, FieldTimeData: FieldTimeDataPlotly, + ModeFieldData: ModeFieldDataPlotly, + ModeData: ModeDataPlotly, } + # get the type of the supplied ``monitor_data``. monitor_data_type = type(monitor_data) - PlotlyDataType = DATA_PLOTLY_MAP.get(monitor_data_type) - if not PlotlyDataType: - raise Tidy3dKeyError(f"could not find the monitor data type: {monitor_data_type}.") + # try to grab the right plotly wrapper and complain (or skip) if not found. + plotly_data_type = data_plotly_map.get(monitor_data_type) + if not plotly_data_type: + log.warning( + f"could not find a plotly wrapper for monitor {monitor_name}" + f"of type {monitor_data_type.__name__}" + ) + return None - return PlotlyDataType(data=monitor_data, monitor_name=monitor_name) + # return the right component + return plotly_data_type(data=monitor_data, monitor_name=monitor_name) -class AbstractFluxDataPlotly(DataPlotly): +class AbstractFluxDataPlotly(DataPlotly, ABC): """Flux data in frequency or time domain.""" - data: Union[FluxData, FluxTimeData] + data: Union[FluxData, FluxTimeData] = pd.Field( + ..., + title="data", + description="A flux data object in freq or time domain.", + ) + + @property + def ft_label_coords_units(self) -> Tuple[str, List[float], str]: + """Get the ``freq`` or ``time`` label, coords, and units depending on the data contents.""" + + if "f" in self.data.data.coords: + ft_label = "freq" + ft_coords = self.data.data.coords["f"].values / TERAHERTZ + ft_units = "THz" + elif "t" in self.data.data.coords: + ft_label = "time" + ft_coords = self.data.data.coords["t"].values / PICOSECOND + ft_units = "ps" + else: + raise Tidy3dKeyError("Neither frequency nor time data found in this data object.") - def make_figure(self): + return ft_label, ft_coords, ft_units + + def make_figure(self) -> PlotlyFig: """Generate plotly figure from the current state of self.""" return self.plotly() - def make_component(self, app) -> dcc.Tab: + def make_component(self, app: Dash) -> dcc.Tab: """Creates the dash component for this montor data.""" - component = dcc.Tab( + # initital setup + fig = self.make_figure() + + # individual components + flux_plot = html.Div( [ - html.Div( - [ - dcc.Graph( - id=self.append_monitor_name("figure"), - figure=self.make_figure(), - ) - ], - style={"padding": 10, "flex": 1}, + dcc.Graph( + id=self.append_monitor_name("figure"), + figure=fig, ) ], + style={"padding": 10, "flex": 1}, + ) + + # define layout + component = dcc.Tab( + [ + html.H1(f"Viewing data for {type(self.data).__name__}: '{self.monitor_name}'"), + flux_plot, + ], label=self.label, ) - def plotly(self): - ft_label, ft_coords = self.ft_label_coords - return px.line(x=ft_coords, y=self.data.data.values) + # return the layout of the component so the app can insert it + return component + + def plotly(self) -> PlotlyFig: + """Generate the plotly figure for this component.""" + + ft_label, ft_coords, ft_units = self.ft_label_coords_units + ft_units = "THz" if "f" in ft_label else "ps" + + fig = go.Figure(go.Scatter(x=ft_coords, y=self.data.data.values)) + + fig.update_layout( + xaxis_title=f"{ft_label} ({ft_units})", + yaxis_title="Flux (normalized)", + ) + return fig class FluxDataPlotly(AbstractFluxDataPlotly): """Flux in frequency domain.""" - data: FluxData + data: FluxData = pd.Field( + ..., + title="data", + description="A flux data object in the frequency domain.", + ) class FluxTimeDataPlotly(AbstractFluxDataPlotly): """Flux in time domain.""" - data: FluxTimeData + data: FluxTimeData = pd.Field( + ..., + title="data", + description="A flux data object in the time domain.", + ) + +class ModeDataPlotly(DataPlotly): + """Flux data in frequency or time domain.""" + + data: ModeData = pd.Field( + ..., + title="data", + description="A mode amplitude data object", + ) + + amps_or_neff: Literal["amps", "neff"] = pd.Field( + "amps", + title="Amps or effective index value", + description="The state of the component's 'amplitude or neff' value.", + ) + + val: str = pd.Field( + "abs^2", + title="Plotting value value", + description="The state of the component's plotting value value.", + ) -class AbstractFieldDataPlotly(DataPlotly): - """Field data in frequency or time domain.""" + dir_val: Direction = pd.Field( + None, title="Direction value", description="The state of the component's direction value." + ) - data: Union[FieldData, FieldTimeData] - field_val: str = "Ex" - cs_axis: Axis = 0 - cs_val: float = None - val: str = "abs" - ft_val: float = None + mode_ind_val: int = pd.Field( + None, title="Mode index value", description="The state of the component's mode index value." + ) @property - def scalar_field_data(self): + def mode_ind_coords(self) -> List[float]: + """Get the mode indices.""" + return self.data.amps.coords["mode_index"].values + + @property + def direction_coords(self) -> List[float]: + """Get the mode indices.""" + return self.data.amps.coords["direction"].values + + @property + def ft_label_coords_units(self) -> Tuple[str, List[float], str]: + """Get the `freq` or `time` label and coords.""" + + ft_label = "freq" + ft_coords = self.data.amps.coords["f"].values / TERAHERTZ + ft_units = "THz" + return ft_label, ft_coords, ft_units + + @property + def dir_dropdown_hidden(self) -> bool: + """Should the dropdown be hidden?""" + return self.amps_or_neff == "neff" + + def make_figure(self) -> PlotlyFig: + """Generate plotly figure from the current state of self.""" + + if self.dir_val is None: + self.dir_val = self.direction_coords[0] + + if self.mode_ind_val is None: + self.mode_ind_val = self.mode_ind_coords[0] + + if self.amps_or_neff == "amps": + return self.plotly_amps(mode_index=self.mode_ind_val, dir_val=self.dir_val) + + return self.plotly_neff(mode_index=self.mode_ind_val) + + def make_component(self, app: Dash) -> dcc.Tab: + """Creates the dash component for this montor data.""" + + # initital setup + fig = self.make_figure() + + # individual components + + # amplitude plot + amp_plot = html.Div( + [ + dcc.Graph( + id=self.append_monitor_name("figure"), + figure=fig, + ) + ], + style={"padding": 10, "flex": 1}, + ) + + # select amps or neff + amps_or_neff_dropdown = dcc.Dropdown( + options=["amps", "neff"], + value=self.amps_or_neff, + id=self.append_monitor_name("amps_or_neff_dropdown"), + ) + + # select real, abs, imag, power + field_value_dropdown = dcc.Dropdown( + options=["real", "imag", "abs^2"], + value=self.val, + id=self.append_monitor_name("val_dropdown"), + ) + + # header for direction dropdown + dir_dropdown_header = html.Div( + [html.H2("Direction: forward (+) or backward (-).")], + hidden=self.dir_dropdown_hidden, + id=self.append_monitor_name("dir_dropdown_header"), + ) + + # select direction + dir_value_dropdown = html.Div( + [ + dcc.Dropdown( + options=list(self.direction_coords), + value=self.dir_val, + id=self.append_monitor_name("dir_dropdown"), + ) + ], + hidden=self.dir_dropdown_hidden, + id=self.append_monitor_name("dir_dropdown_div"), + ) + + # mode index selector + mode_ind_dropdown = html.Div( + [ + dcc.Dropdown( + options=list(self.mode_ind_coords), + value=self.mode_ind_val, + id=self.append_monitor_name("mode_index_selector"), + ), + ] + ) + + # data control panel + panel_children = [ + html.H2("Amplitude or effective index."), + amps_or_neff_dropdown, + html.H2("Value to plot."), + field_value_dropdown, + dir_dropdown_header, + dir_value_dropdown, + html.H2("Mode Index."), + mode_ind_dropdown, + ] + + plot_selections = html.Div(panel_children, style={"padding": 10, "flex": 1}) + + # define layout + component = dcc.Tab( + [ + html.H1(f"Viewing data for {type(self.data).__name__}: '{self.monitor_name}'"), + html.Div( + [ + # left hand side + amp_plot, + # right hand side + plot_selections, + ], + # make elements in above list stack row-wise + style={"display": "flex", "flex-direction": "row"}, + ), + ], + label=self.label, + ) + + # link what happens in the inputs to what gets displayed in the figure + @app.callback( + Output(self.append_monitor_name("figure"), "figure"), + [ + Input(self.append_monitor_name("amps_or_neff_dropdown"), "value"), + Input(self.append_monitor_name("val_dropdown"), "value"), + Input(self.append_monitor_name("dir_dropdown"), "value"), + Input(self.append_monitor_name("mode_index_selector"), "value"), + ], + ) + def set_field(value_amps_or_neff, value_val, value_dir, value_mode_ind): + self.amps_or_neff = str(value_amps_or_neff) + self.val = str(value_val) + self.dir_val = str(value_dir) + + self.mode_ind_val = int(value_mode_ind) if value_mode_ind is not None else None + fig = self.make_figure() + return fig + + @app.callback( + Output(self.append_monitor_name("dir_dropdown_header"), "hidden"), + [ + Input(self.append_monitor_name("amps_or_neff_dropdown"), "value"), + ], + ) + def set_dir_header_visibilty(value_amps_or_neff): + self.amps_or_neff = str(value_amps_or_neff) + return self.dir_dropdown_hidden + + @app.callback( + Output(self.append_monitor_name("dir_dropdown_div"), "hidden"), + [ + Input(self.append_monitor_name("amps_or_neff_dropdown"), "value"), + ], + ) + def set_dir_dropdown_visibilty(value_amps_or_neff): + self.amps_or_neff = str(value_amps_or_neff) + return self.dir_dropdown_hidden + + return component + + def plotly_amps(self, mode_index: int, dir_val: str): + """Make a line chart for the mode amplitudes.""" + + ft_label, ft_coords, ft_units = self.ft_label_coords_units + ft_units = "THz" if "f" in ft_label else "ps" + + amp_val = self.sel_by_val(self.data.amps, val=self.val) + amp_val = amp_val.sel(direction=dir_val, mode_index=mode_index) + fig = go.Figure(go.Scatter(x=ft_coords, y=amp_val)) + + fig.update_layout( + title_text=f"amplitudes of mode w/ index {self.mode_ind_val} in {dir_val} direction.", + title_x=0.5, + xaxis_title=f"{ft_label} ({ft_units})", + yaxis_title="Amplitude (normalized)", + ) + return fig + + def plotly_neff(self, mode_index: int): + """Make a line chart for the mode amplitudes.""" + + ft_label, ft_coords, ft_units = self.ft_label_coords_units + ft_units = "THz" if "f" in ft_label else "ps" + + neff_val = self.sel_by_val(self.data.n_complex, val=self.val) + neff_val = neff_val.sel(mode_index=mode_index) + fig = go.Figure(go.Scatter(x=ft_coords, y=neff_val)) + + fig.update_layout( + title_text=f"effective index of mode w/ index {self.mode_ind_val}", + title_x=0.5, + xaxis_title=f"{ft_label} ({ft_units})", + yaxis_title="Effective index", + ) + return fig + + +class AbstractFieldDataPlotly(DataPlotly, ABC): + """Some kind of field-like data plotted in the app.""" + + data: Union[FieldData, FieldTimeData] = pd.Field( + ..., + title="data", + description="A Field-like data object.", + ) + + field_val: str = pd.Field( + None, title="Field value", description="The component's field component value." + ) + + val: str = pd.Field( + "abs", title="Plot value value", description="The component's plotting value value." + ) + + cs_axis: Axis = pd.Field( + 0, title="Cross section axis value", description="The component's cross section axis value." + ) + + cs_val: float = pd.Field( + None, + title="Cross section position value", + description="The component's cross section position value.", + ) + + ft_val: float = pd.Field( + None, title="Freq or time value", description="The component's frequency or time value." + ) + + mode_ind_val: int = pd.Field( + None, title="Mode index value", description="The component's mode index value value." + ) + + @property + def ft_label_coords_units(self) -> Tuple[str, List[float], str]: + """Get the `freq` or `time` label and coords.""" + + scalar_field_data = self.data.data_dict[self.field_val] + + if "f" in self.scalar_field_data.data.coords: + ft_label = "freq" + ft_coords = scalar_field_data.data.coords["f"].values / TERAHERTZ + ft_units = "THz" + elif "t" in self.scalar_field_data.data.coords: + ft_label = "time" + ft_coords = scalar_field_data.data.coords["t"].values / PICOSECOND + ft_units = "ps" + else: + raise Tidy3dKeyError("neither frequency nor time data found in this data object.") + + return ft_label, ft_coords, ft_units + + @property + def inital_field_val(self): + """The starting field value.""" + field_vals = list(self.data.data_dict.keys()) + if len(field_vals) == 0: + raise ValueError("Data doesn't have any field components stored.") + return field_vals[0] + + @property + def scalar_field_data(self) -> AbstractScalarFieldData: """The current scalar field monitor data.""" + if self.field_val is None: + self.field_val = self.inital_field_val + return self.data.data_dict[self.field_val] @property - def xyz_label_coords(self): + def xyz_label_coords(self) -> Tuple[str, List[float]]: """Get the plane normal direction label and coords.""" xyz_label = "xyz"[self.cs_axis] xyz_coords = self.scalar_field_data.data.coords[xyz_label].values return xyz_label, xyz_coords - def make_figure(self): + @property + def mode_ind_coords(self) -> List[int]: + """Get the mode indices.""" + return self.scalar_field_data.coords["mode_index"].values + + def make_figure(self) -> PlotlyFig: """Generate plotly figure from the current state of self.""" + # if no field specified, use the first one in the fields list + if self.field_val is None: + self.field_val = self.inital_field_val + + # if cross section value, use the average of the coordinates of the current axis xyz_label, xyz_coords = self.xyz_label_coords if self.cs_val is None: self.cs_val = np.mean(xyz_coords) - ft_label, ft_coords = self.ft_label_coords + # if no freq or time value, use the average of the coordinates + ft_label, ft_coords, _ = self.ft_label_coords_units if self.ft_val is None: self.ft_val = np.mean(ft_coords) @@ -170,83 +536,152 @@ def make_figure(self): return self.plotly(**plotly_kwargs) - def make_component(self, app): + def make_component(self, app: Dash) -> dcc.Tab: # pylint:disable=too-many-locals """Creates the dash component.""" # initial setup xyz_label, xyz_coords = self.xyz_label_coords + ft_label, ft_coords, ft_units = self.ft_label_coords_units + fig = self.make_figure() + + # individual components + # plot of the fields + field_plot = html.Div( + [ + dcc.Graph( + id=self.append_monitor_name("figure"), + figure=fig, + ) + ], + style={"padding": 10, "flex": 1}, + ) + + # pick the field component + field_dropdown = dcc.Dropdown( + options=list(self.data.data_dict.keys()), + value=self.field_val, + id=self.append_monitor_name("field_dropdown"), + ) + + # pick the real / imag / abs to plot + field_value_dropdown = dcc.Dropdown( + options=["real", "imag", "abs"], + value=self.val, + id=self.append_monitor_name("val_dropdown"), + ) + + # pick the cross section axis + xyz_dropdown = dcc.Dropdown( + options=["x", "y", "z"], + value=xyz_label, + id=self.append_monitor_name("cs_axis_dropdown"), + ) + + # pick the cross section position + xyz_slider = dcc.Slider( + min=xyz_coords[0], + max=xyz_coords[-1], + value=self.cs_val, + id=self.append_monitor_name("cs_slider"), + ) + + # combine the cross section selection into one component + xyz_selection = html.Div([xyz_dropdown, xyz_slider]) + + # pick the frequency or time + freq_time_slider = html.Div( + [ + dcc.Slider( + min=ft_coords[0], + max=ft_coords[-1], + value=self.ft_val, + id=self.append_monitor_name("ft_slider"), + ), + ] + ) + + # all the controls for adjusting plotted data + plot_selections = html.Div( + [ + html.H2("Field component."), + field_dropdown, + html.H2("Value to plot."), + field_value_dropdown, + html.H2("Cross-section axis and position."), + xyz_selection, + html.H2(f"{ft_label} value ({ft_units})."), + freq_time_slider, + ], + style={"padding": 10, "flex": 1}, + ) + + # add a mode index dropdown to right hand side, if applicable + if self.mode_ind_val is not None: + + # make a mode index dropdown + mode_ind_dropdown = html.Div( + [ + dcc.Dropdown( + options=list(self.mode_ind_coords), + value=self.mode_ind_val, + id=self.append_monitor_name("mode_index_selector"), + ), + ] + ) + + # add this to the plot selections panel component + plot_selections.children.append(mode_ind_dropdown) + + # full layout component = dcc.Tab( [ + # title + html.H1(f"Viewing data for {type(self.data).__name__}: '{self.monitor_name}'"), + # below title html.Div( [ - html.Div( - [ - dcc.Graph( - id=self.append_monitor_name("figure"), - figure=self.make_figure(), - ) - ], - style={"padding": 10, "flex": 1}, - ), - html.Div( - [ - html.H1(f"Viewing data for FieldMonitor: {self.monitor_name}"), - html.H2(f"Field component."), - dcc.Dropdown( - options=list(self.data.data_dict.keys()), - value=self.field_val, - id=self.append_monitor_name("field_dropdown"), - ), - html.H2(f"Value to plot."), - dcc.Dropdown( - options=["real", "imag", "abs"], - value=self.val, - id=self.append_monitor_name("val_dropdown"), - ), - html.Br(), - html.Div( - [ - dcc.Dropdown( - options=["x", "y", "z"], - value=xyz_label, - id=self.append_monitor_name("cs_axis_dropdown"), - ), - dcc.Slider( - min=xyz_coords[0], - max=xyz_coords[-1], - value=np.mean(xyz_coords), - id=self.append_monitor_name("cs_slider"), - ), - ], - ), - ], - style={"padding": 10, "flex": 1}, - ), + # left hand side + field_plot, + # right hand side + plot_selections, ], + # make elements in above list stack row-wise style={"display": "flex", "flex-direction": "row"}, - ) + ), ], + # label for the tab label=self.label, ) - @app.callback( - Output(self.append_monitor_name("figure"), "figure"), - [ - Input(self.append_monitor_name("field_dropdown"), "value"), - Input(self.append_monitor_name("val_dropdown"), "value"), - Input(self.append_monitor_name("cs_axis_dropdown"), "value"), - Input(self.append_monitor_name("cs_slider"), "value"), - ], - ) - def set_field(value_field, value_val, value_cs_axis, value_cs): + # these are the inputs to the callback function which links the buttons to the figure + app_inputs = [ + Input(self.append_monitor_name("field_dropdown"), "value"), + Input(self.append_monitor_name("val_dropdown"), "value"), + Input(self.append_monitor_name("cs_axis_dropdown"), "value"), + Input(self.append_monitor_name("cs_slider"), "value"), + Input(self.append_monitor_name("ft_slider"), "value"), + ] + + # add the mode index dropdown to the app inputs, if defined + if self.mode_ind_val is not None: + app_inputs.append(Input(self.append_monitor_name("mode_index_selector"), "value")) + + # link what happens in the app_inputs to what gets displayed in the figure + @app.callback(Output(self.append_monitor_name("figure"), "figure"), app_inputs) + def set_field( # pylint:disable=too-many-arguments + value_field, value_val, value_cs_axis, value_cs, value_ft, value_mode_ind=None + ): self.field_val = str(value_field) self.val = str(value_val) self.cs_axis = ["x", "y", "z"].index(value_cs_axis) self.cs_val = float(value_cs) + self.ft_val = float(value_ft) + self.mode_ind_val = int(value_mode_ind) if value_mode_ind is not None else None fig = self.make_figure() return fig + # set the minimum of the xyz sliderbar depending on the cross-section axis @app.callback( Output(self.append_monitor_name("cs_slider"), "min"), Input(self.append_monitor_name("cs_axis_dropdown"), "value"), @@ -256,52 +691,123 @@ def set_min(value_cs_axis): _, xyz_coords = self.xyz_label_coords return xyz_coords[0] + # set the xyz slider back to the average if the axis changes. + @app.callback( + Output(self.append_monitor_name("cs_slider"), "value"), + Input(self.append_monitor_name("cs_axis_dropdown"), "value"), + ) + def reset_slider_position(value_cs_axis): + self.cs_axis = ["x", "y", "z"].index(value_cs_axis) + _, xyz_coords = self.xyz_label_coords + self.cs_val = float(np.mean(xyz_coords)) + return self.cs_val + + # set the maximum of the xyz sliderbar depending on the cross-section axis @app.callback( Output(self.append_monitor_name("cs_slider"), "max"), Input(self.append_monitor_name("cs_axis_dropdown"), "value"), ) - def set_max(value_cs_axis): # , value_f): + def set_max(value_cs_axis): self.cs_axis = ["x", "y", "z"].index(value_cs_axis) _, xyz_coords = self.xyz_label_coords return xyz_coords[-1] return component - def plotly( - self, field: str, freq: float, val: Literal["real", "imag", "abs"], x=None, y=None, z=None - ): + def plotly( # pylint:disable=too-many-arguments, too-many-locals + self, + field: str, + val: Literal["real", "imag", "abs"], + freq: float = None, + time: float = None, + x: float = None, + y: float = None, + z: float = None, + mode_index: int = None, + ) -> PlotlyFig: """Creates the plotly figure given some parameters.""" axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z) - scalar_field_data = self.data.data_dict[field] - sel_freq = scalar_field_data.data.sel(f=freq) + + # grab by field name + scalar_field_data = self.scalar_field_data + + # select mode_index, if given + if mode_index is not None: + scalar_field_data = scalar_field_data.sel(mode_index=mode_index) + + # select by frequency, if given + if freq is not None: + freq *= TERAHERTZ + sel_ft = scalar_field_data.data.sel(f=freq, method="nearest") + + # select by time, if given + if time is not None: + time *= PICOSECOND + sel_ft = scalar_field_data.data.sel(f=freq, method="nearest") + + # select the cross sectional plane data xyz_labels = ["x", "y", "z"] xyz_kwargs = {xyz_labels.pop(axis): position} - sel_xyz = sel_freq.interp(**xyz_kwargs) + sel_xyz = sel_ft.sel(**xyz_kwargs, method="nearest") + + # get the correct field value (real, imaginary, abs) sel_val = self.sel_by_val(data=sel_xyz, val=val) - d1 = sel_val.coords[xyz_labels[0]] - d2 = sel_val.coords[xyz_labels[1]] + + # get the correct x and y labels + coords_plot_x = sel_val.coords[xyz_labels[0]] + coords_plot_y = sel_val.coords[xyz_labels[1]] + + # construct the field plot fig = go.Figure( data=go.Heatmap( - x=d1, - y=d2, + x=coords_plot_x, + y=coords_plot_y, z=sel_val.values, transpose=True, type="heatmap", colorscale="magma" if val in "abs" in val else "RdBu", ) ) - fig.update_layout(title=f'{val}[{field}({"xyz"[axis]}={position:.2e}, f={freq:.2e})]') + + # update title and x and y labels. + ft_text = f"f={freq:.2e}" if freq is not None else f"t={time:.2e}" + _, (xlabel, ylabel) = Geometry.pop_axis("xyz", axis=axis) + fig.update_layout( + title_text=f'{val}[{field}({"xyz"[axis]}={position:.2e}, {ft_text})]', + title_x=0.5, + xaxis_title=f"{xlabel} (um)", + yaxis_title=f"{ylabel} (um)", + ) + return fig class FieldDataPlotly(AbstractFieldDataPlotly): - """Plot FieldData in app.""" + """Plot :class:`.FieldData` in app.""" - data: FieldData + data: FieldData = pd.Field( + ..., + title="data", + description="A field data object in the frequency domain", + ) class FieldTimeDataPlotly(AbstractFieldDataPlotly): - """Plot FieldTimeData in app.""" + """Plot :class:`.FieldTimeData` in app.""" + + data: FieldTimeData = pd.Field( + ..., + title="data", + description="A field data object in the time domain.", + ) + + +class ModeFieldDataPlotly(AbstractFieldDataPlotly): + """Plot :class:`.ModeFieldData` in app.""" - data: FieldTimeData + data: ModeFieldData = pd.Field( + ..., + title="data", + description="A mode field object.", + ) diff --git a/tidy3d/plugins/plotly/simulation.py b/tidy3d/plugins/plotly/simulation.py index 0feeca4d53..057311b5fc 100644 --- a/tidy3d/plugins/plotly/simulation.py +++ b/tidy3d/plugins/plotly/simulation.py @@ -1,4 +1,4 @@ -"""Plotly wrapper for tidy3d.""" +"""Simulation and Geometry plotting with plotly.""" # pylint:disable=too-many-arguments, protected-access from typing import Tuple @@ -9,6 +9,7 @@ from .utils import PlotlyFig, add_fig_if_none, equal_aspect_plotly, plot_params_sim_boundary from .component import UIComponent + from ...components.types import Axis from ...components.simulation import Simulation from ...components.structure import Structure @@ -26,7 +27,8 @@ def plotly_shape( ) -> PlotlyFig: """Plot a shape to a figure.""" _shape = Geometry.evaluate_inf_shape(shape) - (xs, ys), _ = Geometry.strip_coords(shape=_shape) + exterior_coords, _ = Geometry.strip_coords(shape=_shape) + xs, ys = list(zip(*exterior_coords)) plotly_trace = go.Scatter( x=xs, y=ys, @@ -120,35 +122,43 @@ def make_figure(self): plotly_kwargs = {xyz_label: self.cs_val} return self.plotly(**plotly_kwargs) - def make_component(self, app): + def make_component(self, app): # pylint: disable=too-many-locals """Creates the dash component.""" xyz_label, (xyz_min, xyz_max) = self.xyz_label_bounds figure = self.make_figure() - xyz_slider = html.Div( - [ - dcc.Dropdown( - options=["x", "y", "z"], - value=xyz_label, - id="simulation_cs_axis_dropdown", - ), - dcc.Slider( - min=xyz_min, - max=xyz_max, - value=self.cs_val, - id="simulation_cs_slider", - ), - ], - style={"padding": 50, "flex": 1}, - ) - graph = html.Div( [dcc.Graph(figure=figure, id="simulation_plot")], style={"padding": 10, "flex": 1} ) + xyz_header = html.H2("Cross-section axis and position.") + + xyz_dropdown = dcc.Dropdown( + options=["x", "y", "z"], + value=xyz_label, + id="simulation_cs_axis_dropdown", + ) + + xyz_slider = dcc.Slider( + min=xyz_min, + max=xyz_max, + value=self.cs_val, + id="simulation_cs_slider", + ) + + xyz_selection = html.Div( + [xyz_header, xyz_dropdown, xyz_slider], + style={"padding": 50, "flex": 1}, + ) + component = dcc.Tab( - [html.Div([graph, xyz_slider], style={"display": "flex", "flex-direction": "row"})], + [ + html.H1("Viewing Simulation."), + html.Div( + [graph, xyz_selection], style={"display": "flex", "flex-direction": "row"} + ), + ], label="Simulation", ) @@ -159,11 +169,22 @@ def make_component(self, app): Input("simulation_cs_slider", "value"), ], ) - def set_xyz_sliderbar(cs_axis_string, cs_val): + def set_fig_from_xyz_sliderbar(cs_axis_string, cs_val): self.cs_axis = ["x", "y", "z"].index(cs_axis_string) self.cs_val = float(cs_val) return self.make_figure() + # set the xyz slider back to the average if the axis changes. + @app.callback( + Output("simulation_cs_slider", "value"), + Input("simulation_cs_axis_dropdown", "value"), + ) + def reset_slider_position(value_cs_axis): + self.cs_axis = ["x", "y", "z"].index(value_cs_axis) + _, (xyz_min, xyz_max) = self.xyz_label_bounds + self.cs_val = float((xyz_min + xyz_max) / 2.0) + return self.cs_val + @app.callback( Output("simulation_cs_slider", "min"), Input("simulation_cs_axis_dropdown", "value"), @@ -309,8 +330,10 @@ def _plotly_shape_structure( self, medium: Medium, shape: ShapelyGeo, fig: PlotlyFig ) -> PlotlyFig: """Plot a structure's cross section shape for a given medium.""" - plot_params_struct = self.simulation._get_structure_plot_params(medium=medium) mat_index = self.simulation.medium_map[medium] + plot_params_struct = self.simulation._get_structure_plot_params( + medium=medium, mat_index=mat_index + ) name = medium.name if medium.name else f"medium[{mat_index}]" fig = plotly_shape(shape=shape, plot_params=plot_params_struct, fig=fig, name=name) return fig @@ -536,9 +559,10 @@ def _plotly_cleanup( _, (xlabel, ylabel) = self.simulation.pop_axis("xyz", axis=normal_axis) fig.update_layout( - title=f'{"xyz"[normal_axis]} = {pos:.2f}', - xaxis_title=rf"{xlabel} ($\mu m$)", - yaxis_title=rf"{ylabel} ($\mu m$)", + title_text=f'Simulation plotted on {"xyz"[normal_axis]} = {pos:.2f} cross-section', + title_x=0.5, + xaxis_title=f"{xlabel} (um)", + yaxis_title=f"{ylabel} (um)", legend_title="Contents", ) @@ -581,7 +605,7 @@ def _plotly_clean_labels(fig: PlotlyFig) -> PlotlyFig: seen = [] for trace in fig["data"]: name = trace["name"] - if name not in seen: + if name is not None and name not in seen: seen.append(name) else: trace["showlegend"] = False diff --git a/tidy3d/plugins/plotly/utils.py b/tidy3d/plugins/plotly/utils.py index ab37e0fb4b..c4409a6581 100644 --- a/tidy3d/plugins/plotly/utils.py +++ b/tidy3d/plugins/plotly/utils.py @@ -1,3 +1,4 @@ +"""Utilities for plotly plotting.""" from functools import wraps import plotly.graph_objects as go diff --git a/tidy3d/plugins/smatrix/smatrix.py b/tidy3d/plugins/smatrix/smatrix.py index 1f46f2fd69..008a2af5b4 100644 --- a/tidy3d/plugins/smatrix/smatrix.py +++ b/tidy3d/plugins/smatrix/smatrix.py @@ -5,8 +5,13 @@ import pydantic as pd import numpy as np -from ... import Simulation, Box, ModeSpec, ModeMonitor, ModeSource, GaussianPulse, SimulationData from ...constants import HERTZ, C_0 +from ...components.simulation import Simulation +from ...components.geometry import Box +from ...components.mode import ModeSpec +from ...components.monitor import ModeMonitor +from ...components.source import ModeSource, GaussianPulse +from ...components.data import SimulationData from ...components.types import Direction, Ax from ...components.viz import add_ax_if_none, equal_aspect from ...components.base import Tidy3dBaseModel @@ -92,7 +97,7 @@ class ComponentModeler(Tidy3dBaseModel): def _sim_has_no_sources(cls, val): """Make sure simulation has no sources as they interfere with tool.""" if len(val.sources) > 0: - raise SetupError(f"Simulation must not have sources.") + raise SetupError("Simulation must not have sources.") return val def _to_monitor(self, port: Port) -> ModeMonitor: @@ -261,5 +266,5 @@ def load(self) -> SMatrixType: """Load an Smatrix from saved BatchData object.""" if self.batch_data is None: - raise SetupError(f"Component modeler has no batch saved. Run .solve() to generate.") + raise SetupError("Component modeler has no batch saved. Run .solve() to generate.") return self._construct_smatrix(batch_data=self.batch_data) From 2cdd0e37aad8b06d32e7ccf87eff3cb9fc1d8e61 Mon Sep 17 00:00:00 2001 From: Shashwat Sharma Date: Mon, 4 Apr 2022 21:41:50 -0400 Subject: [PATCH 3/8] fixed pylint warnings for near2far --- .pylintrc | 2 +- tidy3d/plugins/near2far/near2far.py | 54 +++++++++++++---------------- 2 files changed, 26 insertions(+), 30 deletions(-) diff --git a/.pylintrc b/.pylintrc index 3a29be50e9..a5a9570518 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,6 @@ [MASTER] extension-pkg-allow-list=pydantic -good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz +good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Ay, Az, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz, nk, i, j, k, J, M, _N_th, _N_ph, _L_th, _L_ph, r, Et_array, Ep_array, Er_array, Ht_array, Hp_array, Hr_array, Et, Ep, Er, Ht, Hp, Hr, Ex_data, Ey_data, Ez_data, Hx_data, Hy_data, Hz_data, Ar, Atheta, Aphi [BASIC] diff --git a/tidy3d/plugins/near2far/near2far.py b/tidy3d/plugins/near2far/near2far.py index c0b4ca843d..f746bcc016 100644 --- a/tidy3d/plugins/near2far/near2far.py +++ b/tidy3d/plugins/near2far/near2far.py @@ -49,7 +49,7 @@ def is_plane(cls, val): """Ensures that the monitor is a plane, i.e., its `size` attribute has exactly 1 zero""" size = val.size if size.count(0.0) != 1: - raise ValidationError(f"'{cls.__name__}' object must be planar, given size={val}") + raise ValidationError(f"Monitor '{val.name}' must be planar, given size={size}") return val @@ -402,7 +402,7 @@ def _radiation_vectors_for_surface( Returns ------- - tuple(numpy.ndarray[float], numpy.ndarray[float], numpy.ndarray[float], numpy.ndarray[float]) + tuple(numpy.ndarray[float],numpy.ndarray[float],numpy.ndarray[float],numpy.ndarray[float]) ``N_theta``, ``N_phi``, ``L_theta``, ``L_phi`` radiation vectors for the given surface. """ @@ -426,46 +426,46 @@ def _radiation_vectors_for_surface( J = np.zeros((3, len(theta), len(phi)), dtype=complex) M = np.zeros_like(J) - def integrate_2D(function, pts_u, pts_v): + def integrate_2d(function, pts_u, pts_v): """Trapezoidal integration in two dimensions.""" return np.trapz(np.trapz(function, pts_u, axis=0), pts_v, axis=0) phase = [None] * 3 propagation_factor = -self.phasor_positive_sign * 1j * self.k - def integrate_for_one_theta(i: int): + def integrate_for_one_theta(i_th: int): """Perform integration for a given theta angle index""" - for j in np.arange(len(phi)): + for j_ph in np.arange(len(phi)): - phase[0] = np.exp(propagation_factor * pts[0] * sin_theta[i] * cos_phi[j]) - phase[1] = np.exp(propagation_factor * pts[1] * sin_theta[i] * sin_phi[j]) - phase[2] = np.exp(propagation_factor * pts[2] * cos_theta[i]) + phase[0] = np.exp(propagation_factor * pts[0] * sin_theta[i_th] * cos_phi[j_ph]) + phase[1] = np.exp(propagation_factor * pts[1] * sin_theta[i_th] * sin_phi[j_ph]) + phase[2] = np.exp(propagation_factor * pts[2] * cos_theta[i_th]) phase_ij = phase[idx_u][:, None] * phase[idx_v][None, :] * phase[idx_w] - J[idx_u, i, j] = integrate_2D( + J[idx_u, i_th, j_ph] = integrate_2d( currents["J" + cmp_1].values * phase_ij, pts[idx_u], pts[idx_v] ) - J[idx_v, i, j] = integrate_2D( + J[idx_v, i_th, j_ph] = integrate_2d( currents["J" + cmp_2].values * phase_ij, pts[idx_u], pts[idx_v] ) - M[idx_u, i, j] = integrate_2D( + M[idx_u, i_th, j_ph] = integrate_2d( currents["M" + cmp_1].values * phase_ij, pts[idx_u], pts[idx_v] ) - M[idx_v, i, j] = integrate_2D( + M[idx_v, i_th, j_ph] = integrate_2d( currents["M" + cmp_2].values * phase_ij, pts[idx_u], pts[idx_v] ) if len(theta) < 2: integrate_for_one_theta(0) else: - for i in track( + for i_th in track( np.arange(len(theta)), description=f"Processing surface monitor '{surface.monitor.name}'...", ): - integrate_for_one_theta(i) + integrate_for_one_theta(i_th) cos_th_cos_phi = cos_theta[:, None] * cos_phi[None, :] cos_th_sin_phi = cos_theta[:, None] * sin_phi[None, :] @@ -496,7 +496,7 @@ def _radiation_vectors(self, theta: ArrayLikeN2F, phi: ArrayLikeN2F): Returns ------- - tuple(numpy.ndarray[float], numpy.ndarray[float], numpy.ndarray[float], numpy.ndarray[float]) + tuple(numpy.ndarray[float],numpy.ndarray[float],numpy.ndarray[float],numpy.ndarray[float]) ``N_theta``, ``N_phi``, ``L_theta``, ``L_phi`` radiation vectors. """ @@ -618,18 +618,14 @@ def fields_cartesian(self, x: ArrayLikeN2F, y: ArrayLikeN2F, z: ArrayLikeN2F) -> r, theta, phi = self._car_2_sph(_x, _y, _z) _field_data = self.fields_spherical(r, theta, phi) - Er, Etheta, Ephi = [ - _field_data[comp].values for comp in ["E_r", "E_theta", "E_phi"] - ] - Hr, Htheta, Hphi = [ - _field_data[comp].values for comp in ["H_r", "H_theta", "H_phi"] - ] + Er, Et, Ep = [_field_data[comp].values for comp in ["E_r", "E_theta", "E_phi"]] + Hr, Ht, Hp = [_field_data[comp].values for comp in ["H_r", "H_theta", "H_phi"]] Ex_data[i, j, k], Ey_data[i, j, k], Ez_data[i, j, k] = self._sph_2_car_field( - Er, Etheta, Ephi, theta, phi + Er, Et, Ep, theta, phi ) Hx_data[i, j, k], Hy_data[i, j, k], Hz_data[i, j, k] = self._sph_2_car_field( - Hr, Htheta, Hphi, theta, phi + Hr, Ht, Hp, theta, phi ) dims = ("x", "y", "z") @@ -669,10 +665,10 @@ def power_spherical(self, r: float, theta: ArrayLikeN2F, phi: ArrayLikeN2F) -> x phi = np.atleast_1d(phi) field_data = self.fields_spherical(r, theta, phi) - E_theta, E_phi = [field_data[comp].values for comp in ["E_theta", "E_phi"]] - H_theta, H_phi = [field_data[comp].values for comp in ["H_theta", "H_phi"]] - power_theta = 0.5 * np.real(E_theta * np.conj(H_phi)) - power_phi = 0.5 * np.real(-E_phi * np.conj(H_theta)) + Et, Ep = [field_data[comp].values for comp in ["E_theta", "E_phi"]] + Ht, Hp = [field_data[comp].values for comp in ["H_theta", "H_phi"]] + power_theta = 0.5 * np.real(Et * np.conj(Hp)) + power_phi = 0.5 * np.real(-Ep * np.conj(Ht)) power_data = power_theta + power_phi dims = ("r", "theta", "phi") @@ -751,12 +747,12 @@ def radar_cross_section(self, theta: ArrayLikeN2F, phi: ArrayLikeN2F) -> xr.Data constant = k**2 / (8 * np.pi * eta) term1 = np.abs(L_phi + eta * N_theta) ** 2 term2 = np.abs(L_theta - eta * N_phi) ** 2 - RCS_data = constant * (term1 + term2) + rcs_data = constant * (term1 + term2) dims = ("theta", "phi") coords = {"theta": theta, "phi": phi} - return xr.DataArray(data=RCS_data, coords=coords, dims=dims) + return xr.DataArray(data=rcs_data, coords=coords, dims=dims) @staticmethod def _car_2_sph(x, y, z): From fc91b5bce59037a6d1fddf3d5bc56fc67374fc62 Mon Sep 17 00:00:00 2001 From: momchil Date: Wed, 6 Apr 2022 09:50:02 -0700 Subject: [PATCH 4/8] Linting plugins/mode --- .pylintrc | 2 +- tidy3d/plugins/mode/derivatives.py | 162 +++++++++++++++-------------- tidy3d/plugins/mode/mode_solver.py | 2 + tidy3d/plugins/mode/solver.py | 149 +++++++++++++------------- tidy3d/plugins/mode/transforms.py | 2 +- 5 files changed, 162 insertions(+), 155 deletions(-) diff --git a/.pylintrc b/.pylintrc index a5a9570518..5e2ddca37c 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,6 @@ [MASTER] extension-pkg-allow-list=pydantic -good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Ay, Az, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz, nk, i, j, k, J, M, _N_th, _N_ph, _L_th, _L_ph, r, Et_array, Ep_array, Er_array, Ht_array, Hp_array, Hr_array, Et, Ep, Er, Ht, Hp, Hr, Ex_data, Ey_data, Ez_data, Hx_data, Hy_data, Hz_data, Ar, Atheta, Aphi +good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Ay, Az, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz, nk, i, j, k, J, M, _N_th, _N_ph, _L_th, _L_ph, r, Et_array, Ep_array, Er_array, Ht_array, Hp_array, Hr_array, Et, Ep, Er, Ht, Hp, Hr, Ex_data, Ey_data, Ez_data, Hx_data, Hy_data, Hz_data, Ar, Atheta, Aphi, mu [BASIC] diff --git a/tidy3d/plugins/mode/derivatives.py b/tidy3d/plugins/mode/derivatives.py index 759ea0ecdf..8d55008a9b 100644 --- a/tidy3d/plugins/mode/derivatives.py +++ b/tidy3d/plugins/mode/derivatives.py @@ -1,73 +1,76 @@ +""" Finite-difference derivatives and PML absorption operators expressed as sparse matrices. """ + import numpy as np import scipy.sparse as sp from ...constants import EPSILON_0, ETA_0 -def make_Dxf(dLs, shape): +def make_dxf(dls, shape): """Forward derivative in x.""" Nx, Ny = shape if Nx == 1: return sp.csr_matrix((Ny, Ny)) - Dxf = sp.diags([-1, 1], [0, 1], shape=(Nx, Nx)) - Dxf = sp.diags(1 / dLs).dot(Dxf) - Dxf = sp.kron(Dxf, sp.eye(Ny)) - return Dxf + dxf = sp.diags([-1, 1], [0, 1], shape=(Nx, Nx)) + dxf = sp.diags(1 / dls).dot(dxf) + dxf = sp.kron(dxf, sp.eye(Ny)) + return dxf -def make_Dxb(dLs, shape, pmc): +def make_dxb(dls, shape, pmc): """Backward derivative in x.""" Nx, Ny = shape if Nx == 1: return sp.csr_matrix((Ny, Ny)) - Dxb = sp.diags([1, -1], [0, -1], shape=(Nx, Nx)) - if pmc == True: - Dxb = sp.csr_matrix(Dxb) - Dxb[0, 0] = 2.0 - Dxb = sp.diags(1 / dLs).dot(Dxb) - Dxb = sp.kron(Dxb, sp.eye(Ny)) - return Dxb + dxb = sp.diags([1, -1], [0, -1], shape=(Nx, Nx)) + if pmc: + dxb = sp.csr_matrix(dxb) + dxb[0, 0] = 2.0 + dxb = sp.diags(1 / dls).dot(dxb) + dxb = sp.kron(dxb, sp.eye(Ny)) + return dxb -def make_Dyf(dLs, shape): +def make_dyf(dls, shape): """Forward derivative in y.""" Nx, Ny = shape if Ny == 1: return sp.csr_matrix((Nx, Nx)) - Dyf = sp.diags([-1, 1], [0, 1], shape=(Ny, Ny)) - Dyf = sp.diags(1 / dLs).dot(Dyf) - Dyf = sp.kron(sp.eye(Nx), Dyf) - return Dyf + dyf = sp.diags([-1, 1], [0, 1], shape=(Ny, Ny)) + dyf = sp.diags(1 / dls).dot(dyf) + dyf = sp.kron(sp.eye(Nx), dyf) + return dyf -def make_Dyb(dLs, shape, pmc): +def make_dyb(dls, shape, pmc): """Backward derivative in y.""" Nx, Ny = shape if Ny == 1: return sp.csr_matrix((Nx, Nx)) - Dyb = sp.diags([1, -1], [0, -1], shape=(Ny, Ny)) - if pmc == True: - Dyb = sp.csr_matrix(Dyb) - Dyb[0, 0] = 2.0 - Dyb = sp.diags(1 / dLs).dot(Dyb) - Dyb = sp.kron(sp.eye(Nx), Dyb) - return Dyb + dyb = sp.diags([1, -1], [0, -1], shape=(Ny, Ny)) + if pmc: + dyb = sp.csr_matrix(dyb) + dyb[0, 0] = 2.0 + dyb = sp.diags(1 / dls).dot(dyb) + dyb = sp.kron(sp.eye(Nx), dyb) + return dyb -def create_D_matrices(shape, dLf, dLb, dmin_pmc=(False, False)): +def create_d_matrices(shape, dlf, dlb, dmin_pmc=(False, False)): """Make the derivative matrices without PML. If dmin_pmc is True, the 'backward' derivative in that dimension will be set to implement PMC symmetry.""" - Dxf = make_Dxf(dLf[0], shape) - Dxb = make_Dxb(dLb[0], shape, dmin_pmc[0]) - Dyf = make_Dyf(dLf[1], shape) - Dyb = make_Dyb(dLb[1], shape, dmin_pmc[1]) + dxf = make_dxf(dlf[0], shape) + dxb = make_dxb(dlb[0], shape, dmin_pmc[0]) + dyf = make_dyf(dlf[1], shape) + dyb = make_dyb(dlb[1], shape, dmin_pmc[1]) - return (Dxf, Dxb, Dyf, Dyb) + return (dxf, dxb, dyf, dyb) -def create_S_matrices(omega, shape, npml, dLf, dLb, dmin_pml=(True, True)): +# pylint:disable=too-many-locals, too-many-arguments +def create_s_matrices(omega, shape, npml, dlf, dlb, dmin_pml=(True, True)): """Makes the 'S-matrices'. When dotted with derivative matrices, they add PML. If dmin_pml is set to False, PML will not be applied on the "bottom" side of the domain.""" @@ -75,88 +78,89 @@ def create_S_matrices(omega, shape, npml, dLf, dLb, dmin_pml=(True, True)): # strip out some information needed Nx, Ny = shape N = Nx * Ny - Nx_pml, Ny_pml = npml + nx_pml, ny_pml = npml # Create the sfactor in each direction and for 'f' and 'b' - s_vector_x_f = create_sfactor("f", omega, dLf[0], Nx, Nx_pml, dmin_pml[0]) - s_vector_x_b = create_sfactor("b", omega, dLb[0], Nx, Nx_pml, dmin_pml[0]) - s_vector_y_f = create_sfactor("f", omega, dLf[1], Ny, Ny_pml, dmin_pml[1]) - s_vector_y_b = create_sfactor("b", omega, dLb[1], Ny, Ny_pml, dmin_pml[1]) + s_vector_x_f = create_sfactor("f", omega, dlf[0], Nx, nx_pml, dmin_pml[0]) + s_vector_x_b = create_sfactor("b", omega, dlb[0], Nx, nx_pml, dmin_pml[0]) + s_vector_y_f = create_sfactor("f", omega, dlf[1], Ny, ny_pml, dmin_pml[1]) + s_vector_y_b = create_sfactor("b", omega, dlb[1], Ny, ny_pml, dmin_pml[1]) - # Fill the 2D space with layers of appropriate s-factors - Sx_f_2D = np.zeros(shape, dtype=np.complex128) - Sx_b_2D = np.zeros(shape, dtype=np.complex128) - Sy_f_2D = np.zeros(shape, dtype=np.complex128) - Sy_b_2D = np.zeros(shape, dtype=np.complex128) + # Fill the 2d space with layers of appropriate s-factors + sx_f_2d = np.zeros(shape, dtype=np.complex128) + sx_b_2d = np.zeros(shape, dtype=np.complex128) + sy_f_2d = np.zeros(shape, dtype=np.complex128) + sy_b_2d = np.zeros(shape, dtype=np.complex128) # Insert the cross sections into the S-grids (could be done more elegantly) for i in range(0, Ny): - Sx_f_2D[:, i] = 1 / s_vector_x_f - Sx_b_2D[:, i] = 1 / s_vector_x_b + sx_f_2d[:, i] = 1 / s_vector_x_f + sx_b_2d[:, i] = 1 / s_vector_x_b for i in range(0, Nx): - Sy_f_2D[i, :] = 1 / s_vector_y_f - Sy_b_2D[i, :] = 1 / s_vector_y_b + sy_f_2d[i, :] = 1 / s_vector_y_f + sy_b_2d[i, :] = 1 / s_vector_y_b - # Reshape the 2D s-factors into a 1D s-vecay - Sx_f_vec = Sx_f_2D.flatten() - Sx_b_vec = Sx_b_2D.flatten() - Sy_f_vec = Sy_f_2D.flatten() - Sy_b_vec = Sy_b_2D.flatten() + # Reshape the 2d s-factors into a 1D s-vecay + sx_f_vec = sx_f_2d.flatten() + sx_b_vec = sx_b_2d.flatten() + sy_f_vec = sy_f_2d.flatten() + sy_b_vec = sy_b_2d.flatten() # Construct the 1D total s-vector into a diagonal matrix - Sx_f = sp.spdiags(Sx_f_vec, 0, N, N) - Sx_b = sp.spdiags(Sx_b_vec, 0, N, N) - Sy_f = sp.spdiags(Sy_f_vec, 0, N, N) - Sy_b = sp.spdiags(Sy_b_vec, 0, N, N) + sx_f = sp.spdiags(sx_f_vec, 0, N, N) + sx_b = sp.spdiags(sx_b_vec, 0, N, N) + sy_f = sp.spdiags(sy_f_vec, 0, N, N) + sy_b = sp.spdiags(sy_b_vec, 0, N, N) - return Sx_f, Sx_b, Sy_f, Sy_b + return sx_f, sx_b, sy_f, sy_b -def create_sfactor(direction, omega, dLs, N, N_pml, dmin_pml): +# pylint:disable=too-many-arguments +def create_sfactor(direction, omega, dls, N, n_pml, dmin_pml): """Creates the S-factor cross section needed in the S-matrices""" # For no PNL, this should just be identity matrix. - if N_pml == 0: + if n_pml == 0: return np.ones(N, dtype=np.complex128) # Otherwise, get different profiles for forward and reverse derivatives. if direction == "f": - return create_sfactor_f(omega, dLs, N, N_pml, dmin_pml) - elif direction == "b": - return create_sfactor_b(omega, dLs, N, N_pml, dmin_pml) - else: - raise ValueError("Direction value {} not recognized".format(direction)) + return create_sfactor_f(omega, dls, N, n_pml, dmin_pml) + if direction == "b": + return create_sfactor_b(omega, dls, N, n_pml, dmin_pml) + + raise ValueError(f"Direction value {direction} not recognized") -def create_sfactor_f(omega, dLs, N, N_pml, dmin_pml): +def create_sfactor_f(omega, dls, N, n_pml, dmin_pml): """S-factor profile for forward derivative matrix""" sfactor_array = np.ones(N, dtype=np.complex128) for i in range(N): - if i <= N_pml and dmin_pml: - sfactor_array[i] = s_value(dLs[0], (N_pml - i + 0.5) / N_pml, omega) - elif i > N - N_pml: - sfactor_array[i] = s_value(dLs[-1], (i - (N - N_pml) - 0.5) / N_pml, omega) + if i <= n_pml and dmin_pml: + sfactor_array[i] = s_value(dls[0], (n_pml - i + 0.5) / n_pml, omega) + elif i > N - n_pml: + sfactor_array[i] = s_value(dls[-1], (i - (N - n_pml) - 0.5) / n_pml, omega) return sfactor_array -def create_sfactor_b(omega, dLs, N, N_pml, dmin_pml): +def create_sfactor_b(omega, dls, N, n_pml, dmin_pml): """S-factor profile for backward derivative matrix""" sfactor_array = np.ones(N, dtype=np.complex128) for i in range(N): - if i <= N_pml and dmin_pml: - sfactor_array[i] = s_value(dLs[0], (N_pml - i + 1) / N_pml, omega) - elif i > N - N_pml: - sfactor_array[i] = s_value(dLs[-1], (i - (N - N_pml) - 1) / N_pml, omega) + if i <= n_pml and dmin_pml: + sfactor_array[i] = s_value(dls[0], (n_pml - i + 1) / n_pml, omega) + elif i > N - n_pml: + sfactor_array[i] = s_value(dls[-1], (i - (N - n_pml) - 1) / n_pml, omega) return sfactor_array -def sig_w(dL, step, sorder=3): +def sig_w(dl, step, sorder=3): """Fictional conductivity, note that these values might need tuning""" - sig_max = 0.8 * (sorder + 1) / (ETA_0 * dL) - return sig_max * step**sorder + sig_max = 0.8 * (sorder + 1) / (ETA_0 * dl) + return sig_max * step ** sorder -def s_value(dL, step, omega): +def s_value(dl, step, omega): """S-value to use in the S-matrices""" # print(step) - return 1 - 1j * sig_w(dL, step) / (omega * EPSILON_0) + return 1 - 1j * sig_w(dl, step) / (omega * EPSILON_0) diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index 21bb8248f8..4490815019 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/tidy3d/plugins/mode/mode_solver.py @@ -280,6 +280,7 @@ def plane_sym(self): """Potentially smaller plane if symmetries present in the simulation.""" return self.simulation.min_sym_box(self.plane) + # pylint:disable=too-many-locals def solve(self) -> ModeSolverData: """Finds the modal profile and effective index of the modes. @@ -398,6 +399,7 @@ def rotate_field_coords(self, field): f_rot = np.stack(self.plane.unpop_axis(f_z, (f_x, f_y), axis=self.normal_axis), axis=0) return f_rot + # pylint:disable=too-many-locals def process_fields( self, mode_fields: Array[complex], freq_index: int, mode_index: int ) -> Tuple[FIELD, FIELD]: diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index 4a857ec36e..f8320a24cb 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -7,11 +7,12 @@ from ...components.types import Numpy from ...constants import ETA_0, C_0, fp_eps, pec_val -from .derivatives import create_D_matrices as D_mats -from .derivatives import create_S_matrices as S_mats +from .derivatives import create_d_matrices as d_mats +from .derivatives import create_s_matrices as s_mats from .transforms import radial_transform, angled_transform +# pylint:disable=too-many-statements,too-many-branches,too-many-locals def compute_modes( eps_cross, coords, @@ -136,20 +137,20 @@ def compute_modes( dmin_pmc[1] = True # Primal grid steps for E-field derivatives - dLf = [c[1:] - c[:-1] for c in new_coords] + dl_f = [new_cs[1:] - new_cs[:-1] for new_cs in new_coords] # Dual grid steps for H-field derivatives - dLtmp = [(dL[:-1] + dL[1:]) / 2 for dL in dLf] - dLb = [np.hstack((d1[0], d2)) for d1, d2 in zip(dLf, dLtmp)] + dl_tmp = [(dl[:-1] + dl[1:]) / 2 for dl in dl_f] + dl_b = [np.hstack((d1[0], d2)) for d1, d2 in zip(dl_f, dl_tmp)] # Derivative matrices with PEC boundaries at the far end and optional pmc at the near end - Dmats = D_mats((Nx, Ny), dLf, dLb, dmin_pmc) + der_mats_tmp = d_mats((Nx, Ny), dl_f, dl_b, dmin_pmc) # PML matrices; do not impose PML on the bottom when symmetry present dmin_pml = np.array(symmetry) == 0 - Smats = S_mats(omega, (Nx, Ny), mode_spec.num_pml, dLf, dLb, dmin_pml) + pml_mats = s_mats(omega, (Nx, Ny), mode_spec.num_pml, dl_f, dl_b, dmin_pml) # Add the PML on top of the derivatives; normalize by k0 to match the EM-possible notation - SDmats = [Smat.dot(Dmat) / k0 for Smat, Dmat in zip(Smats, Dmats)] + der_mats = [Smat.dot(Dmat) / k0 for Smat, Dmat in zip(pml_mats, der_mats_tmp)] # Determine initial guess value for the solver in transformed coordinates if mode_spec.target_neff is None: @@ -162,7 +163,7 @@ def compute_modes( target_neff_p = target / np.linalg.norm(kp_to_k) # Solve for the modes - E, H, neff, keff = solver_em(eps_tensor, mu_tensor, SDmats, num_modes, target_neff_p) + E, H, neff, keff = solver_em(eps_tensor, mu_tensor, der_mats, num_modes, target_neff_p) # Reorder if needed if mode_spec.sort_by != "largest_neff": @@ -178,17 +179,17 @@ def compute_modes( # Transform back to original axes, E = J^T E' E = np.sum(jac_e[..., None] * E[:, None, ...], axis=0) - E = E.reshape(3, Nx, Ny, 1, num_modes) + E = E.reshape((3, Nx, Ny, 1, num_modes)) H = np.sum(jac_h[..., None] * H[:, None, ...], axis=0) - H = H.reshape(3, Nx, Ny, 1, num_modes) + H = H.reshape((3, Nx, Ny, 1, num_modes)) neff = neff * np.linalg.norm(kp_to_k) - F = np.stack((E, H), axis=0) + fields = np.stack((E, H), axis=0) - return F, neff + 1j * keff + return fields, neff + 1j * keff -def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): +def solver_em(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess): """Solve for the electromagnetic modes of a system defined by in-plane permittivity and permeability and assuming translational invariance in the normal direction. @@ -198,8 +199,8 @@ def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): Shape (3, 3, N), the permittivity tensor at every point in the plane. mu_tensor : np.ndarray Shape (3, 3, N), the permittivity tensor at every point in the plane. - SDmats : List[scipy.sparse.csr_matrix] - The sparce derivative matrices Dxf, Dxb, Dyf, Dyb, including the PML. + der_mats : List[scipy.sparse.csr_matrix] + The sparce derivative matrices dxf, dxb, dyf, dyb, including the PML. num_modes : int Number of modes to solve for. neff_guess : float @@ -221,12 +222,12 @@ def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): eps_offd = np.abs(eps_tensor[off_diagonals]) mu_offd = np.abs(mu_tensor[off_diagonals]) if np.any(eps_offd > 1e-6) or np.any(mu_offd > 1e-6): - return solver_tensorial(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) - else: - return solver_diagonal(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) + return solver_tensorial(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess) + + return solver_diagonal(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess) -def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): +def solver_diagonal(eps, mu, der_mats, num_modes, neff_guess): """EM eigenmode solver assuming ``eps`` and ``mu`` are diagonal everywhere.""" N = eps.shape[-1] @@ -238,26 +239,26 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): mu_xx = mu[0, 0, :] mu_yy = mu[1, 1, :] mu_zz = mu[2, 2, :] - Dxf, Dxb, Dyf, Dyb = SDmats + dxf, dxb, dyf, dyb = der_mats # Compute the matrix for diagonalization inv_eps_zz = sp.spdiags(1 / eps_zz, [0], N, N) inv_mu_zz = sp.spdiags(1 / mu_zz, [0], N, N) - P11 = -Dxf.dot(inv_eps_zz).dot(Dyb) - P12 = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags(mu_yy, [0], N, N) - P21 = -Dyf.dot(inv_eps_zz).dot(Dyb) - sp.spdiags(mu_xx, [0], N, N) - P22 = Dyf.dot(inv_eps_zz).dot(Dxb) - Q11 = -Dxb.dot(inv_mu_zz).dot(Dyf) - Q12 = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags(eps_yy, [0], N, N) - Q21 = -Dyb.dot(inv_mu_zz).dot(Dyf) - sp.spdiags(eps_xx, [0], N, N) - Q22 = Dyb.dot(inv_mu_zz).dot(Dxf) - - Pmat = sp.bmat([[P11, P12], [P21, P22]]) - Qmat = sp.bmat([[Q11, Q12], [Q21, Q22]]) - A = Pmat.dot(Qmat) + p11 = -dxf.dot(inv_eps_zz).dot(dyb) + p12 = dxf.dot(inv_eps_zz).dot(dxb) + sp.spdiags(mu_yy, [0], N, N) + p21 = -dyf.dot(inv_eps_zz).dot(dyb) - sp.spdiags(mu_xx, [0], N, N) + p22 = dyf.dot(inv_eps_zz).dot(dxb) + q11 = -dxb.dot(inv_mu_zz).dot(dyf) + q12 = dxb.dot(inv_mu_zz).dot(dxf) + sp.spdiags(eps_yy, [0], N, N) + q21 = -dyb.dot(inv_mu_zz).dot(dyf) - sp.spdiags(eps_xx, [0], N, N) + q22 = dyb.dot(inv_mu_zz).dot(dxf) + + pmat = sp.bmat([[p11, p12], [p21, p22]]) + qmat = sp.bmat([[q11, q12], [q21, q22]]) + mat = pmat.dot(qmat) # Call the eigensolver. The eigenvalues are -(neff + 1j * keff)**2 - vals, vecs = solver_eigs(A, num_modes, guess_value=-(neff_guess**2)) + vals, vecs = solver_eigs(mat, num_modes, guess_value=-(neff_guess ** 2)) if vals.size == 0: raise RuntimeError("Could not find any eigenmodes for this waveguide") vre, vim = -np.real(vals), -np.imag(vals) @@ -269,7 +270,7 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): vecs = vecs[:, sort_inds] # Real and imaginary part of the effective index - neff = np.sqrt(vre / 2 + np.sqrt(vre**2 + vim**2) / 2) + neff = np.sqrt(vre / 2 + np.sqrt(vre ** 2 + vim ** 2) / 2) keff = vim / 2 / (neff + 1e-10) # Field components from eigenvectors @@ -277,11 +278,11 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): Ey = vecs[N:, :] # Get the other field components - Hs = Qmat.dot(vecs) - Hx = Hs[:N, :] / (1j * neff - keff) - Hy = Hs[N:, :] / (1j * neff - keff) - Hz = inv_mu_zz.dot((Dxf.dot(Ey) - Dyf.dot(Ex))) - Ez = inv_eps_zz.dot((Dxb.dot(Hy) - Dyb.dot(Hx))) + h_field = qmat.dot(vecs) + Hx = h_field[:N, :] / (1j * neff - keff) + Hy = h_field[N:, :] / (1j * neff - keff) + Hz = inv_mu_zz.dot((dxf.dot(Ey) - dyf.dot(Ex))) + Ez = inv_eps_zz.dot((dxb.dot(Hy) - dyb.dot(Hx))) # Bundle up E = np.stack((Ex, Ey, Ez), axis=0) @@ -293,65 +294,65 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): return E, H, neff, keff -def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): +def solver_tensorial(eps, mu, der_mats, num_modes, neff_guess): """EM eigenmode solver assuming ``eps`` or ``mu`` have off-diagonal elements.""" N = eps.shape[-1] - Dxf, Dxb, Dyf, Dyb = SDmats + dxf, dxb, dyf, dyb = der_mats # Compute all blocks of the matrix for diagonalization inv_eps_zz = sp.spdiags(1 / eps[2, 2, :], [0], N, N) inv_mu_zz = sp.spdiags(1 / mu[2, 2, :], [0], N, N) - axax = -Dxf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + axax = -dxf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( mu[1, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dyf) - axay = -Dxf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + ).dot(dyf) + axay = -dxf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( mu[1, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dxf) - axbx = -Dxf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + ).dot(dxf) + axbx = -dxf.dot(inv_eps_zz).dot(dyb) + sp.spdiags( mu[1, 0, :] - mu[1, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N ) - axby = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + axby = dxf.dot(inv_eps_zz).dot(dxb) + sp.spdiags( mu[1, 1, :] - mu[1, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N ) - ayax = -Dyf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + ayax = -dyf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( mu[0, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dyf) - ayay = -Dyf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + ).dot(dyf) + ayay = -dyf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( mu[0, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dxf) - aybx = -Dyf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + ).dot(dxf) + aybx = -dyf.dot(inv_eps_zz).dot(dyb) + sp.spdiags( -mu[0, 0, :] + mu[0, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N ) - ayby = Dyf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + ayby = dyf.dot(inv_eps_zz).dot(dxb) + sp.spdiags( -mu[0, 1, :] + mu[0, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N ) - bxbx = -Dxb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + bxbx = -dxb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( eps[1, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dyb) - bxby = -Dxb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + ).dot(dyb) + bxby = -dxb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( eps[1, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dxb) - bxax = -Dxb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + ).dot(dxb) + bxax = -dxb.dot(inv_mu_zz).dot(dyf) + sp.spdiags( eps[1, 0, :] - eps[1, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N ) - bxay = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + bxay = dxb.dot(inv_mu_zz).dot(dxf) + sp.spdiags( eps[1, 1, :] - eps[1, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N ) - bybx = -Dyb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + bybx = -dyb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( eps[0, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dyb) - byby = -Dyb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + ).dot(dyb) + byby = -dyb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( eps[0, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dxb) - byax = -Dyb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + ).dot(dxb) + byax = -dyb.dot(inv_mu_zz).dot(dyf) + sp.spdiags( -eps[0, 0, :] + eps[0, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N ) - byay = Dyb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + byay = dyb.dot(inv_mu_zz).dot(dxf) + sp.spdiags( -eps[0, 1, :] + eps[0, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N ) - A = sp.bmat( + mat = sp.bmat( [ [axax, axay, axbx, axby], [ayax, ayay, aybx, ayby], @@ -361,7 +362,7 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): ) # Call the eigensolver. The eigenvalues are 1j * (neff + 1j * keff) - vals, vecs = solver_eigs(A, num_modes, guess_value=1j * neff_guess) + vals, vecs = solver_eigs(mat, num_modes, guess_value=1j * neff_guess) if vals.size == 0: raise RuntimeError("Could not find any eigenmodes for this waveguide") # Real and imaginary part of the effective index @@ -381,9 +382,9 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): # Get the other field components hxy_term = (-mu[2, 0, :] * Hx.T - mu[2, 1, :] * Hy.T).T - Hz = inv_mu_zz.dot(Dxf.dot(Ey) - Dyf.dot(Ex) + hxy_term) + Hz = inv_mu_zz.dot(dxf.dot(Ey) - dyf.dot(Ex) + hxy_term) exy_term = (-eps[2, 0, :] * Ex.T - eps[2, 1, :] * Ey.T).T - Ez = inv_eps_zz.dot(Dxb.dot(Hy) - Dyb.dot(Hx) + exy_term) + Ez = inv_eps_zz.dot(dxb.dot(Hy) - dyb.dot(Hx) + exy_term) # Bundle up E = np.stack((Ex, Ey, Ez), axis=0) @@ -396,17 +397,17 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): return E, H, neff, keff -def solver_eigs(A, num_modes, guess_value=1.0): - """Find ``num_modes`` eigenmodes of ``A`` cloest to ``guess_value``. +def solver_eigs(mat, num_modes, guess_value=1.0): + """Find ``num_modes`` eigenmodes of ``mat`` cloest to ``guess_value``. Parameters ---------- - A : scipy.sparse matrix + mat : scipy.sparse matrix Square matrix for diagonalization. num_modes : int Number of eigenmodes to compute. guess_value : float, optional """ - values, vectors = spl.eigs(A, k=num_modes, sigma=guess_value, tol=fp_eps) + values, vectors = spl.eigs(mat, k=num_modes, sigma=guess_value, tol=fp_eps) return values, vectors diff --git a/tidy3d/plugins/mode/transforms.py b/tidy3d/plugins/mode/transforms.py index c3403e7c6e..44a0b752ee 100644 --- a/tidy3d/plugins/mode/transforms.py +++ b/tidy3d/plugins/mode/transforms.py @@ -53,7 +53,7 @@ def radial_transform(coords, radius, bend_axis): new_coords = (u, v) """The only nontrivial derivative is dwdz and it only depends on the coordinate in the - norm_axis direction (orthogonal to both bend_axis and z). We need to compute that derivative + norm_axis direction (orthogonal to both bend_axis and z). We need to compute that derivative at the En and Hn positions. """ dwdz_e = radius / new_coords[norm_axis][:-1] From d1e1aaa371b82242a718b1d3f0a28cb894a030f3 Mon Sep 17 00:00:00 2001 From: tylerflex Date: Wed, 6 Apr 2022 09:55:24 -0700 Subject: [PATCH 5/8] added weiliangs pylint in --- tidy3d/plugins/dispersion/fit_web.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tidy3d/plugins/dispersion/fit_web.py b/tidy3d/plugins/dispersion/fit_web.py index 10d85dc507..1f323d83cb 100644 --- a/tidy3d/plugins/dispersion/fit_web.py +++ b/tidy3d/plugins/dispersion/fit_web.py @@ -48,7 +48,7 @@ class AdvancedFitterParam(BaseModel): description="Stability constraint: 'hard' constraints are generally recommended since " "they are faster to compute per iteration, and they often require fewer iterations to " "converge since the search space is smaller. But sometimes the search space is " - "so restrictive that all good solutions are missed, then please try the 'soft' constraints " + "so restrictive that all good solutions are missed, then please try the 'soft' constraints " "for larger search space. However, both constraints improve stability equally well.", ) nlopt_maxeval: PositiveInt = Field( @@ -175,7 +175,7 @@ def _setup_server(url_server: str): return headers - def fit( # pylint:disable=arguments-differ + def fit( # pylint:disable=arguments-differ, too-many-locals self, num_poles: PositiveInt = 1, num_tries: PositiveInt = 50, From b0db02eec46fed2bffea1db8cef858fba7447e52 Mon Sep 17 00:00:00 2001 From: tylerflex Date: Wed, 6 Apr 2022 09:56:17 -0700 Subject: [PATCH 6/8] looks good --- tidy3d/plugins/mode/derivatives.py | 2 +- tidy3d/plugins/mode/solver.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tidy3d/plugins/mode/derivatives.py b/tidy3d/plugins/mode/derivatives.py index 8d55008a9b..37019323a9 100644 --- a/tidy3d/plugins/mode/derivatives.py +++ b/tidy3d/plugins/mode/derivatives.py @@ -157,7 +157,7 @@ def create_sfactor_b(omega, dls, N, n_pml, dmin_pml): def sig_w(dl, step, sorder=3): """Fictional conductivity, note that these values might need tuning""" sig_max = 0.8 * (sorder + 1) / (ETA_0 * dl) - return sig_max * step ** sorder + return sig_max * step**sorder def s_value(dl, step, omega): diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index f8320a24cb..68f13febc2 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -258,7 +258,7 @@ def solver_diagonal(eps, mu, der_mats, num_modes, neff_guess): mat = pmat.dot(qmat) # Call the eigensolver. The eigenvalues are -(neff + 1j * keff)**2 - vals, vecs = solver_eigs(mat, num_modes, guess_value=-(neff_guess ** 2)) + vals, vecs = solver_eigs(mat, num_modes, guess_value=-(neff_guess**2)) if vals.size == 0: raise RuntimeError("Could not find any eigenmodes for this waveguide") vre, vim = -np.real(vals), -np.imag(vals) @@ -270,7 +270,7 @@ def solver_diagonal(eps, mu, der_mats, num_modes, neff_guess): vecs = vecs[:, sort_inds] # Real and imaginary part of the effective index - neff = np.sqrt(vre / 2 + np.sqrt(vre ** 2 + vim ** 2) / 2) + neff = np.sqrt(vre / 2 + np.sqrt(vre**2 + vim**2) / 2) keff = vim / 2 / (neff + 1e-10) # Field components from eigenvectors From 6b3b3be2356b0b5b3d6ca5c4c33dd32a11553b6a Mon Sep 17 00:00:00 2001 From: momchil Date: Wed, 6 Apr 2022 09:50:02 -0700 Subject: [PATCH 7/8] Linting plugins/mode --- .pylintrc | 2 +- tidy3d/plugins/mode/derivatives.py | 160 +++++++++++++++-------------- tidy3d/plugins/mode/mode_solver.py | 2 + tidy3d/plugins/mode/solver.py | 147 +++++++++++++------------- tidy3d/plugins/mode/transforms.py | 2 +- 5 files changed, 160 insertions(+), 153 deletions(-) diff --git a/.pylintrc b/.pylintrc index a5a9570518..5e2ddca37c 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,6 +1,6 @@ [MASTER] extension-pkg-allow-list=pydantic -good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Ay, Az, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz, nk, i, j, k, J, M, _N_th, _N_ph, _L_th, _L_ph, r, Et_array, Ep_array, Er_array, Ht_array, Hp_array, Hr_array, Et, Ep, Er, Ht, Hp, Hr, Ex_data, Ey_data, Ez_data, Hx_data, Hy_data, Hz_data, Ar, Atheta, Aphi +good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, x3, x4, y3, y4, xs, ys, zs, Ax, Ay, Az, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi, ux, uy, uz, nk, i, j, k, J, M, _N_th, _N_ph, _L_th, _L_ph, r, Et_array, Ep_array, Er_array, Ht_array, Hp_array, Hr_array, Et, Ep, Er, Ht, Hp, Hr, Ex_data, Ey_data, Ez_data, Hx_data, Hy_data, Hz_data, Ar, Atheta, Aphi, mu [BASIC] diff --git a/tidy3d/plugins/mode/derivatives.py b/tidy3d/plugins/mode/derivatives.py index 759ea0ecdf..37019323a9 100644 --- a/tidy3d/plugins/mode/derivatives.py +++ b/tidy3d/plugins/mode/derivatives.py @@ -1,73 +1,76 @@ +""" Finite-difference derivatives and PML absorption operators expressed as sparse matrices. """ + import numpy as np import scipy.sparse as sp from ...constants import EPSILON_0, ETA_0 -def make_Dxf(dLs, shape): +def make_dxf(dls, shape): """Forward derivative in x.""" Nx, Ny = shape if Nx == 1: return sp.csr_matrix((Ny, Ny)) - Dxf = sp.diags([-1, 1], [0, 1], shape=(Nx, Nx)) - Dxf = sp.diags(1 / dLs).dot(Dxf) - Dxf = sp.kron(Dxf, sp.eye(Ny)) - return Dxf + dxf = sp.diags([-1, 1], [0, 1], shape=(Nx, Nx)) + dxf = sp.diags(1 / dls).dot(dxf) + dxf = sp.kron(dxf, sp.eye(Ny)) + return dxf -def make_Dxb(dLs, shape, pmc): +def make_dxb(dls, shape, pmc): """Backward derivative in x.""" Nx, Ny = shape if Nx == 1: return sp.csr_matrix((Ny, Ny)) - Dxb = sp.diags([1, -1], [0, -1], shape=(Nx, Nx)) - if pmc == True: - Dxb = sp.csr_matrix(Dxb) - Dxb[0, 0] = 2.0 - Dxb = sp.diags(1 / dLs).dot(Dxb) - Dxb = sp.kron(Dxb, sp.eye(Ny)) - return Dxb + dxb = sp.diags([1, -1], [0, -1], shape=(Nx, Nx)) + if pmc: + dxb = sp.csr_matrix(dxb) + dxb[0, 0] = 2.0 + dxb = sp.diags(1 / dls).dot(dxb) + dxb = sp.kron(dxb, sp.eye(Ny)) + return dxb -def make_Dyf(dLs, shape): +def make_dyf(dls, shape): """Forward derivative in y.""" Nx, Ny = shape if Ny == 1: return sp.csr_matrix((Nx, Nx)) - Dyf = sp.diags([-1, 1], [0, 1], shape=(Ny, Ny)) - Dyf = sp.diags(1 / dLs).dot(Dyf) - Dyf = sp.kron(sp.eye(Nx), Dyf) - return Dyf + dyf = sp.diags([-1, 1], [0, 1], shape=(Ny, Ny)) + dyf = sp.diags(1 / dls).dot(dyf) + dyf = sp.kron(sp.eye(Nx), dyf) + return dyf -def make_Dyb(dLs, shape, pmc): +def make_dyb(dls, shape, pmc): """Backward derivative in y.""" Nx, Ny = shape if Ny == 1: return sp.csr_matrix((Nx, Nx)) - Dyb = sp.diags([1, -1], [0, -1], shape=(Ny, Ny)) - if pmc == True: - Dyb = sp.csr_matrix(Dyb) - Dyb[0, 0] = 2.0 - Dyb = sp.diags(1 / dLs).dot(Dyb) - Dyb = sp.kron(sp.eye(Nx), Dyb) - return Dyb + dyb = sp.diags([1, -1], [0, -1], shape=(Ny, Ny)) + if pmc: + dyb = sp.csr_matrix(dyb) + dyb[0, 0] = 2.0 + dyb = sp.diags(1 / dls).dot(dyb) + dyb = sp.kron(sp.eye(Nx), dyb) + return dyb -def create_D_matrices(shape, dLf, dLb, dmin_pmc=(False, False)): +def create_d_matrices(shape, dlf, dlb, dmin_pmc=(False, False)): """Make the derivative matrices without PML. If dmin_pmc is True, the 'backward' derivative in that dimension will be set to implement PMC symmetry.""" - Dxf = make_Dxf(dLf[0], shape) - Dxb = make_Dxb(dLb[0], shape, dmin_pmc[0]) - Dyf = make_Dyf(dLf[1], shape) - Dyb = make_Dyb(dLb[1], shape, dmin_pmc[1]) + dxf = make_dxf(dlf[0], shape) + dxb = make_dxb(dlb[0], shape, dmin_pmc[0]) + dyf = make_dyf(dlf[1], shape) + dyb = make_dyb(dlb[1], shape, dmin_pmc[1]) - return (Dxf, Dxb, Dyf, Dyb) + return (dxf, dxb, dyf, dyb) -def create_S_matrices(omega, shape, npml, dLf, dLb, dmin_pml=(True, True)): +# pylint:disable=too-many-locals, too-many-arguments +def create_s_matrices(omega, shape, npml, dlf, dlb, dmin_pml=(True, True)): """Makes the 'S-matrices'. When dotted with derivative matrices, they add PML. If dmin_pml is set to False, PML will not be applied on the "bottom" side of the domain.""" @@ -75,88 +78,89 @@ def create_S_matrices(omega, shape, npml, dLf, dLb, dmin_pml=(True, True)): # strip out some information needed Nx, Ny = shape N = Nx * Ny - Nx_pml, Ny_pml = npml + nx_pml, ny_pml = npml # Create the sfactor in each direction and for 'f' and 'b' - s_vector_x_f = create_sfactor("f", omega, dLf[0], Nx, Nx_pml, dmin_pml[0]) - s_vector_x_b = create_sfactor("b", omega, dLb[0], Nx, Nx_pml, dmin_pml[0]) - s_vector_y_f = create_sfactor("f", omega, dLf[1], Ny, Ny_pml, dmin_pml[1]) - s_vector_y_b = create_sfactor("b", omega, dLb[1], Ny, Ny_pml, dmin_pml[1]) + s_vector_x_f = create_sfactor("f", omega, dlf[0], Nx, nx_pml, dmin_pml[0]) + s_vector_x_b = create_sfactor("b", omega, dlb[0], Nx, nx_pml, dmin_pml[0]) + s_vector_y_f = create_sfactor("f", omega, dlf[1], Ny, ny_pml, dmin_pml[1]) + s_vector_y_b = create_sfactor("b", omega, dlb[1], Ny, ny_pml, dmin_pml[1]) - # Fill the 2D space with layers of appropriate s-factors - Sx_f_2D = np.zeros(shape, dtype=np.complex128) - Sx_b_2D = np.zeros(shape, dtype=np.complex128) - Sy_f_2D = np.zeros(shape, dtype=np.complex128) - Sy_b_2D = np.zeros(shape, dtype=np.complex128) + # Fill the 2d space with layers of appropriate s-factors + sx_f_2d = np.zeros(shape, dtype=np.complex128) + sx_b_2d = np.zeros(shape, dtype=np.complex128) + sy_f_2d = np.zeros(shape, dtype=np.complex128) + sy_b_2d = np.zeros(shape, dtype=np.complex128) # Insert the cross sections into the S-grids (could be done more elegantly) for i in range(0, Ny): - Sx_f_2D[:, i] = 1 / s_vector_x_f - Sx_b_2D[:, i] = 1 / s_vector_x_b + sx_f_2d[:, i] = 1 / s_vector_x_f + sx_b_2d[:, i] = 1 / s_vector_x_b for i in range(0, Nx): - Sy_f_2D[i, :] = 1 / s_vector_y_f - Sy_b_2D[i, :] = 1 / s_vector_y_b + sy_f_2d[i, :] = 1 / s_vector_y_f + sy_b_2d[i, :] = 1 / s_vector_y_b - # Reshape the 2D s-factors into a 1D s-vecay - Sx_f_vec = Sx_f_2D.flatten() - Sx_b_vec = Sx_b_2D.flatten() - Sy_f_vec = Sy_f_2D.flatten() - Sy_b_vec = Sy_b_2D.flatten() + # Reshape the 2d s-factors into a 1D s-vecay + sx_f_vec = sx_f_2d.flatten() + sx_b_vec = sx_b_2d.flatten() + sy_f_vec = sy_f_2d.flatten() + sy_b_vec = sy_b_2d.flatten() # Construct the 1D total s-vector into a diagonal matrix - Sx_f = sp.spdiags(Sx_f_vec, 0, N, N) - Sx_b = sp.spdiags(Sx_b_vec, 0, N, N) - Sy_f = sp.spdiags(Sy_f_vec, 0, N, N) - Sy_b = sp.spdiags(Sy_b_vec, 0, N, N) + sx_f = sp.spdiags(sx_f_vec, 0, N, N) + sx_b = sp.spdiags(sx_b_vec, 0, N, N) + sy_f = sp.spdiags(sy_f_vec, 0, N, N) + sy_b = sp.spdiags(sy_b_vec, 0, N, N) - return Sx_f, Sx_b, Sy_f, Sy_b + return sx_f, sx_b, sy_f, sy_b -def create_sfactor(direction, omega, dLs, N, N_pml, dmin_pml): +# pylint:disable=too-many-arguments +def create_sfactor(direction, omega, dls, N, n_pml, dmin_pml): """Creates the S-factor cross section needed in the S-matrices""" # For no PNL, this should just be identity matrix. - if N_pml == 0: + if n_pml == 0: return np.ones(N, dtype=np.complex128) # Otherwise, get different profiles for forward and reverse derivatives. if direction == "f": - return create_sfactor_f(omega, dLs, N, N_pml, dmin_pml) - elif direction == "b": - return create_sfactor_b(omega, dLs, N, N_pml, dmin_pml) - else: - raise ValueError("Direction value {} not recognized".format(direction)) + return create_sfactor_f(omega, dls, N, n_pml, dmin_pml) + if direction == "b": + return create_sfactor_b(omega, dls, N, n_pml, dmin_pml) + + raise ValueError(f"Direction value {direction} not recognized") -def create_sfactor_f(omega, dLs, N, N_pml, dmin_pml): +def create_sfactor_f(omega, dls, N, n_pml, dmin_pml): """S-factor profile for forward derivative matrix""" sfactor_array = np.ones(N, dtype=np.complex128) for i in range(N): - if i <= N_pml and dmin_pml: - sfactor_array[i] = s_value(dLs[0], (N_pml - i + 0.5) / N_pml, omega) - elif i > N - N_pml: - sfactor_array[i] = s_value(dLs[-1], (i - (N - N_pml) - 0.5) / N_pml, omega) + if i <= n_pml and dmin_pml: + sfactor_array[i] = s_value(dls[0], (n_pml - i + 0.5) / n_pml, omega) + elif i > N - n_pml: + sfactor_array[i] = s_value(dls[-1], (i - (N - n_pml) - 0.5) / n_pml, omega) return sfactor_array -def create_sfactor_b(omega, dLs, N, N_pml, dmin_pml): +def create_sfactor_b(omega, dls, N, n_pml, dmin_pml): """S-factor profile for backward derivative matrix""" sfactor_array = np.ones(N, dtype=np.complex128) for i in range(N): - if i <= N_pml and dmin_pml: - sfactor_array[i] = s_value(dLs[0], (N_pml - i + 1) / N_pml, omega) - elif i > N - N_pml: - sfactor_array[i] = s_value(dLs[-1], (i - (N - N_pml) - 1) / N_pml, omega) + if i <= n_pml and dmin_pml: + sfactor_array[i] = s_value(dls[0], (n_pml - i + 1) / n_pml, omega) + elif i > N - n_pml: + sfactor_array[i] = s_value(dls[-1], (i - (N - n_pml) - 1) / n_pml, omega) return sfactor_array -def sig_w(dL, step, sorder=3): +def sig_w(dl, step, sorder=3): """Fictional conductivity, note that these values might need tuning""" - sig_max = 0.8 * (sorder + 1) / (ETA_0 * dL) + sig_max = 0.8 * (sorder + 1) / (ETA_0 * dl) return sig_max * step**sorder -def s_value(dL, step, omega): +def s_value(dl, step, omega): """S-value to use in the S-matrices""" # print(step) - return 1 - 1j * sig_w(dL, step) / (omega * EPSILON_0) + return 1 - 1j * sig_w(dl, step) / (omega * EPSILON_0) diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index 21bb8248f8..4490815019 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/tidy3d/plugins/mode/mode_solver.py @@ -280,6 +280,7 @@ def plane_sym(self): """Potentially smaller plane if symmetries present in the simulation.""" return self.simulation.min_sym_box(self.plane) + # pylint:disable=too-many-locals def solve(self) -> ModeSolverData: """Finds the modal profile and effective index of the modes. @@ -398,6 +399,7 @@ def rotate_field_coords(self, field): f_rot = np.stack(self.plane.unpop_axis(f_z, (f_x, f_y), axis=self.normal_axis), axis=0) return f_rot + # pylint:disable=too-many-locals def process_fields( self, mode_fields: Array[complex], freq_index: int, mode_index: int ) -> Tuple[FIELD, FIELD]: diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index 4a857ec36e..68f13febc2 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -7,11 +7,12 @@ from ...components.types import Numpy from ...constants import ETA_0, C_0, fp_eps, pec_val -from .derivatives import create_D_matrices as D_mats -from .derivatives import create_S_matrices as S_mats +from .derivatives import create_d_matrices as d_mats +from .derivatives import create_s_matrices as s_mats from .transforms import radial_transform, angled_transform +# pylint:disable=too-many-statements,too-many-branches,too-many-locals def compute_modes( eps_cross, coords, @@ -136,20 +137,20 @@ def compute_modes( dmin_pmc[1] = True # Primal grid steps for E-field derivatives - dLf = [c[1:] - c[:-1] for c in new_coords] + dl_f = [new_cs[1:] - new_cs[:-1] for new_cs in new_coords] # Dual grid steps for H-field derivatives - dLtmp = [(dL[:-1] + dL[1:]) / 2 for dL in dLf] - dLb = [np.hstack((d1[0], d2)) for d1, d2 in zip(dLf, dLtmp)] + dl_tmp = [(dl[:-1] + dl[1:]) / 2 for dl in dl_f] + dl_b = [np.hstack((d1[0], d2)) for d1, d2 in zip(dl_f, dl_tmp)] # Derivative matrices with PEC boundaries at the far end and optional pmc at the near end - Dmats = D_mats((Nx, Ny), dLf, dLb, dmin_pmc) + der_mats_tmp = d_mats((Nx, Ny), dl_f, dl_b, dmin_pmc) # PML matrices; do not impose PML on the bottom when symmetry present dmin_pml = np.array(symmetry) == 0 - Smats = S_mats(omega, (Nx, Ny), mode_spec.num_pml, dLf, dLb, dmin_pml) + pml_mats = s_mats(omega, (Nx, Ny), mode_spec.num_pml, dl_f, dl_b, dmin_pml) # Add the PML on top of the derivatives; normalize by k0 to match the EM-possible notation - SDmats = [Smat.dot(Dmat) / k0 for Smat, Dmat in zip(Smats, Dmats)] + der_mats = [Smat.dot(Dmat) / k0 for Smat, Dmat in zip(pml_mats, der_mats_tmp)] # Determine initial guess value for the solver in transformed coordinates if mode_spec.target_neff is None: @@ -162,7 +163,7 @@ def compute_modes( target_neff_p = target / np.linalg.norm(kp_to_k) # Solve for the modes - E, H, neff, keff = solver_em(eps_tensor, mu_tensor, SDmats, num_modes, target_neff_p) + E, H, neff, keff = solver_em(eps_tensor, mu_tensor, der_mats, num_modes, target_neff_p) # Reorder if needed if mode_spec.sort_by != "largest_neff": @@ -178,17 +179,17 @@ def compute_modes( # Transform back to original axes, E = J^T E' E = np.sum(jac_e[..., None] * E[:, None, ...], axis=0) - E = E.reshape(3, Nx, Ny, 1, num_modes) + E = E.reshape((3, Nx, Ny, 1, num_modes)) H = np.sum(jac_h[..., None] * H[:, None, ...], axis=0) - H = H.reshape(3, Nx, Ny, 1, num_modes) + H = H.reshape((3, Nx, Ny, 1, num_modes)) neff = neff * np.linalg.norm(kp_to_k) - F = np.stack((E, H), axis=0) + fields = np.stack((E, H), axis=0) - return F, neff + 1j * keff + return fields, neff + 1j * keff -def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): +def solver_em(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess): """Solve for the electromagnetic modes of a system defined by in-plane permittivity and permeability and assuming translational invariance in the normal direction. @@ -198,8 +199,8 @@ def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): Shape (3, 3, N), the permittivity tensor at every point in the plane. mu_tensor : np.ndarray Shape (3, 3, N), the permittivity tensor at every point in the plane. - SDmats : List[scipy.sparse.csr_matrix] - The sparce derivative matrices Dxf, Dxb, Dyf, Dyb, including the PML. + der_mats : List[scipy.sparse.csr_matrix] + The sparce derivative matrices dxf, dxb, dyf, dyb, including the PML. num_modes : int Number of modes to solve for. neff_guess : float @@ -221,12 +222,12 @@ def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): eps_offd = np.abs(eps_tensor[off_diagonals]) mu_offd = np.abs(mu_tensor[off_diagonals]) if np.any(eps_offd > 1e-6) or np.any(mu_offd > 1e-6): - return solver_tensorial(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) - else: - return solver_diagonal(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) + return solver_tensorial(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess) + + return solver_diagonal(eps_tensor, mu_tensor, der_mats, num_modes, neff_guess) -def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): +def solver_diagonal(eps, mu, der_mats, num_modes, neff_guess): """EM eigenmode solver assuming ``eps`` and ``mu`` are diagonal everywhere.""" N = eps.shape[-1] @@ -238,26 +239,26 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): mu_xx = mu[0, 0, :] mu_yy = mu[1, 1, :] mu_zz = mu[2, 2, :] - Dxf, Dxb, Dyf, Dyb = SDmats + dxf, dxb, dyf, dyb = der_mats # Compute the matrix for diagonalization inv_eps_zz = sp.spdiags(1 / eps_zz, [0], N, N) inv_mu_zz = sp.spdiags(1 / mu_zz, [0], N, N) - P11 = -Dxf.dot(inv_eps_zz).dot(Dyb) - P12 = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags(mu_yy, [0], N, N) - P21 = -Dyf.dot(inv_eps_zz).dot(Dyb) - sp.spdiags(mu_xx, [0], N, N) - P22 = Dyf.dot(inv_eps_zz).dot(Dxb) - Q11 = -Dxb.dot(inv_mu_zz).dot(Dyf) - Q12 = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags(eps_yy, [0], N, N) - Q21 = -Dyb.dot(inv_mu_zz).dot(Dyf) - sp.spdiags(eps_xx, [0], N, N) - Q22 = Dyb.dot(inv_mu_zz).dot(Dxf) - - Pmat = sp.bmat([[P11, P12], [P21, P22]]) - Qmat = sp.bmat([[Q11, Q12], [Q21, Q22]]) - A = Pmat.dot(Qmat) + p11 = -dxf.dot(inv_eps_zz).dot(dyb) + p12 = dxf.dot(inv_eps_zz).dot(dxb) + sp.spdiags(mu_yy, [0], N, N) + p21 = -dyf.dot(inv_eps_zz).dot(dyb) - sp.spdiags(mu_xx, [0], N, N) + p22 = dyf.dot(inv_eps_zz).dot(dxb) + q11 = -dxb.dot(inv_mu_zz).dot(dyf) + q12 = dxb.dot(inv_mu_zz).dot(dxf) + sp.spdiags(eps_yy, [0], N, N) + q21 = -dyb.dot(inv_mu_zz).dot(dyf) - sp.spdiags(eps_xx, [0], N, N) + q22 = dyb.dot(inv_mu_zz).dot(dxf) + + pmat = sp.bmat([[p11, p12], [p21, p22]]) + qmat = sp.bmat([[q11, q12], [q21, q22]]) + mat = pmat.dot(qmat) # Call the eigensolver. The eigenvalues are -(neff + 1j * keff)**2 - vals, vecs = solver_eigs(A, num_modes, guess_value=-(neff_guess**2)) + vals, vecs = solver_eigs(mat, num_modes, guess_value=-(neff_guess**2)) if vals.size == 0: raise RuntimeError("Could not find any eigenmodes for this waveguide") vre, vim = -np.real(vals), -np.imag(vals) @@ -277,11 +278,11 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): Ey = vecs[N:, :] # Get the other field components - Hs = Qmat.dot(vecs) - Hx = Hs[:N, :] / (1j * neff - keff) - Hy = Hs[N:, :] / (1j * neff - keff) - Hz = inv_mu_zz.dot((Dxf.dot(Ey) - Dyf.dot(Ex))) - Ez = inv_eps_zz.dot((Dxb.dot(Hy) - Dyb.dot(Hx))) + h_field = qmat.dot(vecs) + Hx = h_field[:N, :] / (1j * neff - keff) + Hy = h_field[N:, :] / (1j * neff - keff) + Hz = inv_mu_zz.dot((dxf.dot(Ey) - dyf.dot(Ex))) + Ez = inv_eps_zz.dot((dxb.dot(Hy) - dyb.dot(Hx))) # Bundle up E = np.stack((Ex, Ey, Ez), axis=0) @@ -293,65 +294,65 @@ def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): return E, H, neff, keff -def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): +def solver_tensorial(eps, mu, der_mats, num_modes, neff_guess): """EM eigenmode solver assuming ``eps`` or ``mu`` have off-diagonal elements.""" N = eps.shape[-1] - Dxf, Dxb, Dyf, Dyb = SDmats + dxf, dxb, dyf, dyb = der_mats # Compute all blocks of the matrix for diagonalization inv_eps_zz = sp.spdiags(1 / eps[2, 2, :], [0], N, N) inv_mu_zz = sp.spdiags(1 / mu[2, 2, :], [0], N, N) - axax = -Dxf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + axax = -dxf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( mu[1, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dyf) - axay = -Dxf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + ).dot(dyf) + axay = -dxf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( mu[1, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dxf) - axbx = -Dxf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + ).dot(dxf) + axbx = -dxf.dot(inv_eps_zz).dot(dyb) + sp.spdiags( mu[1, 0, :] - mu[1, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N ) - axby = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + axby = dxf.dot(inv_eps_zz).dot(dxb) + sp.spdiags( mu[1, 1, :] - mu[1, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N ) - ayax = -Dyf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + ayax = -dyf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( mu[0, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dyf) - ayay = -Dyf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + ).dot(dyf) + ayay = -dyf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( mu[0, 2, :] / mu[2, 2, :], [0], N, N - ).dot(Dxf) - aybx = -Dyf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + ).dot(dxf) + aybx = -dyf.dot(inv_eps_zz).dot(dyb) + sp.spdiags( -mu[0, 0, :] + mu[0, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N ) - ayby = Dyf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + ayby = dyf.dot(inv_eps_zz).dot(dxb) + sp.spdiags( -mu[0, 1, :] + mu[0, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N ) - bxbx = -Dxb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + bxbx = -dxb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( eps[1, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dyb) - bxby = -Dxb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + ).dot(dyb) + bxby = -dxb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( eps[1, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dxb) - bxax = -Dxb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + ).dot(dxb) + bxax = -dxb.dot(inv_mu_zz).dot(dyf) + sp.spdiags( eps[1, 0, :] - eps[1, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N ) - bxay = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + bxay = dxb.dot(inv_mu_zz).dot(dxf) + sp.spdiags( eps[1, 1, :] - eps[1, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N ) - bybx = -Dyb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + bybx = -dyb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( eps[0, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dyb) - byby = -Dyb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + ).dot(dyb) + byby = -dyb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( eps[0, 2, :] / eps[2, 2, :], [0], N, N - ).dot(Dxb) - byax = -Dyb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + ).dot(dxb) + byax = -dyb.dot(inv_mu_zz).dot(dyf) + sp.spdiags( -eps[0, 0, :] + eps[0, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N ) - byay = Dyb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + byay = dyb.dot(inv_mu_zz).dot(dxf) + sp.spdiags( -eps[0, 1, :] + eps[0, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N ) - A = sp.bmat( + mat = sp.bmat( [ [axax, axay, axbx, axby], [ayax, ayay, aybx, ayby], @@ -361,7 +362,7 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): ) # Call the eigensolver. The eigenvalues are 1j * (neff + 1j * keff) - vals, vecs = solver_eigs(A, num_modes, guess_value=1j * neff_guess) + vals, vecs = solver_eigs(mat, num_modes, guess_value=1j * neff_guess) if vals.size == 0: raise RuntimeError("Could not find any eigenmodes for this waveguide") # Real and imaginary part of the effective index @@ -381,9 +382,9 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): # Get the other field components hxy_term = (-mu[2, 0, :] * Hx.T - mu[2, 1, :] * Hy.T).T - Hz = inv_mu_zz.dot(Dxf.dot(Ey) - Dyf.dot(Ex) + hxy_term) + Hz = inv_mu_zz.dot(dxf.dot(Ey) - dyf.dot(Ex) + hxy_term) exy_term = (-eps[2, 0, :] * Ex.T - eps[2, 1, :] * Ey.T).T - Ez = inv_eps_zz.dot(Dxb.dot(Hy) - Dyb.dot(Hx) + exy_term) + Ez = inv_eps_zz.dot(dxb.dot(Hy) - dyb.dot(Hx) + exy_term) # Bundle up E = np.stack((Ex, Ey, Ez), axis=0) @@ -396,17 +397,17 @@ def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): return E, H, neff, keff -def solver_eigs(A, num_modes, guess_value=1.0): - """Find ``num_modes`` eigenmodes of ``A`` cloest to ``guess_value``. +def solver_eigs(mat, num_modes, guess_value=1.0): + """Find ``num_modes`` eigenmodes of ``mat`` cloest to ``guess_value``. Parameters ---------- - A : scipy.sparse matrix + mat : scipy.sparse matrix Square matrix for diagonalization. num_modes : int Number of eigenmodes to compute. guess_value : float, optional """ - values, vectors = spl.eigs(A, k=num_modes, sigma=guess_value, tol=fp_eps) + values, vectors = spl.eigs(mat, k=num_modes, sigma=guess_value, tol=fp_eps) return values, vectors diff --git a/tidy3d/plugins/mode/transforms.py b/tidy3d/plugins/mode/transforms.py index c3403e7c6e..44a0b752ee 100644 --- a/tidy3d/plugins/mode/transforms.py +++ b/tidy3d/plugins/mode/transforms.py @@ -53,7 +53,7 @@ def radial_transform(coords, radius, bend_axis): new_coords = (u, v) """The only nontrivial derivative is dwdz and it only depends on the coordinate in the - norm_axis direction (orthogonal to both bend_axis and z). We need to compute that derivative + norm_axis direction (orthogonal to both bend_axis and z). We need to compute that derivative at the En and Hn positions. """ dwdz_e = radius / new_coords[norm_axis][:-1] From 97c352c3df4ed413097fd60a2092bf6cf423a142 Mon Sep 17 00:00:00 2001 From: tylerflex Date: Wed, 6 Apr 2022 10:47:19 -0700 Subject: [PATCH 8/8] fixed minor thing --- requirements.txt | 2 ++ tidy3d/plugins/plotly/app.py | 8 ++++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/requirements.txt b/requirements.txt index 19cccc1b3d..d599d4b382 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,6 +12,8 @@ PyYAML boto3 requests plotly==5.5.0 +dash +jupyter_dash # required to get xarray to not complain dask diff --git a/tidy3d/plugins/plotly/app.py b/tidy3d/plugins/plotly/app.py index 8b38c3942b..513d413b91 100644 --- a/tidy3d/plugins/plotly/app.py +++ b/tidy3d/plugins/plotly/app.py @@ -105,9 +105,9 @@ def _make_layout(self, app: DASH_APP) -> dcc.Tabs: return layout @classmethod - def from_file(cls, path: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ + def from_file(cls, fname: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ """Load the :class:`.SimulationDataApp` from a tidy3d data file in .hdf5 format.""" - sim_data = SimulationData.from_file(path) + sim_data = SimulationData.from_file(fname) sim_data_normalized = sim_data.normalize() return cls(sim_data=sim_data_normalized, mode=mode) @@ -124,7 +124,7 @@ def _make_layout(self, app: DASH_APP) -> dcc.Tabs: return dcc.Tabs([]) @classmethod - def from_file(cls, path: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ + def from_file(cls, fname: str, mode: AppMode = DEFAULT_MODE): # pylint:disable=arguments-differ """Load the SimulationApp from a tidy3d Simulation file in .json or .yaml format.""" - simulation = Simulation.from_file(path) + simulation = Simulation.from_file(fname) return cls(simulation=simulation, mode=mode)