fancompute / legume Goto Github PK

So far only real/imaginary frequencies work.

In the example of 'Bands of a photonic crystal slab', there is a command 'phc.layers[-1].add_shape(legume.Circle(eps=1.0, r=r))',
Could you please explain the meaning of 'layers[-1]' in the command? Thank you very much!

Currently, the grid of reciprocal lattice vectors used in the PWE/GME expansion is always implemented as all possible combinations of some linspace of b1 vectors and another linspace of b2 vectors, where b1 and b2 are the elementary reciprocal lattice vectors. This makes it easy to work with general lattices, and also ensures that the permittivity matrix is always Toeplitz-block-Toeplitz. However, for the specific case of a hexagonal lattice, this grid is not hexagonal (it looks instead like a parallelogram), and so for example degenerate modes at high-symmetry points will not be exactly degenerate (they will be only in the limit of infinite gmax). It could be useful to have an option to use a hexagonal subspace of reciprocal lattice vectors in order to preserve the symmetry of the lattice.

Would be good to check the computed fields of the eigenmodes vs. some first-principle simulation.

Note that the H-field should be the most "correct" as all components are continuous across all interfaces. The D and especially the E fields, because of the way the expansion works, will have some spurious discontinuities even in cases when they're not supposed to. For example if a patterned and un-patterned slab of the same material sit on top of each other, there is no "interface" in the region where they have the same permittivity, but because the effective permittivity of the patterned slab is lower, there's a discontinuity in the guided modes, and so there will be one in the PhC modes too.

How to debug such error?

TypeError: loop of ufunc does not support argument 0 of type ArrayBox which has no callable exp method

So far I encountered when starting some optimization parameters at zero (unchanged system) and sometimes when reducing the tunable parametric space.

The guided-mode dispersion typically has the strongest variation close to zero wave-vector. It's also hard to resolve the modes very close to zero. We could add an adaptive grid that has more points in that region.

legume.viz.field can be used to plot the field pattern. Sometimes, we need to export the data to do further computation, such as integration. So, how to export the field data of some modes to a txt file or other file types?

In many of the optimization problems we will consider, the change of a layer's average permittivity w.r.t. the in-plane degrees of freedom being optimized could be significant. We should think about how to wrap the call to scipy.optimize.brentq() in slab_modes.py.

Currently, I think this is limiting the performance of the grating optimization I'm testing in opt_grating_polygons.py.

We can use the python logging package for verbose/non-verbose/error/warning reporting: https://docs.python.org/3.7/library/logging.html

This method is available in Nov 2019, but is gone in Feb commits. Didn't find any documentation on it.
What's the updated way to plot this?

Currently, the permittivity FT matrix inversion for the PWE and the GME is done using numpy.linalg.inv. However, the matrix has the structure of a Toeplitz-Block-Toeplitz matrix, which can be used for a faster inversion as per this paper . I also have a MATLAB code for that. It could be implemented at some point (ideally working with the autograd backend too) as this operation could become costly for large elementary cells.

I'm testing the 2D plane wave expansion of legume computing the TM modes of a simple 2D slab waveguide, a problem that has an analytical solution.

This is my code:

import numpy as np
import matplotlib.pyplot as plt

from legume import PlaneWaveExp, Circle, ShapesLayer, Lattice, viz
import legume


eps_waveguide = 3.5**2  
eps_background = 1

d = 1.0  # waveguide width
w = 0.5
# Initialize lattice
lattice = Lattice(np.array([w,0.]) ,np.array( [0.,8.0])  )

# Initialize layer
layer = ShapesLayer(lattice, eps_b = eps_background)
layer.add_shape(legume.Poly(x_edges=[-0.5*w, 0.5*w, 0.5*w, -0.5*w], y_edges= [-0.5*d, -0.5*d, 0.5*d, 0.5*d], eps=eps_waveguide))

kx = np.linspace(-2.0, 2.0, 2)
ky = kx*0
kpoints = np.c_[kx,ky]
pwe = PlaneWaveExp(layer, gmax=7)
viz.eps_ft(pwe, figsize = (3,8))
path = lattice.bz_path(kpoints, [1])

pwe.run(kpoints=path['kpoints'], pol='tm', numeig=2)
freqs_tm = pwe.freqs

print('TM freqs: ', freqs_tm[-1,:])

The code outputs the correct TM frequencies [0.11859563, 0.19938521], which match the ones computed with the analytical solution.

However, when I visualize the Ez and Dz fields, I found the Ez is discontinuous and Dz is continuous at the waveguide boundaries.

