kch3782 / torcwa

GPU-accelerated RCWA with automatic differentiation
Other
104 stars 22 forks source link

Fourier order do't converge for metals #17

Closed MahmutRuzi closed 1 year ago

MahmutRuzi commented 1 year ago

Hello Thanks for the code and quick response to the comments & issues raised. Have you tried any convergence study for metallic structures ? I tried it. Please see the notebook: https://github.com/MahmutRuzi/photonics-sim-examples/blob/main/torcwa_AuNanobar2.ipynb

I tried the Fourier truncation from Nx=Ny=9 up to Nx=Ny=22 for the mentioned example. The calculated reflectance is very erratic with respect to NxNy. I could try increasing the terms further, but not sure even that lead to convergence within reasonable time. The Matlab code RETICOLO converges for NxNy=15 for the sample problem. It seems they implemented the enhanced transmittance method for fast convergence. Can you please comment if the code as it stands now suitable for such simulations? If not, wow it can be modified?
An approach seems to be using adaptive spatial resolution method (https://doi.org/10.1364/JOSAA.28.000238) to sample fine around the dielectric/metal interface.

Another question not so related question: is it possible to import planar structures (GDSII for example) for the code instead of the predefined rectangles, ellipses, etc?

Thank you

kch3782 commented 1 year ago

Currently, TORCWA uses the basic RCWA algorithm, and this algorithm shows slow convergence for metallic structures. Below is the simulation result when the following code is applied. (Transmission and reflection of square gold pillar with respect to Nx, Ny)

lamb0 = torch.tensor(532.,dtype=geo_dtype,device=device) 
inc_ang = 0.*(np.pi/180)
azi_ang = 0.*(np.pi/180)

substrate_eps = 1.46**2
Au_eps = (0.54386+2.2309j)**2

L = [300., 300.]
torcwa.rcwa_geo.dtype = geo_dtype
torcwa.rcwa_geo.device = device
torcwa.rcwa_geo.Lx = L[0]
torcwa.rcwa_geo.Ly = L[1]
torcwa.rcwa_geo.nx = 300
torcwa.rcwa_geo.ny = 300
torcwa.rcwa_geo.grid()
torcwa.rcwa_geo.edge_sharpness = 1000.

layer0_geometry = torcwa.rcwa_geo.square(W=120.,Cx=L[0]/2.,Cy=L[1]/2.)
layer0_eps = layer0_geometry*Au_eps + (1.-layer0_geometry)
layer0_thickness = 300.

txx = []
rxx = []
# Generate and perform simulation
for order_N in range(20):
    order = [order_N,order_N]
    sim = torcwa.rcwa(freq=1/lamb0,order=order,L=L,dtype=sim_dtype,device=device)
    sim.add_output_layer(eps=substrate_eps)
    sim.set_incident_angle(inc_ang=inc_ang,azi_ang=azi_ang)
    sim.add_layer(thickness=layer0_thickness,eps=layer0_eps)
    sim.solve_global_smatrix()
    txx.append(sim.S_parameters(orders=[0,0],direction='forward',port='transmission',polarization='xx',ref_order=[0,0]))
    rxx.append(sim.S_parameters(orders=[0,0],direction='forward',port='reflection',polarization='yy',ref_order=[0,0]))

txx = torch.cat(txx)
rxx = torch.cat(rxx)

T R

Blue lines are the results of TORCWA and red lines are the results of COMSOL simulation. As you can see, the values of COMSOL simulation and TORCWA are close, but slow convergence is seen as the Nx increases, and the fluctuation is greater than that of dielectric (can be found in the paper). Since the basic RCWA algorithm is known to be slow convergence, an algorithm capable of fast convergence for metallic structures have been presented such as factorization and adaptive spatial resolution (like reference https://opg.optica.org/josaa/fulltext.cfm?uri=josaa-3-11-1780). Algorithm additions and improvements are planned in the future, but are not currently implemented (In the current version, it is difficult to modify algorithm outside rcwa class). Therefore, when simulating a metallic structure, set Nxy according to the computer's available resources, but consider that the error may be larger than that of dielectric.

When putting a structure with add_layer, it doesn't matter if structure is in the form of numpy.array or torch.tensor of size [nx, ny]. However, the function that converts files such as GDSII into [nx, ny] array and tensor types is not supported now. Currently, it seems that you need to custom implement a function that converts into an array using KLayout, etc.

Thank you.

MahmutRuzi commented 1 year ago

Thank you for the detailed response. For the specific example I mentioned, it seems to start to converge for Nx=Ny=30. As expected, takes quick some time though.