As an immediate solution, I implemented the coordinate transform method the
following paper to simulate bent waveguides in the MPB modesolver:
https://www.osapublishing.org/ol/abstract.cfm?uri=ol-38-11-1778
Essentially, I define a "get_epsilon" method that returns an epsilon tensor
in the transformed coordinates, then pass that as a "material" into each
block in my geometry. The background index also needs to be transformed, so
in the end every pixel in the mesh needs to evaluate this "get_epsilon"
method.
Unfortunately, the run time of the mode solver becomes extremely long when
passing in a get_epsilon method for the material. Each "initializing
epsilon function" step takes around a minute on a 200x100x1 grid, and this
method seems to be called many times when I call "find_k". Overall the
modesolver takes several minutes to execute. This seems strange, since the
geometry isn't changing. Do you have any suggestions for speeding things
up? I am also not getting the mode solutions I expect, but the slow speed
of the simulations seems more pressing, since it prevents me from debugging
further. Code attached.
Thanks again,
Daniel
On Fri, Aug 14, 2020 at 1:54 PM Steven G. Johnson <[email protected]>
wrote:
> See https://github.com/NanoComp/meep/issues/1319
>
> Note you can also implement any source you want on your own by adding a
> volume source with an amplitude function (e.g. interpolated from a file, or
> your own mode solver).
>
> On Aug 14, 2020, at 12:19 PM, Kramnik Daniel <[email protected]>
> wrote:
>
> Thanks, is there any plan to implement cylindrical coordinates?
>
> On Fri, Aug 14, 2020 at 9:14 AM Ardavan Oskooi <[email protected]>
> wrote:
>
>> Meep's eigenmode decomposition (and eigenmode source) feature only
>> supports Cartesian coordinates.
>>
>> On 8/13/20 01:23, Kramnik Daniel wrote:
>> > I'm trying to design some silicon photonic ring resonator structures
>> > using MEEP/MPB and was wondering if the modesolver can work in
>> > cylindrical coordinates. It looks like the FDTD solver supports
>> > cylindrical coordinates, based on the "ring resonator" examples.
>>
> _______________________________________________
> meep-discuss mailing list
> [email protected]
> http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
>
>
>
# Script that designs wrapped couplers for microring resonators using MEEP's modesolver
# The modesolver does not support cylindrical coordinates, so instead we transform the geometry and permittivity tensor from cylindrical to Cartesian coordinates
# See: "Calculation of bending losses for highly confined modes of optical waveguides with transformation optics", Opt. Lett. 2012
import meep as mp
from meep import mpb
import numpy as np
from matplotlib import pyplot as plt
# Use constant-index while developing
#import SiPh_materials.SiPh_materials as materials
# -----------------------------------------------------------
# Setup a simple technology definition
# -----------------------------------------------------------
tech_params = {
'BOX' : None,
'Si' : None,
'Si3N4' : None,
'BEOL' : None,
}
tech_params[ 'BOX' ] = {
'model' : mp.Medium( index = 1.4 ),
'z_min' : -0.200,
'z_mid' : -0.100,
'z_max' : 0.000,
'thickness' : 0.200,
}
tech_params[ 'BEOL' ] = {
'model' : mp.Medium( index = 1.4 ),
'z_min' : 0.000,
'z_mid' : 2.000,
'z_max' : 4.000,
'thickness' : 4.000,
}
tech_params[ 'Si' ] = {
'model' : mp.Medium( index = 3.5 ),
'z_min' : 0.000,
'z_mid' : 0.040,
'z_max' : 0.080,
'thickness' : 0.080,
}
tech_params[ 'Si3N4' ] = {
'model' : mp.Medium( index = 2.0 ),
'z_min' : 0.080,
'z_mid' : 0.090,
'z_max' : 0.100,
'thickness' : 0.020,
}
# -----------------------------------------------------------
# Define ring and bus geometry parameters
# -----------------------------------------------------------
ring_params = {
'outer_radius' : 5.000,
'wg_width' : 1.000,
}
ring_params[ 'inner_radius' ] = ring_params[ 'outer_radius' ] - ring_params[ 'wg_width' ]
bus_params = {
'gap' : 0.050,
'wg_width' : 0.500,
}
bus_params[ 'inner_radius' ] = ring_params[ 'outer_radius' ] + bus_params[ 'gap' ]
bus_params[ 'outer_radius' ] = bus_params[ 'inner_radius' ] + bus_params[ 'wg_width' ]
# -----------------------------------------------------------
# Define utility functions
# -----------------------------------------------------------
# This is the radius that corresponds to u = 0 in the transformed coordinates
# Hardcoded for now
R_0 = 5
def coord_transform_radius_to_u( r_in, R_0 ):
# r_in : radial location to transform
# R_0 : radial location corresponding to u = 0
return R_0 * np.log( r_in / R_0 )
def get_transformed_epsilon( loc, material, frequency ):
u = loc[ 0 ]
v = loc[ 1 ]
w = loc[ 2 ]
epsilon = material.epsilon( freq = frequency )[ 0, 0 ] # Does not support anisotropic materials
epsilon_transformed = epsilon * mp.Vector3( 1, np.exp( 2 * u / R_0 ), 1 )
return mp.Medium( epsilon_diag = epsilon_transformed )
# -----------------------------------------------------------
# Setup simulation geometry and transformations
# -----------------------------------------------------------
sim_wavelength = 1.310
sim_frequency = 1 / sim_wavelength
geometry_lattice = mp.Lattice( size = mp.Vector3( 4, 2, 0 ) )
# Remember, we also need to transform the background! Mesh this first.
block_vacuum = mp.Block(
center = mp.Vector3( 0, 0, 0 ),
size = mp.Vector3( mp.inf, mp.inf, mp.inf ),
material = lambda loc : get_transformed_epsilon( loc, material = mp.Medium( index = 1 ), frequency = sim_frequency )
)
block_BOX = mp.Block(
center = mp.Vector3(
0,
tech_params[ 'BOX' ][ 'z_mid' ],
0
),
size = mp.Vector3(
mp.inf,
tech_params[ 'BOX' ][ 'thickness' ],
mp.inf
),
material = lambda loc : get_transformed_epsilon( loc, material = tech_params[ 'BOX' ][ 'model' ], frequency = sim_frequency )
)
u_inner = coord_transform_radius_to_u( ring_params[ 'inner_radius' ], R_0 )
u_outer = coord_transform_radius_to_u( ring_params[ 'outer_radius' ], R_0 )
block_ring = mp.Block(
center = mp.Vector3(
( u_outer + u_inner ) / 2,
tech_params[ 'Si' ][ 'z_mid' ],
0
),
size = mp.Vector3(
u_outer - u_inner,
tech_params[ 'Si' ][ 'thickness' ],
mp.inf
),
material = lambda loc : get_transformed_epsilon( loc, material = tech_params[ 'Si' ][ 'model' ], frequency = sim_frequency )
)
u_inner = coord_transform_radius_to_u( bus_params[ 'inner_radius' ], R_0 )
u_outer = coord_transform_radius_to_u( bus_params[ 'outer_radius' ], R_0 )
block_bus = mp.Block(
center = mp.Vector3(
( u_outer + u_inner ) / 2,
tech_params[ 'Si' ][ 'z_mid' ],
0
),
size = mp.Vector3(
u_outer - u_inner,
tech_params[ 'Si' ][ 'thickness' ],
mp.inf
),
material = lambda loc : get_transformed_epsilon( loc, material = tech_params[ 'Si' ][ 'model' ], frequency = sim_frequency )
)
get_transformed_epsilon( loc = [ 1, 0, 0 ], material = tech_params[ 'Si' ][ 'model' ], frequency = sim_frequency )
#print( partial( get_transformed_epsilon, material = tech_params[ 'Si' ][ 'model' ], frequency = sim_frequency ).func( ) )
#block_liner =
#block_BEOL =
# Note: order in list defines mesh order. Put oxides first.
#geometry_device = [ block_vacuum, block_BOX, block_BEOL, block_ring, block_bus, block_liner ]
geometry_device = [ block_vacuum, block_ring, block_bus ]
# -----------------------------------------------------------
# Run ModeSolver and plot some results
# -----------------------------------------------------------
num_bands = 2
resolution = 50
num_k = 20
k_min = 0.1
k_max = 2.0
k_points = mp.interpolate( num_k, [ mp.Vector3( k_min ), mp.Vector3( k_max ) ] )
# Instantiate the ModeSolver
ms = mpb.ModeSolver(
geometry_lattice = geometry_lattice,
geometry = geometry_device,
k_points = k_points,
resolution = resolution,
num_bands = num_bands,
)
ms.run( mpb.display_yparities )
# Expecting to find 2 TE modes for the coupled waveguides
k = ms.find_k(
p = mp.NO_PARITY,
omega = sim_frequency,
band_min = 1,
band_max = 2,
korig_and_kdir = mp.Vector3( 0, 0, 1 ),
tol = 1e-4, # tried relaxing this to make things faster
kmag_guess = sim_frequency * 3.4,
kmag_min = sim_frequency * 0.1,
kmag_max = sim_frequency * 4.0
)
epsilon = ms.get_epsilon( )
for band in [ 0, 1 ]:
fig, axis = plt.subplots( 4 )
E = np.squeeze( ms.get_efield( which_band = band, bloch_phase = True ) )
H = np.squeeze( ms.get_hfield( which_band = band, bloch_phase = True ) )
P = np.squeeze( ms.get_tot_pwr( which_band = band ) )
mag_Ex = np.abs( E[ :, :, 0 ].transpose( ) )
mag_Ey = np.abs( E[ :, :, 0 ].transpose( ) )
mag_Ez = np.abs( E[ :, :, 0 ].transpose( ) )
mag_E = np.sqrt( mag_Ex ** 2 + mag_Ey ** 2 + mag_Ez ** 2 )
mag_Hx = np.abs( H[ :, :, 0 ].transpose( ) )
mag_Hy = np.abs( H[ :, :, 0 ].transpose( ) )
mag_Hz = np.abs( H[ :, :, 0 ].transpose( ) )
mag_H = np.sqrt( mag_Hx ** 2 + mag_Hy ** 2 + mag_Hz ** 2 )
axis[ 0 ].imshow( epsilon.transpose( ), interpolation = None, cmap = 'gist_heat', alpha = 1 )
axis[ 1 ].imshow( mag_E, cmap = 'gist_heat', alpha = 1 )
axis[ 2 ].imshow( mag_H, cmap = 'gist_heat', alpha = 1 )
axis[ 3 ].imshow( np.abs( np.squeeze( P ).transpose( ) ), cmap = 'gist_heat', alpha = 1 )
plt.show( )
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