viz.field(pwe, field='D', kind=-1, mind=0,
            component='z', val='abs', N1=100, N2=100, cbar=True, eps=True);
viz.field(pwe, field='E', kind=-1, mind=0,
            component='z', val='abs', N1=100, N2=100, cbar=True, eps=True);

Since these components are perpendicular to the XY plane where the waveguide is defined, from Maxwell equations, it's expected Ez is continuous since the electric field must fulfill at the waveguide boundaries $\hat{n} \times\left(\overline{\mathrm{E}}_1-\overline{\mathrm{E}}_2\right)=0$. From this relation, it can also be seen that Dz must be discontinuous.

Why am I obtaining the opposite results? Is a problem related with legume or am I misinterpreting some output?

In the section ”How do I incorporate symmetry?“, the author gives the mode index of the quai-TE/TM modes of photonic crystal slabs. Which plane is the reference plane for this classification?

We should have some mechanism for exporting the results of a topology optimization (defined across a polygon grid) into a standard file format.

The motivation for this would be to make it easy to import the structures into other codes. For example, COMSOL can accept STL files and GDS files. Also, way way way down the road it may be desirable to have GDS output for fabrication purposes.

Hello,

I'm trying to add two fourier shapes to the same layer one next to the other with different x_cent and y_cent parameters however, the shapes overlap. Please check the following example:

lattice = legume.Lattice('hexagonal')
def shamrock_phc(d, r0, rd):
    phc = legume.PhotCryst(lattice)
    phc.add_layer(d=d, eps_b=12)
    f_as = np.array([2*r0,   0,  0, rd])
    shamrock = legume.FourierShape(x_cent=0, y_cent=0, f_as=f_as)
    phc.add_shape(shamrock)
    f_as = np.array([2*r0,   0.01,  0, rd])
    shamrock2 = legume.FourierShape(x_cent=0.5, y_cent=0.5, f_as=f_as)
    phc.add_shape(shamrock2)
    return phc
phc = shamrock_phc(0.5, 0.3, 0.05)    
gme = legume.GuidedModeExp(phc, gmax=5)
legume.viz.structure(phc, figsize=2., cbar=False, Nx=200, Ny=300)
legume.viz.eps_ft(gme, figsize=2., cbar=False, Nx=200, Ny=300)

Thanks,
Shualmit

Would be good to have a test verifying the guided mode calculation, especially in the case of multiple coupled slabs.

Currently the guided-mode computation assumes non-degenerate modes. The "resolution" in looking for those is given by gmode_step. The computation becomes really heavy when pushing gmode_step down to values below e.g. 1e-5. In some structures, like two identical (effective) slabs separated by an air gap, the modes can be almost degenerate, and very hard to resolve. There's no obvious solution to that but putting this here as something to think about in the future.

when I use the following command to plot the field profile, I find I can only successfully plot the field when N1=N2, otherwise it will return an error.

"legume.viz.field(gme, 'e', 0, 0, z=D1/2, component='x', val='abs', N1=50, N2=100, cbar=True, eps=True)"

The detailed information about the error is "ValueError: operands could not be broadcast together with shapes (50,) (200,) "

How to solve this problem?

The documentation (FAQ) contains a minor error, saying that : gmode_eps = 'custom' is used for custom eps. It seems to actually be: eps_eff= 'custom'.

Thank you for the great tool.

Hello,
I am trying to implement an optimization whose FOM mixes three differentiaded sub-FOMs, among which the Q-factor; he >problem I am facing is an inconsistency between the autograd and the finite-difference gradient for this quantity:

def of_daisy(params):
DbottomINT = params[0]
DmiddleINT = params[1]
DtopINT = params[2]
DgratingINT = params[3]
DspacerINT = params[4]
fillAINT = params[5]
fillBINT = params[6]
phc = struttura(DbottomINT, DmiddleINT,DtopINT, DgratingINT,DspacerINT,fillAINT,fillBINT)
#path = phc.lattice.bz_path(['G', 0.01*np.array([np.pi, 0])], [100])
neig = 20# number of Bloch bands to store
gmax = 10 # truncation of reciprocal lattice vectors
gme = legume.GuidedModeExp(phc, gmax=gmax)
options = {'gmode_inds': [1], 
       'numeig': neig,
       'verbose': False,
       'gradients': 'exact',
       'compute_im': False
        }
gme.run(kpoints=np.array([[0], [0]]), **options)
(freq_im, _, _) = gme.compute_rad(0, [1])
Qinv = 2*freq_im[0]/gme.freqs[0, 1]
 ....
 ....
 ....
 return npa.divide(Qinv,FOM2)

