Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

tidy3d pluging to allow for symmetries #268

Open
pascualmunoz opened this issue Dec 5, 2023 · 12 comments · Fixed by #269
Open

tidy3d pluging to allow for symmetries #268

pascualmunoz opened this issue Dec 5, 2023 · 12 comments · Fixed by #269
Assignees
Labels
enhancement New feature or request

Comments

@pascualmunoz
Copy link

Is your feature request related to a problem? Please describe.

from collections.abc import Callable
from functools import partial

from pydantic import BaseModel

import gdsfactory as gf
from gdsfactory.add_pins import add_pin_rectangle_inside
from gdsfactory.component import Component
from gdsfactory.config import CONF
from gdsfactory.cross_section import cross_section
from gdsfactory.technology import (
    LayerLevel,
    LayerStack,
    LayerView,
    LayerViews,
    LayerMap,
)
from gdsfactory.typings import Layer
from gdsfactory.config import print_version_pdks, print_version_plugins
from gdsfactory.generic_tech import get_generic_pdk

gf.config.rich_output()
nm = 1e-3

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td

import gplugins.tidy3d as gt
from gplugins import plot
from gplugins.common.config import PATH
import gplugins as gp


mat_Si = td.material_library["cSi"]["Li1993_293K"] # my_fitter(wvl_v, mat.si(wvl_v))
mat_SiO2 = td.material_library["SiO2"]["Horiba"] # my_fitter(wvl_v, mat.sio2(wvl_v))


class LayerMapMY(LayerMap):
    WG: Layer = (1, 0)

LAYER = LayerMapMY()

class MYLayerViews(LayerViews):
    WG: LayerView = LayerView(color="gold")

LAYER_VIEWS = MYLayerViews(layer_map=dict(LAYER))

def get_layer_stack_my(
        thickness_wg: float = 220 * nm,
) -> LayerStack:
    return LayerStack(
        layers=dict(
            strip=LayerLevel(
                layer=LAYER.WG,
                thickness=thickness_wg,
                zmin=0.0,
                material="Si",
            ),
        )
    )

LAYER_STACK = get_layer_stack_my()

WIDTH = 0.5

# Specify a cross_section to use
cs_strip = partial(gf.cross_section.cross_section, width=WIDTH, layer=LAYER.WG)

mmi1x2 = partial(
    gf.components.mmi1x2,
    width=WIDTH,
    width_taper=WIDTH,
    width_mmi=3 * WIDTH,
    cross_section=cs_strip,
)

generic_pdk = get_generic_pdk()

my_pdk = gf.Pdk(
     name="MY_PDK",
     cells=dict(mmi1x2=mmi1x2),
     cross_sections=dict(strip=cs_strip),
     layers=dict(LAYER),
     base_pdk=generic_pdk,
     layer_views=LAYER_VIEWS,
     layer_stack=LAYER_STACK,
)

my_pdk.activate()

sw_test = gf.components.straight(
    length = 10.0,
    width= 0.5,
    cross_section="strip",
)

sw_test.plot()

md = {
    "Si": mat_Si,
    "sio2": mat_SiO2,
}

sp = gt.write_sparameters(sw_test, plot_simulation_layer_name="strip", layer_stack=LAYER_STACK, material_mapping=md)
# lets plot the fundamental input mode
sp = gt.write_sparameters(
    sw_test,  
    plot_mode_port_name="o1", 
    plot_mode_index=0, 
    layer_stack=LAYER_STACK, 
    material_mapping=md,
)
sp = gt.write_sparameters(
    component = sw_test,
    layer_stack = LAYER_STACK,
    material_mapping = md, 
    symmetry = (0,-1,1),
)

Describe the solution you'd like
Result is:

--> 114 sp = gt.write_sparameters(
115 component = sw_test,
116 layer_stack = LAYER_STACK,
117 material_mapping = md,
118 symmetry = (0,-1,1),
119 )

TypeError: write_sparameters() got an unexpected keyword argument 'symmetry'

I have looked into the different .py files for the tidy3d plugin, and to my very poor understanding it seems the native Tidy3D symmetry feature is not included in gplugin/tidy3d.

If that is the case, please consider including it (saves time-money).

