kch3782 / torcwa

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

_LinAlgError for a specific lattice size (px=3000) #21

Closed zhixwang closed 1 year ago

zhixwang commented 1 year ago

Dear Authors,

Hi, thanks a lot for your library. I recently encountered a very strange error using this library. My code is as follows, which is only slightly modified from your Example0.ipynb:

# Import
import numpy as np
import torch
from matplotlib import pyplot as plt

import torcwa

# Hardware
# If GPU support TF32 tensor core, the matmul operation is faster than FP32 but with less precision.
# If you need accurate operation, you have to disable the flag below.
torch.backends.cuda.matmul.allow_tf32 = False
sim_dtype = torch.complex64
geo_dtype = torch.float32
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# Simulation environment
# light
lamb0 = 375.                # nm
azi_ang = 0.*(np.pi/180)    # radian

# material
substrate_eps =  1.4782**2

# geometry
px = 3000
L = [px, px*np.sqrt(3)]            # nm / nm
torcwa.rcwa_geo.device = device
torcwa.rcwa_geo.dtype = geo_dtype
torcwa.rcwa_geo.Lx = L[0]
torcwa.rcwa_geo.Ly = L[1]
res = 30
torcwa.rcwa_geo.nx = int(L[0]/res)
torcwa.rcwa_geo.ny = int(L[1]/res)
torcwa.rcwa_geo.grid()
# torcwa.rcwa_geo.edge_sharpness = 1000.

x_axis = torcwa.rcwa_geo.x.cpu()
y_axis = torcwa.rcwa_geo.y.cpu()

# Generate and perform simulation
order_N = 12
order = [order_N,order_N]
sim = torcwa.rcwa(freq=1/lamb0,order=order,L=L,dtype=sim_dtype,device=device)
sim.add_input_layer(eps=substrate_eps)
sim.add_output_layer(eps=1.**2)
sim.set_incident_angle(inc_ang=0.,azi_ang=0.)

sim.add_layer(thickness=300, eps=1.6**2)
sim.add_layer(thickness=80e3, eps=1.**2)

When I run this snippet (in a Jupyter notebook), I got the following error:

_LinAlgError                              Traceback (most recent call last)
[Example 0 - Fresnel equation_modify.ipynb Cell 1 in 5
     [50](line=49) sim.add_output_layer(eps=1.**2)
     [51](line=50) sim.set_incident_angle(inc_ang=0.,azi_ang=0.)
---> [53](line=52) sim.add_layer(thickness=300, eps=1.6**2)
     [54](line=53) sim.add_layer(thickness=80e3, eps=1.**2)

File [c:\...\anaconda3\envs\torch\lib\site-packages\torcwa\rcwa.py:170](file:///C:/.../anaconda3/envs/torch/lib/site-packages/torcwa/rcwa.py:170), in rcwa.add_layer(self, thickness, eps, mu)
    167 else:
    168     self._eigen_decomposition()
--> 170 self._solve_layer_smatrix()

File [c:\...\anaconda3\envs\torch\lib\site-packages\torcwa\rcwa.py:1271](file:///C:/.../anaconda3/envs/torch/lib/site-packages/torcwa/rcwa.py:1271), in rcwa._solve_layer_smatrix(self)
   1268 Ctmp = torch.hstack((Ctmp1,Ctmp2))
   1270 # Mode coupling coefficients
-> 1271 self.Cf.append(torch.matmul( torch.linalg.inv(Ctmp), torch.vstack((2*torch.eye(2*self.order_N,dtype=self._dtype,device=self._device),
   1272     torch.zeros([2*self.order_N,2*self.order_N],dtype=self._dtype,device=self._device))) ))
   1273 self.Cb.append(torch.matmul( torch.linalg.inv(Ctmp), torch.vstack((torch.zeros([2*self.order_N,2*self.order_N],dtype=self._dtype,device=self._device),
   1274     2*torch.eye(2*self.order_N,dtype=self._dtype,device=self._device))) ))
   1276 self.layer_S11.append(torch.matmul(torch.matmul(self.E_eigvec[-1],phase), self.Cf[-1][:2*self.order_N,:]) + torch.matmul(self.E_eigvec[-1],self.Cf[-1][2*self.order_N:,:]))

_LinAlgError: torch.linalg.inv: The diagonal element 2 is zero, the inversion could not be completed because the input matrix is singular.

I find this error very strange because it only happens when px=3000. If I change the value of px from 3000 to either 2999 or 3001, the error is no longer there.

Do you have any idea why it happens? Is there any wrong way of using your library here?

Many thanks!

Best wishes,

kch3782 commented 1 year ago

Hi.

This error occurs when the lattice period is an integer multiple of the wavelength (3000 = 375 * 8). (In this case, there are diffraction order component in a direction parallel to the layer. It makes singular matrix and error in calculation of inverse matrix.

A way to avoid this problem is to introduce a slight deviation in the wavelength (e.g., set lamb0 to 375.001).

Thank you.

zhixwang commented 1 year ago

Thanks a lot!