The starting point autograd gradient is:
Autograd gradient : [-4.79306617e-17 -6.15112936e-18 -4.55829281e-17 -1.78071921e-17
-5.44515292e-17 -7.66094519e-18 1.08829952e-17]

The finite-difference gradient is:
Numerical gradient: [-3.30385247e-16 -1.63853199e-16 -2.85469258e-18 -6.53521882e-16
-1.30974041e-16 -2.60724278e-16 1.13700734e-16]

Could you help somehow?

I suspect it is either because of multiple layers, or because cladding is not 1. An example:

import numpy as np
import legume
from legume.phc import PhotCryst, Lattice
from legume.shapes import Square, Circle, Poly
from legume.gme import GuidedModeExp
from legume.utils import plot_reciprocal
import matplotlib.pyplot as plt


# Initialize a lattice (can be 'square', 'hexagonal', or defined by primitive vectors)
lattice = Lattice('square')
# Initialize a PhC (optional permittivity in the lower and upper cladding)
phc = PhotCryst(lattice, eps_l=3.2826**2, eps_u=3.2826**2)

phc.add_layer(d=1.007, eps_b=3.5032**2)  # PhC layer

# A triangle
triangle = Poly(x_edges=[0., 0.5256, 0.], y_edges=[0., 0., 0.5256], eps=1.)
phc.layers[-1].add_shape(triangle)

phc.add_layer(d=0.375, eps_b=3.5032**2)  # sch layer
phc.add_layer(d=0.495, eps_b=3.3955**2)  # active layer

# We can plot an overview of what we've built so far
phc.plot_overview()

# Initialize GME and plot what we're simulating
gmax = 5  # plane wave truncation
gme = legume.GuidedModeExp(phc, gmax=gmax)
gme.plot_overview_ft()

n_k = 21
path = phc.lattice.bz_path(['G', np.array([np.pi/2., 0])], [n_k-1])

# See the legume README to understand the options below
neig = 12
options = {'gmode_inds': np.arange(8), 'gmode_npts': 4000, 'numeig': neig, 'verbose': True}


gme.run(kpoints=path.kpoints, **options)

ax = legume.viz.bands(gme)

ks = np.arange(n_k)
freqs_i = []
for kind in ks:
    # kind = 0  # index of the k point
    minds = np.arange(0, neig)  # mode indexes whose loss rates will be computed
    (freqs_im, coup_l, coup_u) = gme.compute_rad(kind=kind, minds=minds)
    freqs_i.append(freqs_im)
    # Note that states below light line still have zero losses
freqs_i = np.array(freqs_i)
qf = gme.freqs/2./freqs_i
# print("Frequencies of first 4 modes at kx = %f are: " % gme.kpoints[0, kind], gme.freqs[kind, minds])
# print("Loss rates of first 4 modes at kx = %f are : " % gme.kpoints[0, kind], freqs_im)

f1 = plt.figure()
ax = f1.add_subplot(111)
ax.plot(np.log10(qf))

Expect: Q factor at gamma to be in the 10^5 range.

Actual outcome: Q factor in 10^2 range.

Is it expected that this mode not be returned? I'm sort of surprised given that we currently solve for all eigenvalues.

Ultimately, if we switch to a "targeted" approach to eigensolving this will be less of an issue. With the way the frequencies are sorted right now though, this missing mode makes defining objective functions around gamma weird because the band indices always shift by one w.r.t. the first k-point.

Hello, I am having problems reproducing the results in Molding the Flow of Light using this code, I don't know how to set the lattice for this structure. According to the current parameters, the results are quite different from the book. How can I reproduce the results correctly?
Here I have uploaded the simulation code（simulation code.zip)

(Chapter07, page128)

（simulation result）

Regarding taking the curl of the H field in real space to get E/D fields: #24 @momchilmm

I'm currently trying to implement the method of taking the curl of the H field in real space to first get the dielectric field and then the electric field - but I'm getting some odd results and wanted to make sure I'm not on a fools errand with the method I'm using.

To get the D field in xz plane at y1, I call get_eps_xz to get the H field in three xz planes at y1, y1+step, y1-step. Then with that data I calculate the curl at all the points in the y1 plane, with the derivatives in the y direction taken across these planes. Is this what you had in mind in #24?

Thanks

The backprop through the guided mode computation is currently a bottleneck in the gradient computation. There are some things that could be tried to speed things up, including vectorizing the computation w.r.t. g.

Hi,

