repository update. Added metaprop_fibermode_doublet_sym.py which implements a constrained optimization enforcing mirror symmetry about the x and y axes. Mirror symmetry can be exploited in this problem since the target mode has such symmetry. Each metasurface has effectively a quarter of the degrees of freedom compared to the unconstrained problem.

2026-10-14 23:06:39 +00:00 · 2026-08-19 17:48:11 +02:00
parent 5ca4f4045a
commit e63c11ea58
37 changed files with 11478 additions and 40 deletions
@@ -0,0 +1,25 @@
+% Model:              cyljnderUnjtCell_solved.mph
+% Versjon:            COMSOL 6.1.0.282
+% Date:               Aug 7 2026, 11:04
+% Table:              Table 3 - Global Evaluatjon 3
+% tPost (nm)                                      R (nm)                                            lambda0 (µm)                                      S21 (1)
+900                                               50                                                1.55                                              -0.845669386929+0.5044362611307165j
+900                                               63.157894736842096                                1.55                                              -0.8251549977113691+0.5387569255954204j
+900                                               76.31578947368419                                 1.55                                              -0.7941914472021001+0.5854009666643545j
+900                                               89.47368421052632                                 1.55                                              -0.7471355423796626+0.6469515347910716j
+900                                               102.6315789473684                                 1.55                                              -0.6739731961216431+0.7260451361240156j
+900                                               115.78947368421052                                1.55                                              -0.5559252634519462+0.8236596878659627j
+900                                               128.94736842105263                                1.55                                              -0.35623973597781955+0.9314716653309569j
+900                                               142.1052631578947                                 1.55                                              -0.005446778866258663+0.9995911107511558j
+900                                               155.26315789473682                                1.55                                              0.5624693373844357+0.8192539801083378j
+900                                               168.42105263157896                                1.55                                              0.9716453344468975-0.026879486169522106j
+900                                               181.57894736842104                                1.55                                              0.33229983484123915-0.9109913324745208j
+900                                               194.73684210526312                                1.55                                              -0.5138804677648452-0.8375112454139356j
+900                                               207.89473684210526                                1.55                                              -0.8979812315659875-0.41599951560137755j
+900                                               221.05263157894737                                1.55                                              -0.8666228661809928-0.35839703156981834j
+900                                               234.21052631578942                                1.55                                              -0.9706534677777011+0.1672691323977627j
+900                                               247.36842105263153                                1.55                                              -0.875683679549086+0.4301009457755804j
+900                                               260.52631578947364                                1.55                                              -0.7186563387658285+0.6394472785795947j
+900                                               273.6842105263158                                 1.55                                              -0.6152456609406248+0.5634675095025711j
+900                                               286.8421052631579                                 1.55                                              -0.3669553952916399+0.7919528013254294j
+900                                               300                                               1.55                                              -0.07785011501821787+0.7863645334693915j
@@ -0,0 +1,31 @@
+nohup: ignoring input
+Rank 1 is processing diameter 0.274 (index 8)
+Rank 2 is processing diameter 0.448 (index 16)
+Rank 0 is processing diameter 0.100 (index 0)
+Rank 0 is processing diameter 0.122 (index 1)
+Rank 1 is processing diameter 0.296 (index 9)
+Rank 0 is processing diameter 0.143 (index 2)
+Rank 1 is processing diameter 0.317 (index 10)
+Rank 2 is processing diameter 0.470 (index 17)
+Rank 0 is processing diameter 0.165 (index 3)
+Rank 1 is processing diameter 0.339 (index 11)
+Rank 0 is processing diameter 0.187 (index 4)
+Rank 1 is processing diameter 0.361 (index 12)
+Rank 2 is processing diameter 0.491 (index 18)
+Rank 0 is processing diameter 0.209 (index 5)
+Rank 1 is processing diameter 0.383 (index 13)
+Rank 2 is processing diameter 0.513 (index 19)
+Rank 0 is processing diameter 0.230 (index 6)
+Rank 1 is processing diameter 0.404 (index 14)
+Rank 2 is processing diameter 0.535 (index 20)
+Rank 0 is processing diameter 0.252 (index 7)
+Rank 1 is processing diameter 0.426 (index 15)
+Rank 2 is processing diameter 0.557 (index 21)
+Rank 2 is processing diameter 0.578 (index 22)
+Rank 2 is processing diameter 0.600 (index 23)
+
+Elapsed run time = 61296.2745 s
+
+Elapsed run time = 61296.2909 s
+
+Elapsed run time = 61299.1873 s
@@ -0,0 +1,31 @@
+nohup: ignoring input
+Rank 2: Running simulation for diameter 0.4478260869565217
+Rank 1: Running simulation for diameter 0.27391304347826084
+Rank 0: Running simulation for diameter 0.1
+Rank 0: Running simulation for diameter 0.12173913043478261
+Rank 1: Running simulation for diameter 0.29565217391304344
+Rank 0: Running simulation for diameter 0.14347826086956522
+Rank 1: Running simulation for diameter 0.3173913043478261
+Rank 2: Running simulation for diameter 0.4695652173913043
+Rank 0: Running simulation for diameter 0.16521739130434782
+Rank 1: Running simulation for diameter 0.33913043478260874
+Rank 0: Running simulation for diameter 0.18695652173913044
+Rank 1: Running simulation for diameter 0.36086956521739133
+Rank 2: Running simulation for diameter 0.4913043478260869
+Rank 0: Running simulation for diameter 0.20869565217391306
+Rank 1: Running simulation for diameter 0.3826086956521739
+Rank 2: Running simulation for diameter 0.5130434782608696
+Rank 0: Running simulation for diameter 0.23043478260869565
+Rank 1: Running simulation for diameter 0.4043478260869565
+Rank 2: Running simulation for diameter 0.5347826086956522
+Rank 0: Running simulation for diameter 0.25217391304347825
+Rank 1: Running simulation for diameter 0.4260869565217391
+Rank 2: Running simulation for diameter 0.5565217391304348
+
+Elapsed run time = 5816.1358 s
+Rank 2: Running simulation for diameter 0.5782608695652174
+Rank 2: Running simulation for diameter 0.6
+
+Elapsed run time = 10194.4932 s
+
+Elapsed run time = 10198.0417 s
@@ -0,0 +1,169 @@
+import meep as mp
+mp.verbosity(0)
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+resolution = 80 # pixels/unit length
+
+wvl_min = 1.5
+wvl_max = 1.6
+wvl_cen = 0.5*(wvl_min+wvl_max)
+fmin = 1/wvl_max
+fmax = 1/wvl_min
+fcen = 1/wvl_cen
+df = fmax-fmin
+nfreq = 1
+
+dpml = 1.0 #PML thickness
+dsub = 1.5 + wvl_cen #substrate thickness
+dpad = dsub #padding between cell and PML
+
+k_point = mp.Vector3(0,0,0) #used for PBC
+
+n_si = 3.45 #silicon refractive index
+si = mp.Medium(index=n_si)
+n_sio2 = 1.45
+sio2 = mp.Medium(index=n_sio2)
+n_air = 1.0 #air refractive index
+air = mp.Medium(index=n_air)
+
+P = 0.7 #unit cell pitch
+h = 0.9 #pillar height
+
+def unitCell_paramsweep(D):
+    sim_x_size = dpml + dsub + dpad + dpml
+    sim_z_size = P
+    sim_y_size = P
+
+    cell_size = mp.Vector3(sim_x_size, sim_y_size, sim_z_size)
+
+    source_position = mp.Vector3(-wvl_cen,0,0) #source position in the middle of the substrate
+    source = mp.Source(
+        mp.GaussianSource(fcen, fwidth=df),
+        component = mp.Ez,
+        center = source_position,
+        size = mp.Vector3(0, sim_y_size, sim_z_size),
+    )
+
+    pml_layers = [mp.PML(dpml,direction=mp.X)]
+
+    monitor_position = mp.Vector3(h + wvl_cen,0,0)
+    obs_point = monitor_position
+    stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, obs_point, 1e-9)
+    #stop_cond = 10
+    sim = mp.Simulation(
+        resolution = resolution,
+        cell_size = cell_size,
+        boundary_layers = pml_layers,
+        k_point = k_point,
+        default_material = sio2,
+        sources = [source],
+    )
+
+    flux_mon = sim.add_flux(
+        [fcen],
+        mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size))
+    )
+
+    sim.run(until_after_sources=stop_cond)
+    input_flux = mp.get_fluxes(flux_mon)
+
+    sim.reset_meep()
+
+    pillar_geo = mp.Cylinder(
+        material = si,
+        radius = D/2,
+        height = h,
+        center = mp.Vector3(h/2,0,0),
+        axis=mp.Vector3(1,0,0),
+    )
+
+    substrate_geo = mp.Block(
+        material = sio2,
+        size = mp.Vector3(dsub+dpml, mp.inf, mp.inf),
+        center = mp.Vector3(-(dsub+dpml)/2,0,0),
+    )
+
+    geometries = [substrate_geo, pillar_geo]
+
+    sim = mp.Simulation(
+        resolution = resolution,
+        cell_size = cell_size,
+        boundary_layers = pml_layers,
+        geometry = geometries,
+        k_point = k_point,
+        sources = [source],
+    )
+
+    mode_mon = sim.add_flux(
+        [fcen],
+        mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size))
+    )
+
+    sim.run(until_after_sources=stop_cond)
+
+    freqs = mp.get_eigenmode_freqs(mode_mon)
+    res = sim.get_eigenmode_coefficients(
+        mode_mon, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y
+    )
+    coeffs = res.alpha
+
+    mode_wvl = [1/freqs[nf] for nf in range(nfreq)]
+    mode_tran = [abs(coeffs[0,nf,0])**2 / input_flux[nf] for nf in range(nfreq)]
+    mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)]
+
+    return mode_wvl, mode_tran, mode_phase
+
+if rank == 0:
+    diameters = np.linspace(0.1, 0.6, 24)
+else:
+    diameters = None
+diameters = comm.bcast(diameters, root=0)
+
+local_diameters = np.array_split(diameters, size)[rank]
+local_indices = np.array_split(np.arange(len(diameters)), size)[rank]
+
+local_mode_tran_arr = np.zeros(len(diameters))
+local_mode_phase_arr = np.zeros(len(diameters))
+
+for local_i, (global_i, D) in enumerate(zip(local_indices, local_diameters)):
+    print(f"Rank {rank} is processing diameter {D:.3f} (index {global_i})")
+    mode_wvl, mode_tran, mode_phase = unitCell_paramsweep(D)
+    local_mode_tran_arr[global_i] = mode_tran[0]
+    local_mode_phase_arr[global_i] = mode_phase[0]
+
+comm.Barrier()
+mode_tran_arr = comm.gather(local_mode_tran_arr, root=0)
+mode_phase_arr = comm.gather(local_mode_phase_arr, root=0)
+final_mode_tran_arr = np.zeros(len(diameters))
+final_mode_phase_arr = np.zeros(len(diameters))
+if rank == 0:
+    for i in range(size):
+        final_mode_tran_arr += mode_tran_arr[i]
+        final_mode_phase_arr += mode_phase_arr[i]
+    
+    
+    np.save("unit_cell_center_wvl_diameters.npy", diameters)
+    np.save("unit_cell_center_wvl_mode_wvl.npy", mode_wvl)
+    np.save("unit_cell_center_wvl_mode_tran.npy", final_mode_tran_arr)
+    np.save("unit_cell_center_wvl_mode_phase.npy", final_mode_phase_arr)
+
+    plt.figure()
+    plt.subplot(1,2,1)
+    plt.plot(diameters,final_mode_tran_arr)
+    plt.xlabel("Pillar Diameter (um)")
+    plt.ylabel("Transmission")
+    plt.subplot(1,2,2)
+    plt.plot(diameters,final_mode_phase_arr)
+    plt.xlabel("Pillar Diameter (um)")
+    plt.ylabel("Phase (rad)")
+    plt.subplots_adjust(wspace=0.5)
+    plt.savefig("unit_cell_fcen.pdf")
+
@@ -0,0 +1,203 @@
+import meep as mp
+mp.verbosity(0)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+resolution = 50 # pixels/unit length
+
+wvl_min = 1.5
+wvl_max = 1.6
+wvl_cen = 0.5*(wvl_min+wvl_max)
+fmin = 1/wvl_max
+fmax = 1/wvl_min
+fcen = 1/wvl_cen
+df = fmax-fmin
+nfreq = 20
+
+dpml = 1.0 #PML thickness
+dsub = 1.5 + wvl_cen #substrate thickness
+dpad = dsub #padding between cell and PML
+
+k_point = mp.Vector3(0,0,0) #used for PBC
+
+n_si = 3.45 #silicon refractive index
+si = mp.Medium(index=n_si)
+n_sio2 = 1.45
+sio2 = mp.Medium(index=n_sio2)
+n_air = 1.0 #air refractive index
+air = mp.Medium(index=n_air)
+
+P = 0.7 #unit cell pitch
+h = 0.9 #pillar height
+
+def unitCell_paramsweep(D):
+    print("Rank {}: Running simulation for diameter {}".format(rank, D))
+    sim_x_size = dpml + dsub + dpad + dpml
+    sim_z_size = P
+    sim_y_size = P
+
+    cell_size = mp.Vector3(sim_x_size, sim_y_size, sim_z_size)
+
+    source_position = mp.Vector3(-wvl_cen,0,0) #source position in the middle of the substrate
+    source = mp.Source(
+        mp.GaussianSource(fcen, fwidth=df),
+        component = mp.Ez,
+        center = source_position,
+        size = mp.Vector3(0, sim_y_size, sim_z_size),
+    )
+
+    pml_layers = [mp.PML(dpml,direction=mp.X)]
+
+    monitor_position = mp.Vector3(h + wvl_cen,0,0)
+    obs_point = monitor_position
+    stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, obs_point, 1e-9)
+    #stop_cond = 10
+    sim = mp.Simulation(
+        resolution = resolution,
+        cell_size = cell_size,
+        boundary_layers = pml_layers,
+        k_point = k_point,
+        default_material = sio2,
+        sources = [source],
+    )
+
+    flux_mon = sim.add_flux(
+        fcen, 
+        df,
+        nfreq, 
+        mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size))
+    )
+
+    sim.run(until_after_sources=stop_cond)
+    input_flux = mp.get_fluxes(flux_mon)
+
+    sim.reset_meep()
+
+    pillar_geo = mp.Cylinder(
+        material = si,
+        radius = D/2,
+        height = h,
+        center = mp.Vector3(h/2,0,0),
+        axis=mp.Vector3(1,0,0),
+    )
+
+    substrate_geo = mp.Block(
+        material = sio2,
+        size = mp.Vector3(dsub+dpml, mp.inf, mp.inf),
+        center = mp.Vector3(-(dsub+dpml)/2,0,0),
+    )
+
+    geometries = [substrate_geo, pillar_geo]
+
+    sim = mp.Simulation(
+        resolution = resolution,
+        cell_size = cell_size,
+        boundary_layers = pml_layers,
+        geometry = geometries,
+        k_point = k_point,
+        sources = [source],
+    )
+
+    mode_mon = sim.add_flux(
+        fcen, df, nfreq, mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size))
+    )
+
+    sim.run(until_after_sources=stop_cond)
+
+    freqs = mp.get_eigenmode_freqs(mode_mon)
+    res = sim.get_eigenmode_coefficients(
+        mode_mon, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y
+    )
+    coeffs = res.alpha
+
+    mode_wvl = [1/freqs[nf] for nf in range(nfreq)]
+    mode_tran = [abs(coeffs[0,nf,0])**2 / input_flux[nf] for nf in range(nfreq)]
+    mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)]
+
+    return mode_wvl, mode_tran, mode_phase
+
+if rank == 0:
+    diameters = np.linspace(0.1, 0.6, 24)
+else:
+    diameters = None
+diameters = comm.bcast(diameters, root=0)
+local_diameters = np.array_split(diameters, size)[rank]
+local_indices = np.array_split(np.arange(len(diameters)), size)[rank]
+
+local_mode_tran_arr = np.zeros((len(diameters), nfreq))
+local_mode_phase_arr = np.zeros((len(diameters), nfreq))
+
+for local_i, (global_i, D) in enumerate(zip(local_indices,local_diameters)):
+    mode_wvl, mode_tran, mode_phase = unitCell_paramsweep(D)
+    local_mode_tran_arr[global_i,:] = mode_tran
+    local_mode_phase_arr[global_i,:] = mode_phase
+
+
+mode_tran_arr = comm.gather(local_mode_tran_arr, root=0)
+mode_phase_arr = comm.gather(local_mode_phase_arr, root=0)
+
+if rank == 0:
+    final_mode_tran_arr = np.zeros((len(diameters), nfreq))
+    final_mode_phase_arr = np.zeros((len(diameters), nfreq))
+    num_diameters_per_rank = len(diameters) // size
+    for i in range(0, size):
+        for j in range(num_diameters_per_rank):
+           global_index = i * num_diameters_per_rank + j
+           final_mode_tran_arr[global_index, :] = mode_tran_arr[i][global_index, :]
+           final_mode_phase_arr[global_index, :] = mode_phase_arr[i][global_index, :]
+
+
+    np.save("unit_cell_param_sweep_diameters.npy", diameters)
+    np.save("unit_cell_param_sweep_mode_wvl.npy", mode_wvl)
+    np.save("unit_cell_param_sweep_mode_tran.npy", final_mode_tran_arr)
+    np.save("unit_cell_param_sweep_mode_phase.npy", final_mode_phase_arr)
+
+    plt.figure()
+    plt.subplot(1,2,1)
+    plt.pcolormesh(
+        mode_wvl,
+        diameters,
+        final_mode_tran_arr,
+        cmap='hot_r',
+        shading='gouraud',
+        vmin=0,
+        vmax=final_mode_tran_arr.max(),
+    )
+
+    plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]])
+    plt.xlabel("Wavelength (um)")
+    plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)])
+    plt.ylabel("Pillar Diameter (um)")
+    plt.title("Transmission")
+    cbar=plt.colorbar()
+    cbar.set_ticks([t for t in np.arange(0, 1.2, 0.2)])
+
+    plt.subplot(1,2,2)
+    plt.pcolormesh(
+        mode_wvl,
+        diameters,
+        final_mode_phase_arr,
+        cmap='RdBu',
+        shading='gouraud',
+        vmin=final_mode_phase_arr.min(),
+        vmax=final_mode_phase_arr.max(),
+    )
+
+    plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]])
+    plt.xlabel("Wavelength (um)")
+    plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)])
+    plt.ylabel("Pillar Diameter (um)")
+    plt.title("phase")
+    cbar=plt.colorbar()
+    cbar.set_ticks([t for t in np.arange(-3, 4, 1)])
+
+
+    plt.subplots_adjust(wspace=0.5)
+    plt.savefig("unit_cell_param_sweep.pdf")
+
@@ -0,0 +1,41 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+#import the parameter sweep arrays for the central wavelength
+diameters = np.load("unit_cell_center_1550_h900nm_res50_wvl_diameters.npy")
+mode_wvl = np.load("unit_cell_center_1550_h900nm_res50_wvl_mode_wvl.npy")
+mode_tran_arr = np.load("unit_cell_center_1550_h900nm_res50_wvl_mode_tran.npy")
+print(mode_tran_arr)
+mode_phase_arr = np.load("unit_cell_center_1550_h900nm_res50_wvl_mode_phase.npy")
+print(mode_wvl)
+unwrap_mode_phase_arr = np.abs(np.unwrap(mode_phase_arr - mode_phase_arr[0], period=2*np.pi))
+#import the reference value array
+ref_results = np.loadtxt("results_S21_p700nm_h900nm.txt", skiprows=5,dtype=complex)
+R_ref = np.real(ref_results[:,1])
+S21_ref = ref_results[:,3]
+S21_ref_mag = np.abs(S21_ref)
+S21_ref_phase = np.abs(np.unwrap(np.angle(S21_ref) - np.angle(S21_ref[0]), period=2*np.pi))
+
+
+
+plt.figure()
+plt.subplot(1,2,1)
+plt.plot(1e3*diameters/2,mode_tran_arr, '^r-.', label='FDTD')
+plt.plot(R_ref,S21_ref_mag, '*k--',label='FEM')
+plt.xlabel("Pillar Radius (um)")
+plt.xticks([t for t in np.linspace(1e3*diameters[0]/2, 1e3*diameters[-1]/2, 5)])
+plt.gca().set_xticklabels([f"{t:.0f}" for t in np.linspace(1e3*diameters[0]/2, 1e3*diameters[-1]/2, 5)])
+plt.ylabel("Transmission")
+plt.legend()
+
+plt.subplot(1,2,2)
+plt.plot(1e3*diameters/2,unwrap_mode_phase_arr, '^r-.', label='FDTD')
+plt.plot(R_ref,S21_ref_phase, '*k--', label='FEM')
+plt.xlabel("Pillar Radius (um)")
+plt.xticks([t for t in np.linspace(1e3*diameters[0]/2, 1e3*diameters[-1]/2, 5)])
+plt.gca().set_xticklabels([f"{t:.0f}" for t in np.linspace(1e3*diameters[0]/2, 1e3*diameters[-1]/2, 5)])
+plt.ylabel("Phase (rad)")
+plt.legend()
+
+plt.subplots_adjust(wspace=0.5)
+plt.savefig("unit_cell_param_sweep_lineplot.pdf")