The differentiable objective quantities for the polarization components of the transmittance spectra: Txx, Txy

[Guizhixing]

Hi, all,
I am trying to use the adjoint solver to inverse design the unit cell structure of a metasurface, with the objective function constructed from the polarization components (Txx, Txy) of the transmittance spectra of planewave at normal incidence. Could you please advise me on how to obtain the differentiable objective quantities: Txx, Txy?

[smartalecH]

You can use an EigenmodeCoefficient to decompose the field into the different polarization states. And then formulate your FOM from that.

[Guizhixing]

Thank you very much for your reply! I am a new user of Meep, and over the past few days, I have created the three-dimensional model of metasurface accroding to the examples. I have placed the cross-sections of medel and the main code here. The region to be optimized is a 64 pixel * 64 pixel area, and I have initialized the epsilon values by reading them from local data. The FOM is the Mean Square Error between the differentiable objective spectra and the ground truth of transmittance spectra. However, after running the program, I did not obtain gradients with respect to the epsilon variables. Could you please help me identify the issues in this code? I mainly have three questions as following:
How to define the planewave with z-polarization at normal incidence？
How to define the periodic boundaries at y- and z- direction？
How to get the z-polarization or y-polarization component of transmission spectrum by EigenmodeCoefficient？

plot2D of model

Main code

#define the unit 
nm = 1e-3
C_0 = 1

# define the frequencies we want to measure
Nfreq = 40
wavelengths = nm * np.linspace(1050, 1400, Nfreq)
freqs = C_0 / wavelengths

# define the frequency center and width of our pulse
freq0 = freqs[len(freqs) // 2]
freqw = freqs[0] - freqs[-1]

# Define material properties
n_SiO2 = 1.46
n_Si = 3.52
SiO2 = mp.Medium(epsilon=n_SiO2**2)
Si = mp.Medium(epsilon=n_Si**2)
Air = mp.Medium(epsilon=1)

resolution = 100  # pixels/um

d_pml = 250*nm  # PML thickness

# geometric parameters
P = 630 * nm  # periodicity in x and y
d_Si = 260 * nm  # thickness of silcon
d_SiO2 = 2 * d_Si  # substrate thickness
d_pad = 0.5 * P  # padding between metasurface and PML

sy = sz = P
sx = d_pml + d_pad + d_Si + d_SiO2 + 2*d_pad + d_pml

cell_size = mp.Vector3(sx, sy, sz)
pml_layers = [mp.PML(thickness=d_pml, direction=mp.X)]

minimum_length = 2.5*nm  # minimum length scale 
eta_i = (
    0.5  # blueprint (or intermediate) design field thresholding point (between 0 and 1)
)
eta_e = 0.55  # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e  # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
print(filter_radius)
def mapping(x, eta, beta):

    # filter
    filtered_field = mpa.conic_filter(
        x,
        filter_radius,
        sy,
        sz,
        resolution,
    )

    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)

    projected_field = (
        npa.rot90(projected_field.T, 2)
    ) / 2  # 90-degree rotation

    # interpolate to actual materials
    return projected_field.flatten()

# define the design region and design variables
Ny = int(resolution * sy) + 1
Nz = int(resolution * sz) + 1

design_variables = mp.MaterialGrid(mp.Vector3(0, Ny, Nz), Air, Si, grid_type="U_MEAN")
design_region = mpa.DesignRegion(
    design_variables,
    volume=mp.Volume(
        center=mp.Vector3(d_Si/2,0,0),
        size=mp.Vector3(d_Si, sy, sz),
    ),
)

geometry = [
    mp.Block(
        material=SiO2,
        size=mp.Vector3(d_SiO2, sy, sz),
        center=mp.Vector3(-d_SiO2/2,0,0),
    ),
    mp.Block(
        center=design_region.center, 
        size=design_region.size, 
        material=design_variables
    ),
]

#How to define the planewave with z-polarization at normal incidence   <-------------Question 1
src_pt = mp.Vector3(-d_SiO2 - d_pad/4, 0, 0) # source center
sources = [
    mp.Source(
        mp.GaussianSource(freq0, fwidth=freqw),
        component=mp.Ez,
        center=src_pt,
        size=mp.Vector3(x=0, y=sy, z=sz),
    )
]
###############################################################################################

# First simulation without the geometry to obtain the input flux
k_point = mp.Vector3(0, 0, 0)
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    boundary_layers=pml_layers, ###  How to define the periodic boundaries at y- and z- direction   <-------------Question 2
    k_point=k_point,            ###  How to define the periodic boundaries at y- and z- direction   <-------------Question 2
    sources=sources,
)

mon_pt = mp.Vector3(d_Si + d_pad/2, 0, 0) # center of monitor
tran_mon = sim.add_mode_monitor(
        freqs,
        mp.ModeRegion(center=mon_pt, size=mp.Vector3(x=0, y=sy, z=sz)),
        yee_grid=True,
    )

# stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-8)
sim.run(until_after_sources=1)

# How to get the z-polarization component of transmission spectrum by EigenmodeCoefficient  <-------------Question 3
eig_parity = mp.ODD_Z
res = sim.get_eigenmode_coefficients(
        tran_mon,
        [1],
        eig_parity=eig_parity,
)

coeffs = res.alpha
input_flux = np.abs(coeffs[0, :, 0]) ** 2
print(input_flux)

sim.reset_meep()
sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    geometry=geometry,
    boundary_layers=pml_layers, ###  How to define the periodic boundaries at y- and z- direction   <-------------Question 2
    k_point=k_point,            ###  How to define the periodic boundaries at y- and z- direction   <-------------Question 2
#     eps_averaging=False,
    sources=sources,
)

# How to get the z-polarization component of transmission spectrum by EigenmodeCoefficient  <-------------Question 3
mode = 1
TE = mpa.EigenmodeCoefficient(
    sim,
    mp.Volume(
        center=mon_pt,
        size=mp.Vector3(x=0, y=sy, z=sz),
    ),
    mode,
    eig_parity=eig_parity,
)
ob_list = [TE]
########################################################################

# read ground truth of transmittance spectra from csv file
import pandas as pd
Tzz = pd.read_csv(r'home\hobart\Meep\Tzz.csv').values[0,0:40]

# The FOM is the Mean Square Error between the differentiable objective quantities and the ground truth of transmittance spectra
def J(alpha):
    Tzz = (npa.abs(alpha) ** 2/input_flux - Tzz)**2
    return npa.mean(Tzz)

opt = mpa.OptimizationProblem(
    simulation=sim,
    objective_functions=J,
    objective_arguments=ob_list,
    design_regions=[design_region],
    frequencies=freqs
)

f0, dJ_du = opt([mapping(eps_vector, eta_i, 1e3)])
print(dJ_du.shape)

[Guizhixing]

Unfortunately, I didn't get the right Txx and Txy components and I'm no longer using MEEP.

…

-----原始邮件----- 发件人:ramji28 ***@***.***> 发送时间:2024-04-20 23:06:12 (星期六) 收件人: NanoComp/meep ***@***.***> 抄送: Hobart ***@***.***>, Author ***@***.***> 主题: Re: [NanoComp/meep] The differentiable objective quantities for the polarization components of the transmittance spectra: Txx, Txy (Discussion #2781) Hello, were you able to get the Txx and Txy components? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>