Enable Eigenmode Features with Dispersive Materials

[smartalecH]

The finished version of #475.

Modifies the get_eigenmode routines such that dispersive materials are now supported. Specifically, I modified the get_chi1inv() function, which now evaluates the (real) inverse permittivity by summing all of the susceptibilities at a given point in a given direction for a specific field component at a specific frequency. It should be backward-compatible, as the default frequency is 0 (i.e. only the DC frequency is evaluated in this case).

As mentioned, it only evaluates the real permittivity, since MPB's iterative solver can (currently) only solve for real eigenvalues.

That being said, this method still enables a lot of important features, one of them being S parameter extraction for dispersive materials (provided the materials are mostly dielectric).

Closes #415.

I've implemented a test that compares the actual permittivity, the return value of chi1(), and the "effective index" of a uniform medium (i.e. just the refractive index). All three should be the same.

I check this for Silicon, SIlver, and Lithium Niobate so that we get a range of susceptibilities.

[smartalecH]

It looks like the serial implementation fails because the new routine doesn't like materials that have been rotated. It's able to pull the epsilon at those values just fine (using sim.fields.get_chi1inv()). But when MPB checks the dielectrics, it throws the "invalid dielectric function" error.

[stevengj]

If MPB is giving that error it must be that your ε matrix is not positive-definite.

[smartalecH]

This latest iteration fixes the original bug, which means that the routine can now use the get_eigenmode() routines with dispersive materials.

In addition, I modified the array_slice and the h5fields routines so that the user can optionally specify a frequency and they will evaluate the permittivity/permeability at that frequency. This is essential for the visualization module and useful for other post-processing applications.

The array_slice routines seem to be working fine (I've implemented additional tests within dispersive_eigenmode.py). The h5 routines fail when executed in parallel. For some reason, there's a communication error with the get_chi1inv_at_pt() function in monitor.cpp. This is the same function used for the get_eigenmode routines and the array_slice routines (with no issues there).

I'm a little stumped, so if anyone else with more experience with parallel h5 generation could look at it, that would be great.

[stevengj]

What's going on here? Why not just res = real(1.0 / eps)?

[smartalecH]

What we really want is real(1/n) ^2, which is not the same as real(1.0 / eps) if eps is complex.

[smartalecH]

Output:

Initializing structure...
Working in 2D dimensions.
Computational cell is 16 x 8 x 0 with resolution 10
     block, center = (0,0,0)
          size (1e+20,1,1)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (1,1,1)
-----------
Initializing structure...
Working in 2D dimensions.
Computational cell is 16 x 8 x 0 with resolution 10
     block, center = (0,0,0)
          size (1e+20,1,1)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (1,1,1)
time for set_epsilon = 0.0751069 s
lorentzian susceptibility: frequency=0.000905797, gamma=0
lorentzian susceptibility: frequency=0.88125, gamma=0
time for set_epsilon = 0.0731909 s
lorentzian susceptibility: frequency=0.000905797, gamma=0
lorentzian susceptibility: frequency=3.31657, gamma=0
lorentzian susceptibility: frequency=0.88125, gamma=0
lorentzian susceptibility: frequency=3.31657, gamma=0
-----------
-----------
creating output file "./test_animation-eps-000000.00.h5"...
creating output file "./test_animation-eps-000000.00.h5"...
HDF5-DIAG: Error detected in HDF5 (1.10.4) thread 0:
  #000: H5F.c line 444 in H5Fcreate(): unable to create file
    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1364 in H5F__create(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #002: H5Fint.c line 1615 in H5F_open(): unable to lock the file
    major: File accessibilty
    minor: Unable to open file
  #003: H5FD.c line 1640 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Can't update object
  #004: H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 35, error message = 'Resource temporarily unavailable'
    major: File accessibilty
    minor: Bad file ID accessed
meep: error on line 455 of ../../src/h5file.cpp: error opening HDF5 output file
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[39506,1],0]
  Exit code:    1
--------------------------------------------------------------------------

[stevengj]

Note that this is not really doing the interpolation from the Yee grid correctly, so it is making a first-order error in principle. In practice, that error only shows up at the boundaries between materials, where it is hard to be accurate anyway, so let's just accept the 1st-order error for now and leave it to another PR to do the interpolation more correctly.

[smartalecH]

The PR should be ready for merging. Here's a summary of the features:

• Enable get_eigenmode routines with dispersive, anisotropic materials (excluding gyromagnetic materials, for now).
• get_array and output_h5 routines now accept an omega parameter and can evaluate the permittivity/permeability at that frequency. Defaults to 0, which performs the previous behavior.
• Visualization module can accept an omega parameter to plot the dielectric at that frequency. Defaults to the center frequency of the first source in the sources list.
• Small bug fixes in the visualization module (with updated tests)
• Added a rotate function to the Medium class.

As discussed, a better interpolation scheme is needed to improve the accuracy around the material boundaries.

[stevengj]

doesn't det == 0.0 work?

[smartalecH]

I was just being safe -- do you want it changed?

[stevengj]

det == 0.0 seems clearer and is safe.

(You can't do det == 0, on the other hand, because C++ doesn't let you mix complex<double> with scalars of other types.)

[smartalecH]

Just fixed with 8590db1

[stevengj]

This line looks off to me — it seems like you are doing an element-wise product rather than a matrix multiplication?

Since σ can only be a diagonal matrix currently, a matrix multiplication would multiply by conductivity[cc][cc] (the diagonal element from row cc), no?

[stevengj]

What are the step_func_args that were added to some functions, like output_epsilon? They don't seem to be used.