@pascualmunoz pascualmunoz added the enhancement New feature or request label Dec 5, 2023
@joamatab
Copy link
Contributor

joamatab commented Dec 5, 2023

yes, it would be great to add it

Also it would be to get support directly from tidy3d team to review support and improve this plugin :)

@tylerflex
@heitzmann
@yaugenst
@momchil-flex

@joamatab
Copy link
Contributor

joamatab commented Dec 5, 2023

Pascual,

can you pip install gplugins --upgrade and try again?

@pascualmunoz
Copy link
Author

Dear all,

I just run the very same script above and now the output to the very last command (having symmetries) is:

14:42:36 CET ERROR: Mode object 'name='o1' type='ModeSource'                    
             center=(-2.3290386521308224, 2.816687638038912e-16, 0.11)          
             size=(0.0, 2.0, 1.5) source_time=GaussianPulse(amplitude=1.0,      
             phase=0.0, type='GaussianPulse', freq0=193710487003340.7,          
             fwidth=25061020522466.062, offset=5.0, remove_dc_component=True)   
             num_freqs=1 direction='+' mode_spec=ModeSpec(num_modes=1,          
             target_neff=None, num_pml=(0, 0), filter_pol='te', angle_theta=0.0,
             angle_phi=0.0, precision='single', bend_radius=None,               
             bend_axis=None, track_freq='central', group_index_step=False,      
             type='ModeSpec') mode_index=0' (at `simulation.sources[0]`) in     
             presence of symmetries must be in the main quadrant, or centered on
             the symmetry axis.                                                 
             ERROR: Mode object 'type='ModeMonitor' center=(-2.3,               
             2.816687638038912e-16, 0.11) size=(0.0, 2.0, 1.5) name='o1'        
             interval_space=(1, 1, 1) colocate=False freqs=(206753419310344.84, 
             205337300000000.0, 203940447619047.62, 202562471621621.62,         
             201202991946308.72, 199861638666666.66, 198538051655629.12,        
             197231880263157.9, 195942783006535.94, 194670427272727.28,         
             193414489032258.06, 192174652564102.56, 190950610191082.78,        
             189742062025316.44, 188548715723270.44, 187370286250000.0,         
             186206495652173.9, 185057072839506.16, 183921753374233.12,         
             182800279268292.66, 181692398787878.78)                            
             apodization=ApodizationSpec(start=None, end=None, width=None,      
             type='ApodizationSpec') mode_spec=ModeSpec(num_modes=1,            
             target_neff=None, num_pml=(0, 0), filter_pol='te', angle_theta=0.0,
             angle_phi=0.0, precision='single', bend_radius=None,               
             bend_axis=None, track_freq='central', group_index_step=False,      
             type='ModeSpec')' (at `simulation.monitors[0]`) in presence of     
             symmetries must be in the main quadrant, or centered on the        
             symmetry axis.                                                     
---------------------------------------------------------------------------
ValidationError                           Traceback (most recent call last)
[/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/04_tidy3D_symmetry_straight_waveguide.ipynb](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/04_tidy3D_symmetry_straight_waveguide.ipynb) Celda 3 line 1
----> 1 sp = gt.write_sparameters(
      2     component = sw_test,
      3     layer_stack = LAYER_STACK,
      4     material_mapping = md,  
      5     symmetry = (0,-1,1),
      6 )

File [~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:597](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:597), in write_sparameters(component, layer_stack, material_mapping, extend_ports, port_offset, pad_xy_inner, pad_xy_outer, pad_z_inner, pad_z_outer, dilation, wavelength, bandwidth, num_freqs, min_steps_per_wvl, center_z, sim_size_z, port_size_mult, run_only, element_mappings, extra_monitors, mode_spec, boundary_spec, symmetry, run_time, shutoff, folder_name, dirpath, verbose, plot_simulation_layer_name, plot_simulation_port_index, plot_simulation_z, plot_simulation_x, plot_mode_index, plot_mode_port_name, plot_epsilon, filepath, overwrite, **kwargs)
    [595](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:595) else:
    [596](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:596)     time.sleep(0.2)
--> [597](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:597)     s = modeler.run()
    [598](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:598)     for port_in in s.port_in.values:
    [599](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:599)         for port_out in s.port_out.values:

File [~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:501](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:501), in ComponentModeler.run(self, path_dir)
    [498](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:498) """Solves for the scattering matrix of the system."""
    [499](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:499) path_dir = self.get_path_dir(path_dir)
--> [501](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:501) batch_data = self._run_sims(path_dir=path_dir)
    [502](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:502) return self._construct_smatrix(batch_data=batch_data)

File [~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:373](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:373), in ComponentModeler._run_sims(self, path_dir)
    [371](https://file+.vscode-resource.vscode-cdn.net/Users/pascual/ownCloud/PASCUAL_GDSFACTORY/20231124-tidy3D/~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:371) def _run_sims(self, path_dir: str = DEFAULT_DATA_DIR) -> BatchData:
...
  Mode object 'type='ModeMonitor' center=(-2.3, 2.816687638038912e-16, 0.11) size=(0.0, 2.0, 1.5) name='o1' interval_space=(1, 1, 1) colocate=False freqs=(206753419310344.84, 205337300000000.0, 203940447619047.62, 202562471621621.62, 201202991946308.72, 199861638666666.66, 198538051655629.12, 197231880263157.9, 195942783006535.94, 194670427272727.28, 193414489032258.06, 192174652564102.56, 190950610191082.78, 189742062025316.44, 188548715723270.44, 187370286250000.0, 186206495652173.9, 185057072839506.16, 183921753374233.12, 182800279268292.66, 181692398787878.78) apodization=ApodizationSpec(start=None, end=None, width=None, type='ApodizationSpec') mode_spec=ModeSpec(num_modes=1, target_neff=None, num_pml=(0, 0), filter_pol='te', angle_theta=0.0, angle_phi=0.0, precision='single', bend_radius=None, bend_axis=None, track_freq='central', group_index_step=False, type='ModeSpec')' (at `simulation.monitors[0]`) in presence of symmetries must be in the main quadrant, or centered on the symmetry axis. (type=value_error.setup)
grid_spec
  'NoneType' object is not iterable (type=type_error)
normalize_index
  object of type 'NoneType' has no len() (type=type_error)

@tylerflex
Copy link

Your mode source and monitors name="o1" must be centered on the symmetry planes, which lie on the dim=0 plane where symmetry != 0. Could you check the simulation geometry without symmetry and make sure that the sources and monitors are centered at 0 along the dimensions where you have symmetry defined? Perhaps you don't want symmetry in z if there's a substrate?

@pascualmunoz
Copy link
Author

Hi, I run the simulation now with:

sp = gt.write_sparameters(
    component = sw_test,
    layer_stack = LAYER_STACK,
    material_mapping = md,  
    symmetry = (0,-1,0),
)

and the result error is exactly the same for the source and monitor. Having a closer look to the y coordinate for both (see above) is 2.816687638038912e-16 ... which to me is zero but, is it for tidy3d?

Best.

@tylerflex
Copy link

Hm, maybe @momchil-flex can comment on the numerical precision handling for symmetry.

@yaugenst
Copy link
Collaborator

Hi all, so as far as I can see there are two separate issues here:

  1. Rounding error: I think it makes sense to round layer centers as well as port locations by default. This is addressed in Round port locations and layer centers to one picometer by default #273 with a default of one picometer.
  2. The other issue is that the SWG layer in @pascualmunoz's example is not centered at zero (in z), see picture:
    screenshot_20231211_083148
    This is because the layer center is at 0.11, and the simulation is centered at 0 (the default). To have the component centered at zero, the component layer's zmin should be at 0 - thickness_wg / 2, i.e. -0.11. The alternative is to center the simulation at the component layer center, which can be done either explicity passing center_z=0.11 or by centering on the layer name, i.e. passing center_z="strip" to write_sparameters(). We can try to implement some logic to do this by default, but I don't see it being very robust as it is hard to guess (in the general case) where the simulation should be centered, as more than one layer can contain geometries. I guess we could implement it for the special case where only a single layer in the layerstack contains non-trivial geometries.

@yaugenst
Copy link
Collaborator

Another point is perhaps that in general, for essentially all components that are not straight waveguides or tapers, symmetries (other than z) are tricky to handle. For example, a 1x2 splitter could exploit a mirror plane only for the simulation of o1, and in that case the ComponentModeler would need to know about the symmetry and "virtually mirror" the output port that is not technically within the simulation. For the other two ports, symmetry would not apply...

@yaugenst yaugenst self-assigned this Dec 11, 2023
@pascualmunoz
Copy link
Author

Hi all, so as far as I can see there are two separate issues here:

  1. Rounding error: I think it makes sense to round layer centers as well as port locations by default. This is addressed in Round port locations and layer centers to one picometer by default #273 with a default of one picometer.
  2. The other issue is that the SWG layer in @pascualmunoz's example is not centered at zero (in z), see picture:
    screenshot_20231211_083148
    This is because the layer center is at 0.11, and the simulation is centered at 0 (the default). To have the component centered at zero, the component layer's zmin should be at 0 - thickness_wg / 2, i.e. -0.11. The alternative is to center the simulation at the component layer center, which can be done either explicity passing center_z=0.11 or by centering on the layer name, i.e. passing center_z="strip" to write_sparameters(). We can try to implement some logic to do this by default, but I don't see it being very robust as it is hard to guess (in the general case) where the simulation should be centered, as more than one layer can contain geometries. I guess we could implement it for the special case where only a single layer in the layerstack contains non-trivial geometries.

Hi, thanks @yaugenst for the comments. I agree it makes little sense for automatic centering at z = 0, and apologies for this not being there but at z = 0.11 µm, it comes from some other example. For the case then z-symmetry doesn't make any sense, so my apologies once more.

@pascualmunoz
Copy link
Author

All, in my quest to save time-money on simulations, I have read the tidy3d documentation on S-parameters:

https://www.flexcompute.com/tidy3d/examples/notebooks/SMatrix/

With slight modifications in the example above:

inp = ("o1", 0)
out = ("o2", 0)

c_ii = (inp, inp)
c_io = (inp, out)
c_oi = (out, inp)
c_oo = (out, out)

map_xx = (c_ii, c_oo, +1)
map_xy = (c_io, c_oi, +1)

element_mappings = (map_xx, map_xy)
sp_port_sym_em_ro = gt.write_sparameters(
    component = sw_test,
    layer_stack = LAYER_STACK,
    material_mapping = md,  
    symmetry = (0,-1,0),
    element_mappings=element_mappings,
    run_only = ( inp, )
)

I was able to reduce the use of FlexCredits (computation time is equivalent to the same without "element_mappings" and without "run_only" though).

However the call to write_sparameters does not end, and the output is the following:

17:30:45 CET loading simulation from                                            
             /Users/pascual/.gdsfactory/sparameters/8b77f32deba9d3328164504fcffc
             774fd72de6843311afd6799ea80caea23a14/fdve-fd324b19-8d0b-451d-9430-5
             b7b93ab5d12.hdf5

{
	"name": "KeyError",
	"message": "\"not all values found in index 'port_in'. Try setting the `method` keyword argument (example: method='nearest').\"",
	"stack": "---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3790, in Index.get_loc(self, key)
   3789 try:
-> 3790     return self._engine.get_loc(casted_key)
   3791 except KeyError as err:

File index.pyx:152, in pandas._libs.index.IndexEngine.get_loc()

File index.pyx:181, in pandas._libs.index.IndexEngine.get_loc()

File pandas/_libs/hashtable_class_helper.pxi:7080, in pandas._libs.hashtable.PyObjectHashTable.get_item()

File pandas/_libs/hashtable_class_helper.pxi:7088, in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'o2'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
File ~/anaconda3/lib/python3.11/site-packages/xarray/core/indexes.py:486, in PandasIndex.sel(self, labels, method, tolerance)
    485 try:
--> 486     indexer = self.index.get_loc(label_value)
    487 except KeyError as e:

File ~/anaconda3/lib/python3.11/site-packages/pandas/core/indexes/base.py:3797, in Index.get_loc(self, key)
   3796         raise InvalidIndexError(key)
-> 3797     raise KeyError(key) from err
   3798 except TypeError:
   3799     # If we have a listlike key, _check_indexing_error will raise
   3800     #  InvalidIndexError. Otherwise we fall through and re-raise
   3801     #  the TypeError.

KeyError: 'o2'

The above exception was the direct cause of the following exception:

[... omitted ...]

File ~/anaconda3/lib/python3.11/site-packages/gplugins/tidy3d/component.py:606, in write_sparameters(component, layer_stack, material_mapping, extend_ports, port_offset, pad_xy_inner, pad_xy_outer, pad_z_inner, pad_z_outer, dilation, wavelength, bandwidth, num_freqs, min_steps_per_wvl, center_z, sim_size_z, port_size_mult, run_only, element_mappings, extra_monitors, mode_spec, boundary_spec, symmetry, run_time, shutoff, folder_name, dirpath, verbose, plot_simulation_layer_name, plot_simulation_port_index, plot_simulation_z, plot_simulation_x, plot_mode_index, plot_mode_port_name, plot_epsilon, filepath, overwrite, **kwargs)
    604 else:
    605     time.sleep(0.2)
--> 606     s = modeler.run()
    607     for port_in in s.port_in.values:
    608         for port_out in s.port_out.values:

File ~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:502, in ComponentModeler.run(self, path_dir)
    499 path_dir = self.get_path_dir(path_dir)
    501 batch_data = self._run_sims(path_dir=path_dir)
--> 502 return self._construct_smatrix(batch_data=batch_data)

File ~/anaconda3/lib/python3.11/site-packages/tidy3d/plugins/smatrix/smatrix.py:493, in ComponentModeler._construct_smatrix(self, batch_data)
    486     port_in_to, mode_index_in_to = col_out
    487     coords_to = dict(
    488         port_in=port_in_to,
    489         mode_index_in=mode_index_in_to,
    490         port_out=port_out_to,
    491         mode_index_out=mode_index_out_to,
    492     )
--> 493     s_matrix.loc[coords_to] = mult_by * s_matrix.loc[coords_from].values
    495 return s_matrix

File ~/anaconda3/lib/python3.11/site-packages/xarray/core/dataarray.py:222, in _LocIndexer.__setitem__(self, key, value)
    219     labels = indexing.expanded_indexer(key, self.data_array.ndim)
    220     key = dict(zip(self.data_array.dims, labels))
--> 222 dim_indexers = map_index_queries(self.data_array, key).dim_indexers
    223 self.data_array[dim_indexers] = value

File ~/anaconda3/lib/python3.11/site-packages/xarray/core/indexing.py:190, in map_index_queries(obj, indexers, method, tolerance, **indexers_kwargs)
    188         results.append(IndexSelResult(labels))
    189     else:
--> 190         results.append(index.sel(labels, **options))
    192 merged = merge_sel_results(results)
    194 # drop dimension coordinates found in dimension indexers
    195 # (also drop multi-index if any)
    196 # (.sel() already ensures alignment)

File ~/anaconda3/lib/python3.11/site-packages/xarray/core/indexes.py:488, in PandasIndex.sel(self, labels, method, tolerance)
    486                 indexer = self.index.get_loc(label_value)
    487             except KeyError as e:
--> 488                 raise KeyError(
    489                     f\"not all values found in index {coord_name!r}. \"
    490                     \"Try setting the `method` keyword argument (example: method='nearest').\"
    491                 ) from e
    493 elif label_array.dtype.kind == \"b\":
    494     indexer = label_array

KeyError: \"not all values found in index 'port_in'. Try setting the `method` keyword argument (example: method='nearest').\""
}         

What do you suggest?

Thanks once more.

@yaugenst
Copy link
Collaborator

It looks like it's not finding the data for port o2, which makes sense since you specified run_only=(("o1", 0),). I don't think you need to use both run_only and an element_mapping in this case, the mapping should already take care of the ports that are needed to run.

@pascualmunoz
Copy link
Author

It looks like it's not finding the data for port o2, which makes sense since you specified run_only=(("o1", 0),). I don't think you need to use both run_only and an element_mapping in this case, the mapping should already take care of the ports that are needed to run.

Hi @yaugenst thanks for the comment. I assumed the same but I had run several simulations, first one without element mapping and second one with the mapping. They took the same (FlexCredits). So I added the run_only for one port and then the FlexCredits halved, the simulation completed, but the retrieving / assembling of results fails as detailed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants