fancompute/ceviche

:shrimp: Electromagnetic Simulation + Automatic Differentiation

ceviche

Electromagnetic Simulation Tools + Automatic Differentiation. Code for paper Forward-Mode Differentiation of Maxwell's Equations.

What is ceviche?

ceviche provides two core electromagnetic simulation tools for solving Maxwell's equations:

• finite-difference frequency-domain (FDFD)

• finite-difference time-domain (FDTD)

Both are written in numpy / scipy and are compatible with the HIPS autograd package, supporting forward-mode and reverse-mode automatic differentiation.

This allows you to write code to solve your E&M problem, and then use automatic differentiation on your results.

As a result, you can do gradient-based optimization, sensitivity analysis, or plug your E&M solver into a machine learning model without having to go through the tedious process of deriving your derivatives by hand.

Examples

There is a comprehensive ceviche tutorial available at this link with several ipython notebook examples:

1. Running FDFD simulations in ceviche.
2. Performing inverse design of a mode converter.
3. Adding fabrication constraints and device parameterizations.
4. Inverse design of a wavelength-division multiplexer and advanced topics.

There are also a few examples in the examples/* directory.

Installation

There are many ways to install ceviche.

The easiest is by

pip install ceviche

But to install from a local copy, one can instead do

git clone https://github.com/twhughes/ceviche.git
pip install -e ceviche
pip install -r ceviche/requirements.txt

from the main directory.

Alternatively, just download it:

git clone https://github.com/twhughes/ceviche.git

and then import the package from within your python script

import sys
sys.path.append('path/to/ceviche')

Package Structure

Ceviche

The ceviche directory contains everything needed.

To get the FDFD and FDTD simulators, import directly from ceviche import fdtd, fdfd_ez, fdfd_hz

To get the differentiation, import from ceviche import jacobian.

constants.py contains some constants EPSILON_0, C_0, ETA_0, Q_E, which are needed throughout the package

utils.py contains a few useful functions for plotting, autogradding, and various other things.

optimizers.py contains optimizer functions for doing inverse design.

viz.py are functions that help with plotting fields and sructures.

modes.py contains a mode sorter (WIP) that can be used to create waveguide mode profiles for the simulation, for example.

Examples

There are many demos in the examples directory, which will give you a good sense of how to use the package.

Tests

Tests are located in tests. To run, cd into tests and

python -m unittest

to run all or

python specific_test.py

to run a specific one. Some of these tests involve visual inspection of the field plots rather than error checking on values.

To run all of the gradient checking functions, run

chmod +x test/test_all_gradients.sh
tests/test_all_gradients.sh

Credits

If you use this for your research or work, please cite

@article{hughes2019forward,
  title={Forward-Mode Differentiation of Maxwell’s Equations},
  author={Hughes, Tyler W and Williamson, Ian AD and Minkov, Momchil and Fan, Shanhui},
  journal={ACS Photonics},
  volume={6},
  number={11},
  pages={3010--3016},
  year={2019},
  publisher={ACS Publications}
}

Our logo was created by @nagilmer

ceviche's Issues

Total Field - Scattered Field

Is your feature request related to a problem? Please describe.
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]

Describe the solution you'd like
A clear and concise description of what you want to happen.

Describe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.

Additional context
Add any other context or screenshots about the feature request here.

Question about circular polarization

Hello !

Thank you for your amazing work!

I was wondering if it is possible to have a circularly polarized source? Also, is the 3D implementation still in your plans? I see in the code that it is in progress.

Have a nice day.

Questions on spatial (grid) convergence

When I was looking into a simple scenario involving planewave excited by a line source, I noticed that Ceviche somehow requires higher spatial resolution to generate a visibly good planewave. I consider the simplest scenario -- the whole space is vacuum with PML on the y-direction, and periodic in the x-direction. To obtain a good planewave, the pixel size needs to be <lambda/50 for Ceviche, while ~lambda/10 is already good enough for another FDFD software based on Yee-grid without epsilon-averaging. (lambda being the wavelength).

The code is as follows:
`import numpy as np

import matplotlib as mpl
mpl.rcParams['figure.dpi']=100
import matplotlib.pylab as plt

import ceviche
from ceviche.constants import EPSILON_0, C_0, ETA_0

Source angular frequency in Hz

lam = 1e-6
omega = 2np.piC_0/lam

Resolution in nm

dl = lam/20

Simulation domain size in number of pixels

Nx = int(10.*lam/dl)
Ny = Nx

Initialize relative permittivity of the domain

epsr = np.ones((Nx, Ny))

Initialize the source position and amplitude

src_x = np.arange(0,Nx-1)
src_y = int(7.5*lam/dl) * np.ones(src_x.shape, dtype=int)
source = np.zeros((Nx, Ny), dtype=np.complex)
source[src_x, src_y] = 1

Create the simulation object for 'Ez' (TE) polarization

simulation = ceviche.fdfd_ez(omega, dl, epsr, [0,20],bloch_phases=[0.,0.])

Run the simulation with the given source

Hx, Hy, Ez = simulation.solve(source)

Visualize the real value of the Ez field component

ax = ceviche.viz.real(Ez)

Plot the source position

plt.show()`

For dl=lam/20 as above, the planewave will appear as the attached:

Is there a particular reason why Ceviche converges slower? Thanks.

Eigenmode sources for Hz polarization

Currently get_modes() and insert_modes() only support the Ez polarization. The calculation should be extended with a user toggle to select either the Hz or the Ez polarization.

Question on the Mode Overlap function

Hello,

I am a student researcher trying to use your software to design dielectric waveguides. I was hoping someone could please share some insight on the mode overlap function, as I think I may be using it incorrectly. I understand the mode overlap to be the coupling efficiency.

Specifically, I had a concern when simulating a straight waveguide. It is clear that the mode propagates across the waveguide, but the dB10 value of the mode overlap is about -50dB. This value practically stays the same, even as I move the probe (output port) closer to the source. I am using the functions defined in the example notebooks, (e.g., 02_Invdes_intro.ipynb).

The design:

The efficiency as I move the probe closer to the source:

The mode overlap function:

def mode_overlap(E1, E2):
    """Defines an overlap integral between the simulated field and desired field
    """
    return npa.abs(npa.sum(npa.conj(E1)*E2))

What is wrong?

Thanks,
Ricardo Sendrea

Robust PEC / PMC

Right now set with eps_r -> infty. Allow option to zero out matrix elements directly based on where PEC is.

Complex permittivity support?

Hello, I'm not sure if the simulator allows complex permittivity definitions. This feature of material allows to include non linear material as are observers.

If you consider a simple definition in yours examples such as "eps_r=12" to "eps_r=12+2j". Have you got any example that contemplate this feature in other way?? It's because doing this change produces an error

Thank you in advance,

Publish to PyPi

Is your feature request related to a problem? Please describe.
The newest bugfixes/compatibility upgrades present in this repository are not present in the current pypi release (v0.1.2) published August 2020.

Describe the solution you'd like
Could jou please release the current version of ceviche to pypi

Describe alternatives you've considered
Installing directly from github works in most cases with some exceptions. For example publishing a package on pypi with ceviche as dependency.

Thanks for maintaining the project still!!
Regards
JD

3D FDFD

Allow for setting of MuR

Would be useful to set permeability

Contribution with large improvement

Hello,

As a great improvement it would be possible new non-linear materials incorporating your previous simulator "angler", really nice contribution as well!!!

Thank you in advance,

References:
ceviche
Forward-Mode Differentiation of Maxwell’s Equations
Tyler W. Hughes, Ian A. D. Williamson, Momchil Minkov, and Shanhui Fan
Department of Electrical Engineering and Ginzton Laboratory, Stanford University, Stanford, California 94305, United States
ACS Photonics 2019, 6, 11, 3010–3016
Publication Date: October 21, 2019

angler
Adjoint Method and Inverse Design for Nonlinear Nanophotonic Devices
Tyler W. Hughes, Momchil Minkov,‡ Ian A. D. Williamson, and Shanhui Fan
Department of Electrical Engineering and Ginzton Laboratory, Stanford University, Stanford, California 94305, United States
ACS Photonics 2018, 5, 4781−4787
Publication Date: December 3, 2018

---

ceviche.viz.abs Plotting error

Describe the bug
While practicing the inverse design on my workstation where i had pip installed ceviche packages in anaconda3 base environment, I am stuck at a visualization error against the following command referred to first program of workshop.

ceviche.viz.abs(eps_r,cbar=True)

---

ValueError Traceback (most recent call last)
in
----> 1 ceviche.viz.abs(eps_r,cbar=True)

E:\anaconda3\lib\site-packages\ceviche\viz.py in abs(val, outline, ax, cbar, cmap, outline_alpha, outline_val)
32
33 vmax = np.abs(val).max()
---> 34 h = ax.imshow(np.abs(val.T), cmap=cmap, origin='lower left', vmin=0, vmax=vmax)
35
36 if outline_val is None and outline is not None: outline_val = 0.5*(outline.min()+outline.max())

E:\anaconda3\lib\site-packages\matplotlib_init_.py in inner(ax, data, *args, **kwargs)
1436 def inner(ax, *args, data=None, **kwargs):
1437 if data is None:
-> 1438 return func(ax, *map(sanitize_sequence, args), **kwargs)
1439
1440 bound = new_sig.bind(ax, *args, **kwargs)

E:\anaconda3\lib\site-packages\matplotlib\axes_axes.py in imshow(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, filternorm, filterrad, resample, url, **kwargs)
5517 aspect = rcParams['image.aspect']
5518 self.set_aspect(aspect)
-> 5519 im = mimage.AxesImage(self, cmap, norm, interpolation, origin, extent,
5520 filternorm=filternorm, filterrad=filterrad,
5521 resample=resample, **kwargs)

E:\anaconda3\lib\site-packages\matplotlib\image.py in init(self, ax, cmap, norm, interpolation, origin, extent, filternorm, filterrad, resample, **kwargs)
898 self._extent = extent
899
--> 900 super().init(
901 ax,
902 cmap=cmap,

E:\anaconda3\lib\site-packages\matplotlib\image.py in init(self, ax, cmap, norm, interpolation, origin, filternorm, filterrad, resample, **kwargs)
244 if origin is None:
245 origin = mpl.rcParams['image.origin']
--> 246 cbook._check_in_list(["upper", "lower"], origin=origin)
247 self.origin = origin
248 self.set_filternorm(filternorm)

E:\anaconda3\lib\site-packages\matplotlib\cbook_init_.py in _check_in_list(_values, **kwargs)
2264 for k, v in kwargs.items():
2265 if v not in values:
-> 2266 raise ValueError(
2267 "{!r} is not a valid value for {}; supported values are {}"
2268 .format(v, k, ', '.join(map(repr, values))))

ValueError: 'lower left' is not a valid value for origin; supported values are 'upper', 'lower'

To Reproduce
Steps to reproduce the behavior:

1. Go to '...'
2. Click on '....'
3. Scroll down to '....'
4. See error

Expected behavior
A clear and concise description of what you expected to happen.

Desktop (please complete the following information):

• OS: [Windows10]
• Browser [chrome]
• Version [beta]

Question about eigenmodes of a dielectric waveguide

Hi, I would like to use Ceviche to design a dielectric waveguide.
Its inner layer is dielectric, whereas its outer layer is a conductor, shown in the figure below.
I set it to 10 for illustration purpose.

What is the value permittivity that I should set for a conductor ?
In addition, how should I visualize the mode profile after solving it with ceviche.modes.get_modes ?
Thank you very much.

import numpy as np
import matplotlib.pyplot as plt
import ceviche

freq = 300E9
wavelength = 299792458 / freq
omega = 2*np.pi*freq
dl = wavelength / 100
Nx = Ny = 101
Npml = 60
permitivity = 4.41

# Initialize relative permittivity of the domain
epsr = np.ones((Nx, Ny))

xx = np.linspace(-int(Nx/2), int(Nx/2), Nx) * dl
yy = np.linspace(-int(Ny/2), int(Ny/2), Ny) * dl

radius_inner = 450E-6
radius_outer = 500E-6

for idx in range(Nx):
    for idy in range(Ny):
        rr = np.sqrt(xx[idx]**2 + yy[idy]**2)
        if (rr < radius_outer) and (rr >= radius_inner):
            epsr[idx, idy] = permitivity
        if (rr <= radius_conductor) and (rr > radius_outer):
            epsr[idx, idy] = 10

vals, vecs = ceviche.modes.get_modes(epsr, omega, dl, npml=Npml, m=1)

Allow Complex Permittivity

And include a few materials to choose from

Missing 1/jomega in calculation of Hx and Hy for Ez polarization

The computation of Hx and Hy from Ez are missing a 1/1j/omega prefactor. One consequence of this is that a subsequent Poynting vector calculation is incorrect.

I believe the issue occurs on the following lines:
https://github.com/twhughes/ceviche/blob/master/ceviche/fdfd.py#L248-L256
Though, I'm not sure exactly what is stored in info_dict and if omega is available there.

Example in readme doesn't work

The example included near the top of the readme doesn't work and/or is incomplete.

I think I removed it in a previous commit but it now seems to be back through some sort of merging mistake. I've had several people copy/run it and then ask what was wrong, so I think we should either remove it completely, update it to be a complete example, or just link to one of the examples in the sub-folder.

Grid Averaging

Describe the bug
Need to implement grid averaging of Epsilon_r for Hz polarization.