I am trying to diagnose a NaN error I am getting from my Near2Far function.

[fumblerful]

Github Upload.zip

I am running a script I got from my a colleague but I have modified it a lot. Previous modifications worked fine, but now I am running into an issue where farfield calculations produce NaN. In my diagnosing print commands, it looks like part of my nearfield box is to blame because the first roughly fifth of results are NaN but the rest are fine. I do not know why this is the case because the normal vector for each section of the box looks like it is pointing the right way.

Thank you!

import meep as mp
import matplotlib.pyplot as plt
import numpy as np
from datetime import datetime

gdsII_file = 'dipole.gds'
CELL_LAYER = 1
LEFT_SOURCE_LAYER = 21
RIGHT_SOURCE_LAYER = 15
POLE_RIGHT = 2
POLE_LEFT = 3
PLA_BOTTOM = 4
PLA_MIDDLE_LEFT = 5
PLA_TOP = 6
WIRE_LEFT = 7
WIRE_RIGHT = 8
WIRE_RIGHT_DI_MID_LEFT = 9
WIRE_RIGHT_DI_SANDWICH = 10
WIRE_LEFT_DI_MID_LEFT = 11
WIRE_LEFT_DI_SANDWICH = 12
WIRE_LEFT_CAP = 13
WIRE_RIGHT_CAP = 14
VSWR_WIRE_LEFT = 16
VSWR_BOTTOM_LEFT = 17
VSWR_RIGHT = 18
VSWR_LEFT = 19
VSWR_TOP = 20
VSWR_BOTTOM_RIGHT = 22
VSWR_WIRE_RIGHT = 23
WIRE_RIGHT_DI_MID_RIGHT = 24
WIRE_LEFT_DI_MID_RIGHT = 25
PLA_MIDDLE_RIGHT = 26
PLA_MIDDLE_RIGHT_2 = 27
PLA_MIDDLE_RIGHT_3 = 28
PLA_MIDDLE_LEFT_2 = 29
PLA_MIDDLE_LEFT_3 = 30

pt=mp.Vector3(146,30,0)
dpml = 1
dielectric=mp.Medium(epsilon=2)
pla = mp.Medium (epsilon = 3)
pla_all = 0
t_all = 0
wire= 0
global_time_norm = 20
global_time_reg = mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3)
global_time_radiation = mp.stop_when_fields_decayed(50, mp.Ez, pt, 1e-3) #4*global_time_norm
resolution = 20

def pulse(pulse_w):
    return lambda t: np.exp(-0.5*(t-5*pulse_w)*(t-5*pulse_w)/pulse_w/pulse_w)
def neg_pulse(pulse_w):
    return lambda t: -1.0*np.exp(-0.5*(t-5*pulse_w)*(t-5*pulse_w)/pulse_w/pulse_w)

