diff --git a/gplugins/gmeep/get_meep_geometry.py b/gplugins/gmeep/get_meep_geometry.py index fba89e7e..af9c5324 100644 --- a/gplugins/gmeep/get_meep_geometry.py +++ b/gplugins/gmeep/get_meep_geometry.py @@ -16,7 +16,6 @@ def get_meep_geometry_from_component( layer_stack: LayerStack | None = None, material_name_to_meep: dict[str, str | float] | None = None, wavelength: float = 1.55, - is_3d: bool = False, dispersive: bool = False, exclude_layers: LayerSpecs | None = None, **kwargs, @@ -28,7 +27,6 @@ def get_meep_geometry_from_component( layer_stack: for material layers. material_name_to_meep: maps layer_stack name to meep material name. wavelength: in um. - is_3d: renders in 3D. dispersive: add dispersion. exclude_layers: these layers are ignored during geometry creation. kwargs: settings. @@ -62,10 +60,9 @@ def get_meep_geometry_from_component( if layer_index in exclude_layers or layer_index not in polygons_per_layer: continue - zmin_um = level.zmin if is_3d else 0 - sw_angle = np.pi * level.sidewall_angle / 180 if is_3d else 0 + sw_angle = np.pi * level.sidewall_angle / 180 for polygon in polygons_per_layer[layer_index]: - vertices = [mp.Vector3(p[0], p[1], zmin_um) for p in polygon] + vertices = [mp.Vector3(p[0], p[1], level.zmin) for p in polygon] material_name = level.material if material_name: @@ -79,7 +76,7 @@ def get_meep_geometry_from_component( mp.Prism( vertices=vertices, height=level.thickness, - sidewall_angle=sw_angle, + sidewall_angle=sw_angle, # TODO: libctl has issues with slanted prisms -> incorrect geometry material=material, ) ) diff --git a/gplugins/gmeep/get_simulation.py b/gplugins/gmeep/get_simulation.py index 38f2f247..63f9b3fe 100644 --- a/gplugins/gmeep/get_simulation.py +++ b/gplugins/gmeep/get_simulation.py @@ -4,7 +4,7 @@ import inspect import warnings -from typing import Any +from typing import Any, Literal import gdsfactory as gf import meep as mp @@ -12,7 +12,7 @@ from gdsfactory.component import Component from gdsfactory.pdk import get_layer_stack from gdsfactory.technology import LayerStack -from gdsfactory.typings import LayerSpecs +from gdsfactory.typings import LayerSpecs, Float3 from gplugins.common.base_models.component import move_polar_rad_copy from gplugins.gmeep.get_material import get_material @@ -26,6 +26,19 @@ settings_meep = set(sig.parameters.keys()) +def is_point_in_plane( + test_point: Float3, + plane_support: Float3, + plane_normal: Float3, + tolerance: float = 1e-6, # 1 pm for coordinates in µm +): + a, b, c = plane_normal + xt, yt, zt = test_point + x0, y0, z0 = plane_support + distance = (a*(xt-x0) + b*(yt-y0) + c*(zt-z0)) / np.linalg.norm(plane_normal) + return bool(abs(distance) <= tolerance) + + def get_simulation( component: Component, resolution: int = 30, @@ -36,6 +49,8 @@ def get_simulation( tpml: float = 1.5, clad_material: str = "SiO2", is_3d: bool = False, + normal_2d: Literal['X', 'Y', 'Z'] = 'Z', + point_2d: tuple[float, float, float] = (0, 0, 0), wavelength_start: float = 1.5, wavelength_stop: float = 1.6, wavelength_points: int = 50, @@ -105,6 +120,8 @@ def get_simulation( tpml: PML thickness (um). clad_material: material for cladding. is_3d: if True runs in 3D. + normal_2d: specified normal of 2D simulation plane + point_2d: specifies support point for 2D simulation plane wavelength_start: wavelength min (um). wavelength_stop: wavelength max (um). wavelength_points: wavelength steps. @@ -141,23 +158,19 @@ def get_simulation( for setting in settings: if setting not in settings_meep: raise ValueError(f"{setting!r} not in {sorted(settings_meep)}") + normal_2d = normal_2d.upper() + normal_vec = {'X': (1, 0, 0), 'Y': (0, 1, 0), 'Z': (0, 0, 1)}[normal_2d] layer_stack = layer_stack or get_layer_stack() layer_to_thickness = layer_stack.get_layer_to_thickness() - dummy_component = gf.Component() - component_ref = dummy_component << component - component_ref.x = 0 - component_ref.y = 0 - wavelength = (wavelength_start + wavelength_stop) / 2 - wavelengths = np.linspace(wavelength_start, wavelength_stop, wavelength_points) - port_names = [port.name for port in component_ref.ports] + port_names = [port.name for port in component.ports] if port_source_name not in port_names: warnings.warn(f"port_source_name={port_source_name!r} not in {port_names}") - port_source = component_ref.ports[0] + port_source = component.ports[0] port_source_name = port_source.name warnings.warn(f"Selecting port_source_name={port_source_name!r} instead.") @@ -165,14 +178,7 @@ def get_simulation( f"component needs to be a gf.Component, got Type {type(component)}" ) - component_extended = ( - gf.c.extend_ports( - component=component, length=extend_ports_length, centered=True - ) - if extend_ports_length - else component - ) - + component_extended = gf.c.extend_ports(component=component, length=extend_ports_length) component_extended = component_extended.copy() component_extended.flatten() @@ -191,12 +197,12 @@ def get_simulation( t_core = sum( layers_thickness ) # This isn't exactly what we want but I think it's better than max - cell_thickness = tpml + zmargin_bot + t_core + zmargin_top + tpml if is_3d else 0 + cell_thickness = tpml + zmargin_bot + t_core + zmargin_top + tpml cell_size = mp.Vector3( - component.xsize + 2 * tpml, - component.ysize + 2 * tpml, - cell_thickness, + 0 if normal_2d == 'X' and not is_3d else component.xsize + 2 * tpml, + 0 if normal_2d == 'Y' and not is_3d else component.ysize + 2 * tpml, + 0 if normal_2d == 'Z' and not is_3d else cell_thickness, ) geometry = get_meep_geometry_from_component( @@ -204,7 +210,6 @@ def get_simulation( layer_stack=layer_stack, material_name_to_meep=material_name_to_meep, wavelength=wavelength, - is_3d=is_3d, dispersive=dispersive, exclude_layers=exclude_layers, ) @@ -214,19 +219,23 @@ def get_simulation( frequency_width = dfcen * fcen # Add source - port = component_ref.ports[port_source_name] + port = component.ports[port_source_name] angle_rad = np.radians(port.orientation) width = port.width + 2 * port_margin size_x = width * abs(np.sin(angle_rad)) size_y = width * abs(np.cos(angle_rad)) size_x = 0 if size_x < 0.001 else size_x size_y = 0 if size_y < 0.001 else size_y - size_z = cell_thickness - 2 * tpml if is_3d else 20 - size = [size_x, size_y, size_z] + size_z = cell_thickness - 2 * tpml + size = [ + 0 if normal_2d == 'X' and not is_3d else size_x, + 0 if normal_2d == 'Y' and not is_3d else size_y, + 0 if normal_2d == 'Z' and not is_3d else size_z, + ] xy_shifted = move_polar_rad_copy( np.array(port.center), angle=angle_rad, length=port_source_offset ) - center = xy_shifted.tolist() + [0] # (x, y, z=0) + center = xy_shifted.round(6).tolist() + [0] # (x, y, z=0) if np.isclose(port.orientation, 0): direction = mp.X @@ -242,6 +251,11 @@ def get_simulation( "not 0, 90, 180, 270 degrees" ) + if not (is_3d or is_point_in_plane(center, point_2d, normal_vec)): + raise ValueError( + f"Source '{port_source_name}' (center={center}) is not in {normal_2d}-normal 2D simulation domain around {point_2d}." + ) + sources = [ mp.EigenModeSource( src=mp.ContinuousSource(fcen) @@ -250,15 +264,21 @@ def get_simulation( size=size, center=center, eig_band=port_source_mode + 1, - eig_parity=mp.NO_PARITY if is_3d else mp.EVEN_Y + mp.ODD_Z, + eig_parity=mp.NO_PARITY, eig_match_freq=True, eig_kpoint=-1 * mp.Vector3(x=1).rotate(mp.Vector3(z=1), angle_rad), direction=direction, ) ] + sim_center = mp.Vector3( + point_2d[0] if normal_2d == 'X' and not is_3d else component.x, + point_2d[1] if normal_2d == 'Y' and not is_3d else component.y, + point_2d[2] if normal_2d == 'Z' and not is_3d else 0, + ) sim = mp.Simulation( cell_size=cell_size, + geometry_center=sim_center, boundary_layers=[mp.PML(tpml)], sources=sources, geometry=geometry, @@ -273,7 +293,7 @@ def get_simulation( # Add port monitors dict monitors = {} - for port in component_ref.ports: + for port in component.ports: port_name = port.name angle_rad = np.radians(port.orientation) width = port.width + 2 * port_margin @@ -281,8 +301,11 @@ def get_simulation( size_y = width * abs(np.cos(angle_rad)) size_x = 0 if size_x < 0.001 else size_x size_y = 0 if size_y < 0.001 else size_y - size = mp.Vector3(size_x, size_y, size_z) - size = [size_x, size_y, size_z] + size = mp.Vector3( + 0 if normal_2d == 'X' and not is_3d else size_x, + 0 if normal_2d == 'Y' and not is_3d else size_y, + 0 if normal_2d == 'Z' and not is_3d else size_z, + ) # if monitor has a source move monitor inwards length = ( @@ -293,10 +316,14 @@ def get_simulation( xy_shifted = move_polar_rad_copy( np.array(port.center), angle=angle_rad, length=length ) - center = xy_shifted.tolist() + [0] # (x, y, z=0) - m = sim.add_mode_monitor(freqs, mp.ModeRegion(center=center, size=size)) - m.z = 0 - monitors[port_name] = m + center = xy_shifted.round(6).tolist() + [0] # (x, y, z=0) + if is_3d or is_point_in_plane(center, point_2d, normal_vec): + m = sim.add_mode_monitor(freqs, mp.ModeRegion(center=center, size=size)) + m.z = 0 + monitors[port_name] = m + else: + warnings.warn(f"Monitor at port '{port_name}' ignored, " + f"because it is not in the {normal_2d}-normal 2D simulation domain around {point_2d}.") return dict( sim=sim, cell_size=cell_size, diff --git a/gplugins/gmeep/test_materials.py b/gplugins/gmeep/test_materials.py index 543ace76..afa65c2f 100644 --- a/gplugins/gmeep/test_materials.py +++ b/gplugins/gmeep/test_materials.py @@ -3,6 +3,8 @@ from gplugins.common.utils import optical_constants from gplugins.gmeep.write_sparameters_meep import write_sparameters_meep +gf.gpdk.PDK.activate() + def test_materials_override() -> None: """Checks that materials are properly overridden if index is provided.""" diff --git a/gplugins/gmeep/test_write_sparameters_meep.py b/gplugins/gmeep/test_write_sparameters_meep.py index ec9a36ed..c43400aa 100644 --- a/gplugins/gmeep/test_write_sparameters_meep.py +++ b/gplugins/gmeep/test_write_sparameters_meep.py @@ -8,6 +8,7 @@ import gplugins as sim import gplugins.gmeep as gm +gf.gpdk.PDK.activate() simulation_settings = dict(resolution=20, is_3d=False) diff --git a/gplugins/gmeep/write_sparameters_meep.py b/gplugins/gmeep/write_sparameters_meep.py index 366b189a..746de668 100644 --- a/gplugins/gmeep/write_sparameters_meep.py +++ b/gplugins/gmeep/write_sparameters_meep.py @@ -19,7 +19,7 @@ from gdsfactory.pdk import get_layer_stack from gdsfactory.serialization import clean_value_json from gdsfactory.technology import LayerStack -from gdsfactory.typings import ComponentSpec, PathType, Port, PortSymmetries, LayerSpec +from gdsfactory.typings import ComponentSpec, PathType, Ports, PortSymmetries, LayerSpec from tqdm.auto import tqdm from gplugins.common.utils import port_symmetries @@ -52,14 +52,15 @@ def remove_simulation_kwargs(d: dict[str, Any]) -> dict[str, Any]: def parse_port_eigenmode_coeff( - port_name: str, ports: dict[str, Port], sim_dict: dict, port_mode: int = 0 + port_name: str, ports: Ports, sim_dict: dict, port_mode: int = 0 ): """Returns the coefficients relative to whether the wavevector is entering or exiting simulation. Args: - port_index: index of port. + port_name: name of port. ports: component_ref.ports. sim_dict: simulation dict. + port_mode: mode to get coefficient for """ if port_name not in ports: @@ -384,7 +385,7 @@ def write_sparameters_meep( sim.plot2D( output_plane=mp.Volume( size=mp.Vector3(sim.cell_size.x, sim.cell_size.y, 0), - center=mp.Vector3(0, 0, z), + center=mp.Vector3(component.x, component.y, z), ), **plot_args, ) @@ -406,7 +407,6 @@ def sparameter_calculation( port_source_name: str, component: Component, port_symmetries: PortSymmetries | None = port_symmetries, - port_names: list[str] = port_names, port_source_mode: int = 0, wavelength_start: float = wavelength_start, wavelength_stop: float = wavelength_stop, @@ -438,11 +438,8 @@ def sparameter_calculation( ) sim = sim_dict["sim"] - # freqs = sim_dict["freqs"] - # wavelengths = 1 / freqs - # print(sim.resolution) - # Terminate when the area in the whole area decayed + # Terminate when the energy in the whole area decayed termination = [mp.stop_when_energy_decayed(dt=50, decay_by=decay_by)] if animate: @@ -493,7 +490,7 @@ def sparameter_calculation( port_source_name, component.ports, sim_dict, port_mode=port_source_mode ) # Get coefficients - for port_name in port_names: + for port_name in sim_dict['monitors'].keys(): for port_mode in port_modes: _, monitor_exiting = parse_port_eigenmode_coeff( port_name, @@ -535,7 +532,6 @@ def sparameter_calculation( wavelength_stop=wavelength_stop, wavelength_points=wavelength_points, animate=animate, - port_names=port_names, **settings, ) # Synchronize dicts @@ -568,7 +564,6 @@ def sparameter_calculation( wavelength_stop=wavelength_stop, wavelength_points=wavelength_points, animate=animate, - port_names=port_names, **settings, ) )