NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.17k stars 598 forks source link

Point dft monitors report different results from a planar dft monitor #2621

Open FuGuiLee0807 opened 10 months ago

FuGuiLee0807 commented 10 months ago

In my project, I need to calculate the near fields at some specific points, so I added dft point monitors to get the fields. However, when the results were compared with that achieved from a planar dft monitor, a significant difference is appeared.

The verification steps are described below:

  1. Add a planar dft monitor (set _monitortype in the follow code as "plane"), run the simulation, get the dft fields E1 and the coordinates of points where fields are recorded.
  2. Add point dft monitors by two for-loop (set _monitortype in the follow code as "point"),, the center of the point monitors are set as the same to coordinates got from step1. Run the simulation, and get the dft fields E2.
  3. Compare E1 and E2.

There are two main problems:

  1. E2 contains A lot of NAN value.
  2. At some point monitors, the difference between E1 and E2 is unaccepable.

The code for a test case (plane wave propogation in free spapce) ard results are shown below.

# Plane wave propogation in free space. 

import meep as mp
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path

#-----------------Parameter` definition-------------------
PML_thickness = 0.5 # um
wvl = 0.4    # um
dx = 0.02  # grid size, um

# simulation domain size (without PML)
domain_size_x = 0.8
domain_size_y = 1.6
domain_size_z = 0

# simulation domain size (with PML)
cell_size_x = domain_size_x
cell_size_y = domain_size_y + 2*PML_thickness
cell_size_z = domain_size_z

# source size and center point position
source_size_x = cell_size_x
source_size_y = 0
source_size_z = 0
source_center_x = 0
source_center_y = 0.75
source_center_z = 0

# monitor size and center point position [used if monitor_type = "plane"]
monitor_size_x = domain_size_x
monitor_size_y = domain_size_y - 0.1
monitor_size_z = domain_size_z
monitor_center_x = 0
monitor_center_y = 0
monitor_center_z = 0

# point monitor coordinates [used if monitor_type = "point"]
pt_monitor_x = np.loadtxt("x.txt") 
pt_monitor_y = np.loadtxt("y.txt")

# monitor type
monitor_type = "plane" # "point" or "plane"

#-----------------Define simulation domain and resolution--------------
pml_layers = [mp.PML(thickness=PML_thickness, direction=mp.Y)]
cell_size = mp.Vector3(cell_size_x, cell_size_y, cell_size_z) 
resolution = 1/dx  # pixels/um

#-----------------Define soruce-------------
fcen = 1/wvl
source_center = mp.Vector3(source_center_x, source_center_y, source_center_z)
source_size = mp.Vector3(source_size_x, source_size_y, source_size_z)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=0.5*fcen, is_integrated=True),
                     center=source_center,
                     size=source_size,
                     component=mp.Ez)]

#-----------------Parameter configuration-----------------
sim = mp.Simulation(resolution= resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    sources=sources,
                    k_point=mp.Vector3())

if monitor_type == "plane":
    #-----------------Add monitor---------------
    monitor_center = mp.Vector3(monitor_center_x, monitor_center_y, monitor_center_z)
    monitor_size = mp.Vector3(monitor_size_x, monitor_size_y, monitor_size_z)
    field_monitor = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=monitor_center, size=monitor_size)

    #-----------------Run simulation----------------
    sim.run(until_after_sources=100)  

    #-----------------Get fields------------------
    fields = sim.get_dft_array(field_monitor, mp.Ez, 0)
    field_amplitude = np.abs(fields)    
    field_real = np.real(fields)
    field_imag = np.imag(fields)
    field_real_normlized = field_real / field_amplitude  
    field_imag_normlized = field_imag / field_amplitude  
    amplitude_normlized = np.sqrt(np.power(field_real_normlized,2) + np.power(field_imag_normlized,2))

    #-----------------save data------------------
    [x,y,z,w] = sim.get_array_metadata(dft_cell=field_monitor)  
    file_name1 = 'Ez_real.txt'
    file_name2 = 'Ez_imag.txt'
    file_name3 = 'x.txt'
    file_name4 = 'y.txt'
    np.savetxt(file_name1, field_real_normlized)
    np.savetxt(file_name2, field_imag_normlized)
    np.savetxt(file_name3,x)
    np.savetxt(file_name4,y)

elif monitor_type == "point":

    #-----------------Add monitor-----------------
    pt_monitor_x_num = np.size(pt_monitor_x)
    pt_monitor_y_num = np.size(pt_monitor_y)
    read_x = np.zeros([pt_monitor_x_num, 1]) 
    read_y = np.zeros([pt_monitor_y_num, 1]) 
    field_monitor = []
    for index_x in range(pt_monitor_x_num):
        field_monitor.append([])
        for index_y in range(pt_monitor_y_num):
            pt_monitor_center = mp.Vector3(pt_monitor_x[index_x], pt_monitor_y[index_y], 0)
            field_monitor[index_x].append(sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=pt_monitor_center, size=mp.Vector3()))

    #-----------------Run simulation----------------
    sim.run(until_after_sources=100)  

    #-----------------Get fields------------------
    field_real_normlized = np.zeros([pt_monitor_x_num, pt_monitor_y_num])
    field_imag_normlized = np.zeros([pt_monitor_x_num, pt_monitor_y_num]) 
    field_amplitude = np.zeros([pt_monitor_x_num, pt_monitor_y_num]) 
    for index_x in range(pt_monitor_x_num):
        for index_y in range(pt_monitor_y_num): 
            fields = sim.get_dft_array(field_monitor[index_x][index_y], mp.Ez, 0)
            field_amplitude[index_x][index_y] = np.abs(fields)
            field_real = np.real(fields)
            field_imag = np.imag(fields)            
            field_real_normlized[index_x][index_y] = field_real / field_amplitude[index_x][index_y]  
            field_imag_normlized[index_x][index_y]= field_imag / field_amplitude[index_x][index_y]  

    #-----------------save data------------------
    file_name1 ='Ez_real_2.txt'
    file_name2 ='Ez_imag_2.txt'
    np.savetxt(file_name1, field_real_normlized)
    np.savetxt(file_name2, field_imag_normlized)

E1 E2 E3