def main(pulse_width,frequency,delta_f,n_freq,file_title1,file_title2,file_title3,write_out):
    cell = mp.GDSII_vol(gdsII_file,CELL_LAYER,-t_all,t_all)

    l_src_vol = mp.GDSII_vol(gdsII_file,LEFT_SOURCE_LAYER,-wire,wire)
    r_src_vol = mp.GDSII_vol(gdsII_file,RIGHT_SOURCE_LAYER,-wire,wire)

    pole_right = mp.get_GDSII_prisms(mp.metal,gdsII_file,POLE_RIGHT,-t_all,t_all)
    pole_left = mp.get_GDSII_prisms(mp.metal,gdsII_file,POLE_LEFT,-t_all,t_all)
    
    pla_bottom = mp.get_GDSII_prisms(pla,gdsII_file,PLA_BOTTOM,(-t_all-pla_all),-t_all)
    pla_middle_left = mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_LEFT,-t_all,t_all)+mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_LEFT_2,-t_all,t_all)+mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_LEFT_3,-t_all,t_all)
    pla_middle_right = mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_RIGHT,-t_all,t_all)+mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_RIGHT_2,-t_all,t_all)+mp.get_GDSII_prisms(pla,gdsII_file,PLA_MIDDLE_RIGHT_3,-t_all,t_all)
    pla_top = mp.get_GDSII_prisms(pla,gdsII_file,PLA_TOP,t_all,(t_all+pla_all))
    
    wire_left = mp.get_GDSII_prisms(mp.metal,gdsII_file,WIRE_LEFT,-wire,wire)
    wire_right = mp.get_GDSII_prisms(mp.metal,gdsII_file,WIRE_RIGHT,-wire,wire)
    
    wire_right_di_mid_left = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_RIGHT_DI_MID_LEFT,-wire,wire)
    wire_right_di_mid_right = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_RIGHT_DI_MID_RIGHT,-wire,wire)
    wire_left_di_mid_left = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_LEFT_DI_MID_LEFT,-wire,wire)
    wire_left_di_mid_right = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_LEFT_DI_MID_RIGHT,-wire,wire)
    wire_right_di_sandwich_bottom = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_RIGHT_DI_SANDWICH,-3*wire,-wire)
    wire_right_di_sandwich_top = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_RIGHT_DI_SANDWICH,wire,3*wire)
    wire_left_di_sandwich_bottom = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_LEFT_DI_SANDWICH,-3*wire,-wire)
    wire_left_di_sandwich_top = mp.get_GDSII_prisms(dielectric,gdsII_file,WIRE_LEFT_DI_SANDWICH,wire,3*wire)
    
    wire_left_cap = mp.get_GDSII_prisms(mp.metal,gdsII_file,WIRE_LEFT_CAP,3*-wire,3*wire)
    wire_right_cap = mp.get_GDSII_prisms(mp.metal,gdsII_file,WIRE_RIGHT_CAP,3*-wire,3*wire)
    
    
    region_1 = mp.GDSII_vol(gdsII_file,VSWR_WIRE_LEFT,-t_all,t_all)
    #region_2 = mp.GDSII_vol(gdsII_file,VSWR_WIRE_RIGHT,-t_all,t_all)
    near2far_region_bottom_left = mp.GDSII_vol(gdsII_file,VSWR_BOTTOM_LEFT,-t_all,t_all)
    near2far_region_bottom_right = mp.GDSII_vol(gdsII_file,VSWR_BOTTOM_RIGHT,-t_all,t_all)
    near2far_region_right = mp.GDSII_vol(gdsII_file,VSWR_RIGHT,-t_all,t_all)
    near2far_region_left = mp.GDSII_vol(gdsII_file,VSWR_LEFT,-t_all,t_all)
    near2far_region_top = mp.GDSII_vol(gdsII_file,VSWR_TOP,-t_all,t_all)

    
    #Normalization run: no short at the end of the cable
    geometry = wire_left_cap + wire_right_cap + pole_right + pole_left + pla_middle_left + pla_middle_right + wire_left + wire_right + wire_right_di_mid_left+ wire_right_di_mid_right + wire_left_di_mid_left + wire_left_di_mid_right
    sources = []
    sources.append(
        mp.Source(
            mp.GaussianSource(frequency, width=pulse_width, fwidth=delta_f),
            component=mp.Ex,
            volume=l_src_vol,
            amplitude=1
        )
    )
    sources.append(
        mp.Source(
            mp.GaussianSource(frequency, width=pulse_width, fwidth=delta_f),
            component=mp.Ex,
            volume=r_src_vol,
            amplitude=1
        )
    )
    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell.size,
                        boundary_layers=[mp.PML(dpml)],
                        sources=sources,
                        geometry=geometry)
    flux_r1 = mp.FluxRegion(volume=region_1,direction=mp.Y)
    flux_1 = sim.add_flux(frequency,delta_f,n_freq,flux_r1)
    if(write_out):
        sim.run(mp.to_appended(file_title1,mp.at_every(0.5,mp.output_efield_x)),until=global_time_norm)
    else:
        sim.run(until=global_time_norm)
    #Incident fields
    fields_1_norm = sim.get_flux_data(flux_1)
    #Incident fluxes
    flux_1_norm = mp.get_fluxes(flux_1)
    
    sim.reset_meep()
    
    #Data Run: there is a metal short at the end of the cable
    
    geometry = wire_left_cap + wire_right_cap + pole_right + pole_left + pla_middle_left + pla_middle_right + wire_left + wire_right + wire_right_di_mid_left+ wire_right_di_mid_right + wire_left_di_mid_left + wire_left_di_mid_right
    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell.size,
                        boundary_layers=[mp.PML(dpml)],
                        sources=sources,
                        geometry=geometry)
    flux_1 = sim.add_flux(frequency,delta_f,n_freq,flux_r1)
    sim.load_minus_flux_data(flux_1,fields_1_norm)
    
    #Adding radiation pattern:
    nearfield_box = sim.add_near2far(frequency,delta_f,n_freq,
                                     mp.Near2FarRegion(volume=near2far_region_bottom_left,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_bottom_right,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_right,direction=mp.X,weight=+1),
                                     mp.Near2FarRegion(volume=near2far_region_left,direction=mp.X,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_top,direction=mp.Y,weight=+1),
                                    )
    '''
    nearfield_box = sim.add_near2far(frequency,delta_f,n_freq,
                                     mp.Near2FarRegion(volume=near2far_region_bottom_left,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_bottom_right,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_right,direction=mp.X,weight=+1),
                                     mp.Near2FarRegion(volume=near2far_region_left,direction=mp.X,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_top,direction=mp.Y,weight=+1),
                                    )
    '''
    if(write_out):
        sim.run(mp.to_appended(file_title2,mp.at_every(0.5,mp.output_efield_x)),until=global_time_reg)
    else:
        sim.run(until=global_time_reg)
    #Incident fluxes
    flux_1_reg = mp.get_fluxes(flux_1)
    #Frequencies
    flux_freqs = mp.get_flux_freqs(flux_1)
    
    #Show geometry
    '''
    fig = plt.figure(dpi=200)
    sim.plot2D(ax=fig.gca(),eps_parameters={'cmap': 'gray'})
    plt.savefig('geometry.png')
    plt.close
    '''
    R_1 = []
    VSWR_1 = []
    real_freqs = []
    print('--------------' + str(datetime.now()) + '--------------')
    for i in range(n_freq):
        print('VSWR %:' + str(i/n_freq*100))
        if flux_freqs[i]>0 and flux_1_norm[i]!=0:
            lambda_meters = 1.0/flux_freqs[i]/100000.0
            real_f_GHz = 0.3/lambda_meters
            gamma_1 = -flux_1_reg[i]/flux_1_norm[i]
            #print(gamma_1)
            R_1 = np.append(R_1,gamma_1)
            current_1 = (1+gamma_1)/(1-gamma_1)
            VSWR_1 = np.append(VSWR_1,(1+gamma_1)/(1-gamma_1))
            real_freqs = np.append(real_freqs,real_f_GHz)
    
    plt.figure(figsize=(8, 6), dpi=150)
    
    plt.plot(real_freqs,VSWR_1,color='black',linewidth=3.0)
    plt.xlabel("Frequency (MHz)",fontsize=22)
    plt.ylabel("VSWR",fontsize=22)
    plt.xticks(fontsize=22)
    plt.yticks(fontsize=22)
    #plt.axis([15,20,-1,1])
    plt.grid()
    plt.tight_layout()
    plt.savefig('VSWR.png')
    #print(VSWR_1)
    #print(real_freqs)
    
    sim.reset_meep()
    geometry = wire_left_cap + wire_right_cap + pole_right + pole_left + pla_middle_left + pla_middle_right + wire_left + wire_right + wire_right_di_mid_left+ wire_right_di_mid_right + wire_left_di_mid_left + wire_left_di_mid_right
    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell.size,
                        boundary_layers=[mp.PML(dpml)],
                        sources=sources,
                        geometry=geometry)
    flux_1 = sim.add_flux(frequency,delta_f,n_freq,flux_r1)
    sim.load_minus_flux_data(flux_1,fields_1_norm)
    
    #Adding radiation pattern:
    nearfield_box = sim.add_near2far(frequency,delta_f,n_freq,
                                     mp.Near2FarRegion(volume=near2far_region_bottom_left,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_bottom_right,direction=mp.Y,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_left,direction=mp.X,weight=-1),
                                     mp.Near2FarRegion(volume=near2far_region_top,direction=mp.Y,weight=+1),
                                     mp.Near2FarRegion(volume=near2far_region_right,direction=mp.X,weight=+1),
                                    )
    
    if(write_out):
        sim.run(mp.to_appended(file_title3,mp.at_every(0.5,mp.output_efield_x)),until=global_time_radiation)
    else:
        sim.run(until=global_time_radiation)
    
    #Incident fluxes
    flux_1_reg = mp.get_fluxes(flux_1)
    #Frequencies
    flux_freqs = mp.get_flux_freqs(flux_1)
    
    #Show geometry
    '''
    fig = plt.figure(dpi=200)
    sim.plot2D(ax=fig.gca(),eps_parameters={'cmap': 'gray'})
    plt.savefig('geometry.png')
    plt.close
    '''
    R_1 = []
    VSWR_1 = []
    real_freqs = []
    print('--------------' + str(datetime.now()) + '--------------')
    for i in range(n_freq):
        print('Radiation %:' + str(i/n_freq*100))
        if 3>0:#flux_freqs[i]>0 and flux_1_norm[i]!=0:
            lambda_meters = 1.0/(flux_freqs[i]*100000.0)
            real_f_GHz = 0.3/lambda_meters
            gamma_1 = -flux_1_reg[i]/flux_1_norm[i]
            R_1 = np.append(R_1,gamma_1)
            current_1 = (1+gamma_1)/(1-gamma_1)
            VSWR_1 = np.append(VSWR_1,(1+gamma_1)/(1-gamma_1))
            real_freqs = np.append(real_freqs,real_f_GHz)
    
        
    
    
    y_mark = pla_bottom[0].center[1]
    
    npts = 360
    angles = 2*np.pi/npts*np.arange(npts)
    r = 1000
    E = np.zeros((npts,3),dtype=np.complex128)
    H = np.zeros((npts,3),dtype=np.complex128)
    which_freq = 16
    freq_marker_1 = (which_freq-1)*6
    freq_marker_2 = freq_marker_1+3
    print('--------------' + str(datetime.now()) + '--------------')
    for n in range(npts):
        print('npts %:' + str(n/npts*100))
        ff = sim.get_farfield(nearfield_box,mp.Vector3(r*np.cos(angles[n]),r*np.sin(angles[n])+y_mark,0.0))
        print(ff)
        E[n,:] = [np.conj(ff[j]) for j in range(freq_marker_1,freq_marker_2)]
        H[n,:] = [ff[j+3] for j in range(freq_marker_1,freq_marker_2)]
        print(E[n,:])
        print(H[n,:])

    Px = np.real(E[:,1]*H[:,2]-E[:,2]*H[:,1])
    Py = np.real(E[:,2]*H[:,0]-E[:,0]*H[:,2])
    Pz = np.real(E[:,0]*H[:,1]-E[:,1]*H[:,0])
    Pr = np.sqrt(np.square(Px)+np.square(Py)+np.square(Pz))
    directivity = 10.0*np.log10(Pr/max(Pr))
    print(directivity)
    
    plt.figure(2)
    plt.polar(angles,directivity,color='blue')
    
    ax = plt.gca()
    ax.set_rlim(-30,0)
    ax.set_rticks([-30, -15, -3, 0])
    ax.grid(True)
    ax.set_rlabel_position(180)
    ax.tick_params(labelsize=18)
    plt.title(real_f_GHz)
    plt.tight_layout()
    plt.savefig('radiation.png')
    
    
    
