Unexpected RuntimeError 3D Near-to-Far Fields with a k-point complains of "dft_near2far requires surfaces in a homogeneous medium"

[gummyequation]

Unexpected RuntimeError when computing Near-to-Far Fields

I am simulating a square lattice nanowire in 3-D using MEEP, using the tutorial “Near-to-Far Field Spectra” as a reference. I’ve defined the parameters for a square lattice nanowire, such as the geometry, PMLs, source, and k-point, near2far_box, and finally the computation to get the far-fields without a k-point.

I have successfully obtained the near fields with DFT and PadeDFT methods, and the plots agree with simulations done in COMSOL with the same simulation properties.

The near-to-far-fields work as expected without a k point, (besides vispy render) however with a k-point added to sim.Simulation( ) this error emerges:“RuntimeError: meep: dft_near2far requires surfaces in a homogeneous medium” .

Unless I understand the contents of the tutorial and the information in the documentation, or I am badly implementing them the error does not make sense. The error seems to stem from the “ _meep.fields_add_dft_near2far(self, *args)” — as the traceback suggests — yet the parameters of the near2far_box are what they should be.

Visualization with VisPy has a few issues. The first being the apparent orientation of the near2far planes. Additionally, the direction arguments of the near2far object change nothing.

Either vispy provides accurate information of the parameters of the simulation or, there are issues with how vispy renders the provided data.

I find it difficult to learn how to obtain the fields far from the near field with my parameters. I am trying to solve for the fields and not the flux and the documentation is a bit vague about computing the far fields. I understand that the near to far box needs to be a in homogenous medium, the expressions for coordinates satisfy that condition, yet the runtimerror still occurs.

An additional note on vispy, the second issue is when the dpml is increased beyond 0.3 the unit cell expands as well, even if the thickness of the PML only changes in one direction.

Please let me know if anyone has run into a similar issue or has any advice to give.

import meep as mp
import numpy as np
import vispy

Define simulation parameters
resolution = 25 # [pixels/um]

# PML
dpml = 0.200 #thickness of pml
pml_layers = [mp.PML(dpml, direction=mp.Z)] # PML only on the top and bottom of the Nanorod

# Define the nanowire
nw_d = 0.173
nw_r = nw_d/2 # radius of the nanowire
height_nw = 0.500
eps = 5.35

# Create the nanowire geometry
geometry = [mp.Cylinder(radius=nw_r,height = height_nw, material=mp.Medium(epsilon=eps),center=mp.Vector3(0,0,0))]

# Cell dimensions
a = 0.200
sy = a
sx = a
gap = 0.400 # equivalent to a maximum wavelength
sz = height_nw + 2*dpml + 2*gap # functional dependence of sz
sz = np.round(sz,4)
cell_size = mp.Vector3(sx,sy,sz)

# Define the source
fcen = 1.00 # Pulse center frequency
df = 3.00 # Pulse width
# Source
source = [
    mp.Source(mp.GaussianSource(fcen, fwidth=df),
    component=mp.Ez,center=mp.Vector3(0, nw_r/6, 0))] # a bit off-center

# Near-field box surrounding the nanowire in homogeneous medium
delta = 0.001 # Small shift to place near2far surfaces just outside the PML and nanowire (but within the cell)

# The height of the cell minus that should be occupied by the near2box
height_to_cover = sz/2 + 2*dpml + 2*gap - delta
distance_from_top = sz/2 - dpml/2 - gap/2 + delta

# Region between the unit cell and cell of the simulation
distance_from_center = np.round(sx/2 - nw_r +delta,4)

near2far_box =[
    # Z-side
    mp.Near2FarRegion(center=mp.Vector3(0, 0, distance_from_top),direction = mp.Y, size=mp.Vector3(sx, sy, 0), weight= 1 ),  # Top z side
    mp.Near2FarRegion(center=mp.Vector3(0, 0, -distance_from_top),direction = mp.Y, size=mp.Vector3(sx, sy, 0),weight= -1), # Bottom z side
    
    # Sides, x-side (upwards)
    mp.Near2FarRegion(center=mp.Vector3(distance_from_center, 0 , 0),direction = mp.Z, size=mp.Vector3(sx, 0, height_to_cover),weight= 1), # +x side
    mp.Near2FarRegion(center=mp.Vector3(-distance_from_center, 0 , 0),direction = mp.Z, size=mp.Vector3(sx, 0,height_to_cover),weight= -1), # -x side
    
    # Sides, y-side (Flat)
    mp.Near2FarRegion(center=mp.Vector3(0, distance_from_center, 0), direction = mp.Z, size=mp.Vector3(0, sy, height_to_cover),weight= 1), # +y side
    mp.Near2FarRegion(center=mp.Vector3(0, -distance_from_center, 0), direction = mp.Z, size=mp.Vector3(0, sy, height_to_cover),weight= -1)  # -y side
    ]

