diff --git a/hdf5_vis.py b/hdf5_vis.py new file mode 100644 index 0000000..eda9917 --- /dev/null +++ b/hdf5_vis.py @@ -0,0 +1,49 @@ +#visualize fields from hdf5 files +import h5py +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import glob +import numpy as np + +input_files = sorted(glob.glob("metapeg_v2-ex-*.h5")) + +data = [] +times = [] + +for filename in input_files: + with h5py.File(filename, "r") as f: + Ex = f["ex"][:] + # Extract time from filename (e.g., "metapeg_v2-ex-000050.00.h5" → 50.00) + time = float(filename.split("-")[-1].split(".h5")[0]) + times.append(time) + data.append(Ex) + print(f"Loaded {filename}, time={time}") +print(len(data)) +# Create animation +fig, ax = plt.subplots() +x = np.linspace(-4, 4) #assuming x is the last dimension +z = np.linspace(-3, 3) #assuming z is the first dimension +max_clim = np.max([np.abs(data[i][:,0,:]) for i in range(len(data))]) +min_clim = np.min([np.abs(data[i][:,0,:]) for i in range(len(data))]) +print(max_clim) +def animate(frame): + ax.clear() + + im = ax.imshow(np.rot90(np.abs(data[frame][:,0,:]),k=1), origin="lower", + extent=[z.min(), z.max(),x.min(), x.max(),], + + cmap="viridis", + vmin=min_clim, vmax=max_clim) + + ax.set_xlabel("x (um)") + ax.set_ylabel("z (um)") + ax.set_title(f"Electric field (Ex) at t={times[frame]:.2f}") + #plt.colorbar(im, ax=ax, label="Ex (a.u.)") + return im, + +anim = animation.FuncAnimation(fig, animate, frames=len(data), + interval=100, blit=True, repeat=True) +anim.save("results/Ex_animation.gif", writer="ffmpeg", fps=2) +plt.close() +print("Animation saved to results/Ex_animation.gif") + diff --git a/metapeg_v2.py b/metapeg_v2.py index c6395aa..e8e2652 100644 --- a/metapeg_v2.py +++ b/metapeg_v2.py @@ -2,6 +2,7 @@ from typing import List, NamedTuple, Tuple from autograd import numpy as npa, tensor_jacobian_product, grad import matplotlib.pyplot as plt +import matplotlib.animation as animation import meep as mp import meep.adjoint as mpa import nlopt @@ -10,25 +11,25 @@ import os import xarray as xr from PIL import Image -mp.verbosity(1) +mp.verbosity(3) if not os.path.exists("results"): os.makedirs("results") ###### PHYSICS SETUP ###### -RESOLUTION = 30 #50 #pixels per unit length +RESOLUTION = 18 #pixels per meep unit length (1um) MAX_RUN_TIME = 1100 #meep units WAVELENGTH_MIN_UM = 1.50 WAVELENGTH_MAX_UM = 1.60 PML_UM = 1.0 #PML thickness -BUFFER = 3.0 +BUFFER = 2.0 DESIGN_WAVELENGTHS_UM = [1.50,1.55,1.60] #wavelengths at which to optimize the device -DESIGN_REGION_UM = mp.Vector3(5,5,1) +DESIGN_REGION_UM = mp.Vector3(2,2,1) DESIGN_REGION_RESOLUTION = int(2 * RESOLUTION) DESIGN_REGION_CENTER = mp.Vector3(0, 0, DESIGN_REGION_UM.z / 2) NX_DESIGN_GRID = int(DESIGN_REGION_UM.x * DESIGN_REGION_RESOLUTION) + 1 NY_DESIGN_GRID = int(DESIGN_REGION_UM.y * DESIGN_REGION_RESOLUTION) + 1 -MIN_LENGTH_UM = 0.15 #150 nm minimum length scale +MIN_LENGTH_UM = 0.5 #500 nm minimum length scale SILICON = mp.Medium(index = 3.5) SILICON_DIOXIDE = mp.Medium(index = 1.5) AIR = mp.Medium(index = 1.0) @@ -63,10 +64,10 @@ frequency_center = 0.5 * (frequency_min + frequency_max) wavelength_center = 1 / frequency_center #remember that c = 1 in meep units frequency_width = frequency_max - frequency_min -SOURCE_CENTER = mp.Vector3(0,0,-2) +SOURCE_CENTER = mp.Vector3(0,0,-BUFFER/2) SOURCE_SIZE = mp.Vector3(physical_domains_size.x,physical_domains_size.y,0) -NEAR_REGION_MONITOR_CENTER = mp.Vector3(0,0, DESIGN_REGION_UM.z + BUFFER/2) +NEAR_REGION_MONITOR_CENTER = mp.Vector3(0,0, DESIGN_REGION_UM.z) NEAR_REGION_MONITOR_SIZE = mp.Vector3(physical_domains_size.x,physical_domains_size.y,0) ff_monitor_center = mp.Vector3(0,0,DESIGN_REGION_UM.z + BUFFER) @@ -75,7 +76,7 @@ xs = ys = np.linspace(-ff_monitor_size.x/2, ff_monitor_size.x/2, NX_DESIGN_GRID) #collection of points in the far-field monitor ff_points = [mp.Vector3(x_p,y_p,ff_monitor_center.z) for x_p in xs for y_p in ys] -stop_cond = mp.stop_when_fields_decayed(50, mp.Ex, NEAR_REGION_MONITOR_CENTER, 1e-6) +stop_cond = mp.stop_when_fields_decayed(25, mp.Ex, NEAR_REGION_MONITOR_CENTER, 1e-6) pml_layers = [mp.PML(PML_UM)] frequencies = [1/wavelength for wavelength in DESIGN_WAVELENGTHS_UM] @@ -343,6 +344,15 @@ def intensity_desired_fn_pattern( return pattern_to_return + +def record_fields(sim): + plt.figure() + ax = plt.gca() + sim.plot2D(output_plane=view_2D_plane, fields=mp.Ex, ax=ax) + plt.savefig(f"results/fields_t{sim.timestep()}.pdf") + plt.close() + return True + def normalization_sim() -> np.ndarray: ''' Computes the Far-Field pattern at the monitor location and without the design region. Used for normalization purposes. @@ -354,7 +364,9 @@ def normalization_sim() -> np.ndarray: sources = [ mp.Source( - src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + src=mp.GaussianSource( + frequency=frequency_center, + fwidth=frequency_width/2), component=mp.Ex, center=SOURCE_CENTER, size=SOURCE_SIZE, @@ -379,6 +391,7 @@ def normalization_sim() -> np.ndarray: geometry=geometry, boundary_layers=pml_layers, k_point=mp.Vector3(), + split_chunks_evenly=False, ) plt.figure() @@ -395,7 +408,21 @@ def normalization_sim() -> np.ndarray: norm_near2far = norm_sim.add_near2far(frequencies, *NearRegions) - norm_sim.run(until_after_sources=stop_cond) + #norm_sim.run(until_after_sources=stop_cond) + norm_sim.run( + mp.at_every(20,mp.in_volume( + mp.Volume(center=mp.Vector3(), size=physical_domains_size), + mp.output_efield_x + )), + #mp.at_every(20,record_fields), + until_after_sources=stop_cond + ) + #make a plot of the fields + plt.figure() + ax = plt.gca() + norm_sim.plot2D(output_plane=view_2D_plane, fields=mp.Ex, ax=ax) + plt.savefig("results/NORM_FIELDS.pdf") + ref_fields = np.array( [norm_sim.get_farfield(norm_near2far, point) for point in ff_points] @@ -507,6 +534,7 @@ def intensity_optimization( geometry=geometry, boundary_layers=pml_layers, k_point=mp.Vector3(), + split_chunks_evenly=False, ) plt.figure()