main(0.01,0.006,0.03,512,"ex_prop_norm","ex_prop_reg","radiation",0)
#def main(pulse_width,frequency,delta_f,n_freq,file_title1,file_title2,file_title3,write_out):


#The mains below here are from really good runs
'''
This one was the first really good run.  Side spikes are there, but very much less than before.  Center beam is clear and like a balloon.  Still losing a lot in the cavity, though.
main(1.5,0.5,0.1,32,"ex_prop_norm","ex_prop_reg",1)
Second good one.  Increasing Frequncy slightly made beam narrower without affecting quality.
main(1.5,0.75,0.1,32,"ex_prop_norm","ex_prop_reg",1)
'''

[smartalecH]

it looks like part of my nearfield box is to blame because the first roughly fifth of results are NaN but the rest are fine

As in, the near fields themselves are NaN?

[fumblerful]

What I figured was something was causing the near fields to interact oddly with the box I designed and that was causing the NaN. However, I do not know if it is my signal to blame.
When I print my get_farfield of my nearfield box I get a block of nan+nanj and then normal values. I attached a shortened version of that.
[(nan+nanj), (nan+nanj), 0j, 0j, 0j, (nan+nanj), (nan+nanj), (nan+nanj), 0j, 0j, 0j, ... (-7.931952217684334e-08-2.2202470229319037e-08j), 0j, 0j, 0j, (-3.242523036998343e-08+2.714339438704945e-08j)

[smartalecH]

First, check if the near fields themselves are nan. Meep has some checks it enforces every N time steps, but maybe you're between that.

That way you can narrow it down to the near2far routine itself.

[smartalecH]

I'd also check your fields you're using to normalize.