# The cell coordinates for reference
print(cell_size)
# To print the coordinates of the farfield quickly
print("values of the centers")
for i in range(len(near2far_box)):
    # print('values of the centers:')
    print(near2far_box[i].center)
print("values of the sizes")
for i in range(len(near2far_box)):
    # print('values of the size:')
    print(near2far_box[i].size)

# Initialize the simuion
sim = mp.Simulation(
    cell_size=cell_size,
    geometry=geometry,
    k_point = mp.Vector3(0,0.5/a,0), # Xm-point
    sources=source,
    resolution=resolution,
    boundary_layers=pml_layers)

# Add near2far regions to the simulation
nfreqs = 10
nprds = 10
nf = sim.add_near2far(fcen, df, nfreqs,*near2far_box, nperiods = nprds)

# 3-D Rendering
sim.plot3D()

# Run the simulation
sim_time = 200
sim.run(until = sim_time) # *issue occurs here*

# 'similar to the Antenna-radiation.py tutorial'
far_x = np.linspace(-20, 20, 100)  # far-field observation points
far_y = np.linspace(-20, 20, 100) # we increased this to get the farther far fields
far_z =  10 # observation in the xy-plane
X, Y = np.meshgrid(far_x, far_y)

#_Method I for finding far_field_E_x_
def get_far_fields(
        near2far,
        far_a:np.ndarray=None,
        far_b:np.ndarray=None,
        far_c:np.ndarray= 10,
        component:mp.Vector3Type =None
        ):
    nf = near2far
    if far_a is None or far_b is None or far_c is None:
        raise ValueError("far_a, far_b, and far_c must be provided.")
    
    valid_components = [mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz]
    if component not in valid_components:
        raise ValueError(f"Invalid component. Must be one of {valid_components}.")
    get_far_fields.component = component
    
    # Initialize an array to store the far field components
    far_field = np.zeros((len(far_a), len(far_b)))

    for i,x in enumerate(far_a):
        for j,y in enumerate(far_b):
            r = mp.Vector3(x, y, far_c)
            far_field_components = sim.get_farfield(nf, r)
            far_field[i, j] = far_field_components[component]
           return far_field

far_field_Ex = get_far_fields(nf,far_x,far_y,far_z,mp.Ex)
far_field_Ey = get_far_fields(nf,far_x,far_y,far_z,mp.Ey)
far_field_Ez = get_far_fields(nf,far_x,far_y,far_z,mp.Ez)

# Method II as per the tutorial
npts = 100 # number of points in [0,2*pi) range of angles
angles = 2*np.pi/npts*np.arange(npts)
r = 1000 / fcen
E = np.zeros((npts,3),dtype=np.complex128)

for n in range(npts): #freq ??
    ff = sim.get_farfield(nf, mp.Vector3(r * np.cos(angles[n]), r * np.sin(angles[n])))
    E[n,:] = [ff[j] for j in range(3)]

# PLOTTING
# Extract the real and imaginary parts of the far-field components
# ignore code for plotting since I never get that far 
# for the sake of consistency let
far_field = far_field_Ex 
U = np.real(far_field)
V = np.imag(far_field)

# Normalize U and V for better visualization
U_norm = U / np.max(np.abs(U)) # norm of real components 
V_norm = V / np.max(np.abs(V)) # norm of imag components

far_field_magnitude = np.sqrt(U_norm**2 + V_norm**2)

fig, ax = plt.subplots(figsize= (8,8))

# Create the quiver plot
quiver = ax.quiver(X, Y, U_norm, V_norm,far_field_magnitude, cmap='plasma', scale= 40)

# Add a color bar
cbar = fig.colorbar(quiver, ax=ax)
cbar.set_label('|E|')

# Set labels and title
ax.set_xlabel('x (µm)') # position
ax.set_ylabel('y (µm)') # position 
ax.set_title(f'Far-field Quiver Plot at z={far_z}')

# Save the plot
plt.savefig("far_field_quiver.png")
plt.show()

Vispy or MEEP insists those coordinates are what they should be ....