I have been tinkering around with the Legume software to optimize the design of an L3 photonic crystal cavity to improve its quality factor (Q). In order to validate the design, I am comparing the quality factor of the fundamental cavity mode obtained via Legume and an FDTD solver (Tidy3D). With the initial design, Legume gave me Q ~111k using only 1 k-point and ~170k when averaging over a 2x2 k-point grid (gmax = 2.5). Whereas I got a Q ~ 150k using Tidy3D.

Running the optimization routine, the quality factor increased to Q ~ 1 million (averaging over 2x2 k-point grid, gmax = 2.5), whereas the FDTD solver gave me a Q ~ 255k, which is a huge difference.

I understand that Legume is an approximate solver and that there are several parameters at play here. But is there an intuition on why there could be such a huge difference between the obtained value of Q, particularly for the optimized cavity, and how it can be mitigated?

Recently, I read a PRB paper with the title "Bound states in the continuum in symmetric and asymmetric photonic crystal slabs". I try to repeat the result in Fig. 2(b) and (c). I found the results given by legume do not include the second BIC mode. The paper and my program are attached. Is there any problem in my program?
PhysRevB.101.155303.pdf
BIC in 1D grating.zip

Ideally, all plot_ methods should go in viz.py. The phc, gme, and pwe methods could instead have get_eps methods that the plotting functions use. This would allow for the user to also have access to the raw data and use their own plotting methods.

Hello,
I have checked the examples, but I am wondering if it is possible to calculate the bands for a structure with a bottom mirror. For example a slab with a pattern with air above and below, but as boundary condition or adding an additional layer with a the bottom a mirror. Something like this,

Open/Air/Slab/Air/Mirror

Thank you,

Hello,

I'm trying to set a figure and gridspec using the function eps_ft, upon which I get the error "NameError: name 'ars' is not defined" from within the function.

gs = mpl.gridspec.GridSpecFromSubplotSpec(len(ars), 1, gridspec)
NameError: name 'ars' is not defined

Looking at the eps_ft function itself, it seems 'ar' is defined at line 680 but never 'ars'.

Thanks.

Not sure if it's just me, can someone check? Here is what I did:
download source zip. unzip. in console, enter the directory, and type pip install . The package installs fine, but the sub-directiory gme/ is missing in the installed packaage.

Hello,

When I was installing Legume my "pip" attempts failed as it was unable to find the "requirements.txt" file. When I download the .tar.gz version of the project from github this file is also missing from that. Adding this file in manually solved all the problems and it has now installed fine. Perhaps this file is somehow missing from the tar.gz version of the project.

Thanks for the great tool, now its installed looking forward to using it.

Add an option to either plot the periodic parts of the fields only, or the full fields with the exp(i k rho) dependence.

Hi,

I tried to visualize the y component of the E field in the example "Guided mode expansion with Autograd". But the field is not confined in the cavity. What should I do to get the correct field distribution?

Thanks in advance.

When running the current version of opt_grating_polygon.py, I get the following error

~/drive/Code/legume/legume/gme.py in construct_mat(self, k)
    415                                                                         As2, Bs2, chis2, indmode1, oms1, As1, Bs1,
    416 									chis1, pq, -1j)
--> 417                                         mat_block = np.conj(np.transpose(mat_block))
    418
    419                                 mat_blocks[im1].append(mat_block)

AttributeError: 'ArrayBox' object has no attribute 'conjugate'

The variables coup_l and coup_u could be useful, but a bit more work needs to be done to return something meaningful.

Hi, I've been using Legume for a few days and toying with the exemples a bit. But I'm facing an issue with the GuidedModeExp class.

In the band of a photonic crystal slab exemple, when I change the line :
path = lattice.bz_path(['G', 'M', 'K', 'G'], [15, 10, 15])

to :
path = lattice.bz_path(['G', 'K', 'M', 'G'], [15, 10, 15])

I obtain the following error :

