WIP: IIR Filter Susceptibilites

[smartalecH]

Preliminary attempt that addresses #398. By enabling a susceptibility profile that just takes a ratio of polynomials, we can fit some of the more exotic dispersion profiles and even save on computation in some cases.

My methodology is to let the user specify the coefficients of the zeros (numerator polynomial) and the poles (denominator polynomial) in the S domain. Then, when the susceptibility is initialized, it transforms the coefficients to the Z domain (currently with a bilinear transform).

I feel it's important that MEEP performs the transformation (rather than the user a priori) since it's dependent on the time step size (which is dependent on the resolution and Courant factor).

Once the coefficients are calculated, a custom structure will hold the right number of feedforward (W) and feedback (P) coefficients. The rest of the routine can then follow the Lorentzian routine closely.

Once this is implemented, it can be easy to create other susceptibilities (i.e. conjugate pole, Debye, etc.) that map to this.

Any thoughts?

[stevengj]

See #970 on how to rebase.

[smartalecH]

I've gone through the code pretty thoroughly, but I'm a little stuck. The new IIR implementation seems to be unstable for a simple Lorentzian (this is my first benchmark).

I'm using the IIR update equation used by both Matlab and SciPy (Transposed Direct II). I convert the S domain to the Z domain using a bilinear transform (an algorithm used by SciPy). I print out the Z domain coefficients and test them in Matlab. They are stable. I'm even able to simulate a step response using the same filter algorithm. I'm assuming that there's additional feedback I'm not aware of somewhere that is making the filter explode within meep.

The meat of the algorithm is within update_P() inside of iir_susceptibility.cpp.

Here's a simple test script:

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


# Default lorentzian for glass (pulled from materials library)
min_wavelength = 0.25
max_wavelength = 1.77
um_scale = 1.0
SiO2_range = mp.FreqRange(min=um_scale/max_wavelength, max=um_scale/min_wavelength)
SiO2_frq1 = 1/(0.103320160833333*um_scale)
SiO2_gam1 = 1/(12.3984193000000*um_scale)
SiO2_sig1 = 1.12
SiO2_susc = [mp.LorentzianSusceptibility(frequency=SiO2_frq1, gamma=SiO2_gam1, sigma=SiO2_sig1)]
SiO2 = mp.Medium(epsilon=1.0, E_susceptibilities=SiO2_susc, valid_freq_range=SiO2_range)

# Equivalent sus for glass using IIR filter
num = [(2*np.pi*SiO2_frq1)**2]
den = [1,2*np.pi*SiO2_gam1,(2*np.pi*SiO2_frq1)**2]
SiO2_iir_susc = [mp.IIR_Susceptibility(num, den, sigma=SiO2_sig1)]
SiO2_iir = mp.Medium(epsilon=1.0, E_susceptibilities=SiO2_iir_susc, valid_freq_range=SiO2_range)

# -------------------------------------------------- #
# Calibration run
# -------------------------------------------------- #
print("Simulating a calibration run:")
cell_size = mp.Vector3(z=16)
resolution = 50
fcen = 0.5 * (1/max_wavelength + 1/min_wavelength)
df = 0.2*fcen
src = [mp.Source(
    src=mp.GaussianSource(fcen,fwidth=df),
    component=mp.Ex,
    center = mp.Vector3(z=-4)
)]
pml_layers = [mp.Absorber(thickness=2.0)]
sim = mp.Simulation(
    resolution = resolution,
    boundary_layers=pml_layers,
    cell_size=cell_size,
    dimensions = 1,
    sources=src
)

FR = mp.FluxRegion(center=mp.Vector3(z=4),size=mp.Vector3())
transmission_cal = sim.add_flux(fcen,df,100,FR)

sim.run(until_after_sources=400)

freqs = np.squeeze(mp.get_flux_freqs(transmission_cal))
flux_cal_data = sim.get_flux_data(transmission_cal)
flux_cal = np.squeeze(mp.get_fluxes(transmission_cal))

# -------------------------------------------------- #
# Lorentzian simulation
# -------------------------------------------------- #
print("Simulating the Lorentzian:")
def med_func(pt):
    if pt.z > 0:
        return SiO2
    else:
        return mp.Medium(epsilon=1)

geom = [mp.Block(center=mp.Vector3(z=0),size=mp.Vector3(z=mp.inf), material=med_func)]

sim.reset_meep()
sim = mp.Simulation(
    resolution = resolution,
    boundary_layers=pml_layers,
    geometry = geom,
    cell_size=cell_size,
    dimensions = 1,
    sources=src,
    extra_materials = [SiO2]
)

transmission_lor = sim.add_flux(fcen,df,100,FR)
sim.run(until_after_sources=400)
flux_lor = np.squeeze(mp.get_fluxes(transmission_lor))
T_lor = flux_lor / flux_cal

# -------------------------------------------------- #
# IIR simulation
# -------------------------------------------------- #
print("Simulating the IIR filter:")
def med_func(pt):
    if pt.z > 0:
        return SiO2_iir
    else:
        return mp.Medium(epsilon=1)

geom = [mp.Block(center=mp.Vector3(z=0),size=mp.Vector3(z=mp.inf), material=med_func)]

sim.reset_meep()
sim = mp.Simulation(
    resolution = resolution,
    boundary_layers=pml_layers,
    geometry = geom,
    cell_size=cell_size,
    dimensions = 1,
    sources=src,
    extra_materials = [SiO2_iir]
)

transmission_iir = sim.add_flux(fcen,df,100,FR)
sim.run(until_after_sources=400)
flux_iir = np.squeeze(mp.get_fluxes(transmission_iir))
T_iir = flux_iir / flux_cal


# -------------------------------------------------- #
# Compare results
# -------------------------------------------------- #

plt.figure()
plt.plot(1/freqs,T_lor,label='Lorentzian')
plt.plot(1/freqs,T_iir,label='IIR')
plt.legend()
plt.show()

[oskooi]

For consistency, note the Meep convention of specifying the ordinary frequency rather than the angular frequency in the interface resulting in:

num = [(SiO2_frq1)**2]
den = [1,SiO2_gam1,(SiO2_frq1)**2]

With this change and a resolution of 1000, the results are at least similar in magnitude:

[stevengj]

try replacing w_cur with white noise to see if this filter implementation still blows up (in particular, you can check both whether the resulting colored-noise polarization blows up, and whether the resulting electromagnetic fields blow up)