File` "c:\python\python37\lib\site-packages\legume\gme\gme.py", line 909, in compute_rad
count+self.modes_numg[kind][im1]][:, np.newaxis])

ValueError: operands could not be broadcast together with shapes (121,0) (120,1)

Is it normal that I can't change the path in the k-space ?

Paul

Just a note for the future when this package is more complete.

I think the standard python package directory structure is to have all package code (currently just core/) under a folder named after the package itself. Here, I guess this would amount to having a folder called pygme/ at the base of the repo. Then, core/ would become pygme/core/. You would also want a pygme/__init__.py.

Then, in the notebooks within this repo, you would have:

import sys
sys.path.append('../')
from pygme.core.pwe import PlaneWaveExp

This way, as long as the repo folder is in PYTHONPATH or if it gets installed as a pip package, you can also access the package externally through the pygme namespace.

It might be possible to significantly improve the memory usage especially when computing gradients, but how exactly requires some thought.

I am not 100% sure if it's related to multi-layer. I've tried single layer it seems to work. Here is a multi-layer example where it raises error:

import numpy as np
import legume
from legume.phc import PhotCryst, Lattice
from legume.shapes import Square, Circle, Poly
from legume.gme import GuidedModeExp
from legume.utils import plot_reciprocal
import matplotlib.pyplot as plt


# Initialize a lattice (can be 'square', 'hexagonal', or defined by primitive vectors)
lattice = Lattice('square')
# Initialize a PhC (optional permittivity in the lower and upper cladding)
phc = PhotCryst(lattice, eps_l=3.2826**2, eps_u=3.2826**2)

# phc.add_layer(d=3.208, eps_b=3.2826**2)  # upper cladding
phc.add_layer(d=1.007, eps_b=3.5032**2)  # PhC layer

# A triangle
triangle = Poly(x_edges=[0., 0.5256, 0.], y_edges=[0., 0., 0.5256], eps=1.)
phc.layers[-1].add_shape(triangle)

phc.add_layer(d=0.375, eps_b=3.5032**2)  # sch layer
phc.add_layer(d=0.495, eps_b=3.3955**2)  # active layer
# phc.add_layer(d=4.6416, eps_b=3.2826**2)  # lower cladding

# We can plot an overview of what we've built so far
phc.plot_overview()

# Initialize GME and plot what we're simulating
gmax = 5  # plane wave truncation
gme = legume.GuidedModeExp(phc, gmax=gmax)
gme.plot_overview_ft()

n_k = 21
path = phc.lattice.bz_path(['G', np.array([np.pi/2., 0])], [n_k-1])

# See the legume README to understand the options below
neig = 12
options = {'gmode_inds': np.arange(12), 'gmode_npts': 4000, 'numeig': neig, 'verbose': True}

gme.run(kpoints=path.kpoints, **options)

Raises error:

Traceback (most recent call last):
  File "T:\legume\legume\gme.py", line 320, in run
    mat = self.construct_mat(k=k)
  File "T:\legume\legume\gme.py", line 397, in construct_mat
    (indmode2, oms2, As2, Bs2, chis2) = self._get_guided(gk, mode2)
  File "T:\legume\legume\gme.py", line 138, in _get_guided
    mode//2, self.omegas_te, self.coeffs_te)
  File "T:\legume\legume\gme.py", line 121, in interp_guided
    gs = self.g_array[-len(omegas[im]):]
IndexError: list index out of range

Currently it's up to the user to make sure that shapes added to a ShapesLayer are not overlapping. If they are, the permittivity of both those shapes is added in that region, which is probably not expected behavior.

It might be useful to set the permittivity to the one of the last shape, but it's not clear how exactly since the shapes are defined directly through their Fourier components (no real-space grid). Or at least we could raise a warning if something like this is happening.

Also this could be important for optimizations. If overlap can't be handled properly (see above), it could be useful to somehow implement constraints that allow you to enforce that the shapes are not crossing as the optimization progresses.

MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
linewidth=markeredgewidth)

How to solve it?

Currently, the eigenvalues are computed using numpy.linalg.eigh, which returns all eigenvalues of a matrix. For large matrices, it could be much more efficient to use instead something like scipy.sparse.linalg.eig, where only a specified number of eigenvalues is computed. To make this compatible with the 'autograd' backend, an autograd primitive should be defined for that function, which shouldn't be too hard starting from the primitive for numpy.linalg.eigh, which can already be found in pygme.primitives.py.

Hey @momchil-flex,
I have been trying to figure out a way to differentiate over eigenproblem solvers, that use sparse matrices and only return a set of k eigenvalues. I found your implementation here, that uses Steven Johnsons method.
I have two questions:
This is equation 7 from Steven Johnsons writeup:

If I understand it correctly the latter part is the vjp:
$$v^T x_p = - \lambda_0^T A_p x + g_\alpha x^T A_p x$$
How does $A_p$ get incorporated in your implementation?

It seems as if $A_p$ is left out and $x^Tx$ is replaced with $xx^T$. I most probably don't fully understand the underlying AD mechanics.

You also stated in a JAX issue (google/jax#3112) that you modified the formulation to account for Hermitian matrices. What modification did you make?

Sorry to bother you...
Regards JD