CHANGELOG

Released on Januart 16, 2023
----------------------------
1. Transitioned license from GPLv3 to the much more permissive BSD-3.
2. Updated boost download URL in install.py

Released on May 6, 2019
-----------------------

This is a major update. Please read through each of these changes carefully and
update your code as needed. Code built upon previous version of EMopt are not
guaranteed to run without errors.

1. New versioning format based on YYYY.MM.DD.

2. Changed names of AdjointMethodFM and AdjointMethodPNF to AdjointMethodFM2D
and AdjointMethodPNF2D. This is to distinguish it from the new
AdjointMethodPNF3D implementation. Note: This will break your code, but is
easily fixed by tacking '2D' onto AdjointMethodPNF.

3. Introduced AdjointMethodPNF3D class. This makes it easier to implement
figures of merit which are normalized with respect to power in 3D simulations.
The implementation is very similar to its 2D counterpart. See the
MMI_splitter_3D exmaple for a demonstration of how to use it.

4. Added a gaussian_mode function which calculates the electric and magnetic
field of an approximate fiber mode. This is useful for grating couplers and the
like.

5. ModeFullVector.get_field and ModeFullVector.get_field_interp have been
updated to return non-squeezed fields (i.e. 3D arrays with one dimension equal
to 1). This reduces the need for additional reshaping code when using mode
fields for figures of merit. In order to plot the result of a ModeFullVector
calculation, you can supply the optional squeeze argument to get_field_interp.
See wg_modes_3D.py for an example.

6. Fixed a bug in  the ModeTM solver which caused symmetries to not work
properly.

7. Improved performance of FDFD_TE.update and FDFD_TM.update. In some cases,
the update can be 2x faster. This will speed up gradient calculations in 2D.

8. Changed update box format for 2D solvers to be consistent with 3D solvers
(and to be generally easier to work with). Now, a tuple of real-space bounds
are supplied instead of grid indices. See Silicon_Grating_Coupler for an
example. Note: this will break existing codes built on top of EMopt, but it is
typically a very easy fix.

9. Added emopt.misc.n_SiO2 and emopt.misc.n_Si3N4 material data functions for
convenience.

10. Fixed improper normalization of magnetic fields in gaussian_fields
function.

11. Add support for python 3.x (using futurize). All of the examples have been
tested using python 2.7.5 and python 3.6.6. If you encounter problems with
other version of python, please post a bug report.

12. Updated install.py script for newest SLEPc version.

13. Added complex material support to FDTD solver. In order to use complex
materials, you must set the optional argument complex_eps=True when
instantiating your FDTD object. After doing so, you may assign complex
permittivities to structures. For example, given a rectangle called 'rect', we
can assign it a lossy material value with 'rect.material_value = 2.0 + 1j*0.1'.

Important: currently the source power calculation does not account for the
complex permittivity.

Important 2: currently, negative (real) permittivities are not supported. This is a
result of the fact that dispersion is not currently implemented which is
necessary to make simulation of negative (real) permittivities stable.

Released on March 19, 2019
--------------------------

1. Added function emopt.fomutils.dist_to_edges which calculates gap and bridge
widths in a polygon. This makes it straightforward to impose optimization
constraints which limit the formation of narrow pieces of materials and narrow
gaps.

2. Fixed a problem which has resulted from an mpi4py update which lead to
certain gather operations in adjoint_method throwing an error. This is a
temporary fix and should be revisited in the future.

Released on March 6, 2019
-------------------------

1. Added points_per_bend argument to emopt.geometry.fillet function. This allows you
to explicitly specify a fixed number of points for the rounded corner. This is
important when changing the angle of a corner in an optimization. In retrospect,
specifying the number of points in the fillet is probably preferable to specifying
the number of points per 90 degrees. Consider points_per_90 to be depricated.

2. Changed behavior of FDFD_TE.set_sources() and FDFD_TM.set_sources such that
they accept a mode solver object as an alternative to explicitly specifying the
current density arrays. The FDFD objects will automatically extract the current
density arrays. This is consistent with how the 3D solvers work.

3. Modified 3D examples to use FDTD solver instead of FDFD solver.

4. Big documentation updates: Readthedocs now has tutorials. In particular,
tutorials for:

     a) Setting up and running 2D simulations
     b) Setting up and running 1D mode solutions
     d) Setting up and running 2D mode solutions
     e) Setting up and performing 2D optimizations

have been added. This documentation is still a work in progress and 3D
tutorials will be added soon!

Released on November 16, 2018
------------------------------

1. Bug fixes
   a. Fixed _showwarning definition
   b. Fixed round-off error in DomainCoordinates.copy which led to errors when
   setting the adjoint sources for very particular choices of grid spacing.
   c. Tweaked 1D mode solvers such that build() function actually updates the
   permittivity and permeability. Before this was only being done upon instantiation.

2. Added NURBS class to emopt.geometry. This is a wrapper for the NURBS-python
library which makes it quick and easy to use. NURBS are a good choice of spline for
parameterizing boundaries. In order to use the NURBS class, NURBS-python must be
installed (pip install geomdl).

3. Minor performance increases to gradient calculation. In theory, it can leverage
parallelism a bit better now, but speed improvements may not be significant.

4. Significant performance increase to permittivity grid smoothing calculation in 3D.
For 3D slab structures, a lot of redundant calculation was being performed.
Permittivity values are now cached which can significantly improve speed when
calculating permittivity inside of volumes. This also translates to speedups in
gradient calculations (depending on the problem, speedups could range from ~3x to
~20x).

Released on September 17, 2018
------------------------------

1. Minor bug fixes.

2. Significant speed improvements in 3D FDTD solver.

A number of PML parameters are now precomputed which significantly speeds up the 3D
FDTD solver. You should expect to see an improvement of 2-3x in runtime depending on
the number of cores you use. Currently, this puts EMopt within about a factor of two
or so of certain popular commercial FDTD solvers.

3. Changes to how edge values are shared between processes in the 3D solver. This
change should enable the implementation of periodic/Bloch periodic boundary
conditions in the future. In fact, periodic boundary conditions might already work.

4. Minor documentation improvements to the C++ FDTD code.

Released on July 30, 2018
-------------------------

1. Minor bug fixes

2. "improved" setup.py which helps ease the installation process.

The setup file now looks for a file called '.emopt_deps' in your home directory. This
file defines environment variables which point to the install paths of eigen, boost,
petsc, and slepc. The file has the following format:

    EIGEN_DIR=/path/to/eigen/includes
    BOOST_DIR=/path/to/boost/includes
    PETSC_DIR=/path/to/petsc/installation
    SLEPC_DIR=/path/to/slepc/installation

If you already have emopt installed, simply create the .emopt_deps file in your home
directory:

    $ touch ~/.emopt_deps

And populate the file with the appropriate paths. For example:

    $ echo 'EIGEN_DIR=/opt/include' >> ~/.emopt_deps
    $ echo 'BOOST_DIR=/opt/include' >> ~/.emopt_deps
    $ echo 'PETSC_DIR=/opt/petsc/petsc-3.8.0' >> ~/.emopt_deps
    $ echo 'SLEPC_DIR=/opt/slepc/slepc-3.8.1' >> ~/.emopt_deps

After creating this file, running setup.py should run as it did previously.

The reason for this change is that it has allowed us to introduce an "easy"
installation script which downloads, compiles, and installs the emopt dependencies.
This new script is discussed in the next changelog entry!

3. Added installation script for retrieving and installing the emopt dependencies.

Until now, installing emopt has been a bit of a pain, especially if you are not very
linux savvy. In order to reduce the difficulty of installing emopt, we have
introduced an install.py script which downloads, compiles, and installs the
non-python dependencies.

In order to use this script, a few things may need to be installed first using your
package manager. Once these prerequisites are satisfied, the installation script may
be run:

    $ python install.py

By default, this creates a directory ~/.emopt which contains all of the compiled and
installed dependencies as well as a filed ~/.emopt_deps which tells the emopt
setup.py where to find the installation. Optionally, this script can be run with an
alternative directory in which the dependencies are installed. For example:

    $ python install.py --prefix=~/local/

In this case, the dependencies will be installed to ~/local instead of ~/.emopt. In
either case, the ~/.emopt_deps file will be created (this is currently not something
that can be changed).

The installation documentation has been updated to match this new process. Note that
this is still somewhat experimental, so post an issue if you encounter any errors.

4. Added emopt.geometry which will (hopefully) evolve into a useful library of
subroutines for creating and manipulating geometry. Consider also looking into the
shapely module for additional geometry tools.

Currently, emopt.geometry contains a simple fillet function and a function which
populates existing line segments with more points.

Released on June 15, 2018
-------------------------

1. Introduced new 3D CW-FDTD solver.

This new solver scales to larger/higher resolution problems much better than the
current 3D FDFD solver. This makes it useful for optimizing larger problems which
require high simulation resolutions (i.e. anytime you desire high coupling
efficiencies).

The implementation of this new solver differs a bit from a typical FDTD solver.
Specifically, it uses a ramped CW source to solve for  the frequency domain fields at
a single frequency. This has a number of advantages for using it to do optimization
and has the added benefit of being a bit simpler to implement.

In order to use this solver for optimization, see
examples/MMI_splitter_3D/mmi_1x2_splitter_3D_fdtd.py. The key things to note are that
square grid cells MUST be used (i.e. dx=dy=dz) and symmetry planes can lead to
increased gradient error (in particular, it seems as if fdtd+adjoint method+z
symmetry really struggles to produce accurate results). Depending on your problem,
these constraints may not be a significant issue.

In addition to these optimization constraints, the current CW-FDTD solver does not
support complex materials. This functionality should be added soon, however.

Finally, it is worth noting that this new CW-FDTD solver is not yet totally
optimized. As such, it may run noticeably slower than available commercial solvers.
This is an ongoing effort.

2. Small changes to power normalization in FDFD to make it consistent with new
CW-FDTD solver.

3. Documentation Updates

4. Makefile improvements

Released on May 21, 2018
------------------------

1. Added code needed for most 3D optimizations!

3D optimizations can now be run. For a demonstration see examples/MMI_splitter_3D.
The process for optimizing in 3D is very similar to 2D. There are some minor
differences due to the additional spatial dimension, but philosophically everything
is the same.

There is one caveat: currently AdjointMethodPNF does not support 3D simulations (I
need to come up with a cleaner dimension-agnostic way to handle absorption in the
source power calculation). For non-complex materials, this same functionality can be
achieved using functions provided by emopt.fomutils. See the above examples for
details.

2. Added a function emopt.grid.Polygon.populate_lines which will insert evenly-spaced
points into lists of line segments. This is useful when setting up more complicated
shape optimizations in which you want to manipulate a complicated smooth boundary.

3. Bug fixes, 3D solver performance improvements.

Released on May 11, 2018
------------------------

1. Small bug fixes.

2. Added add_primitives function to StructuredMaterial2D to reduce verbosity. This
function takes a list of Primitive objects.

3. Simple structures can now be imported from klayout! A function called
emopt.io.load_gds_txt has been introduced which allows you to load very simple
text-format GDS files exported from klayout. This provides users with a way to define
their geometry using a more visual process. An example has been added which shows how
to import GDS files.

In order to import GDS files, they must be exported from klayout using File > Save as
and then select Files of type "GDS2 Text files". It is encouraged that you define
different shapes on unique layers so that they are easier to separate out in code.

4. Fixed bugs in ModeFullVector related to generating mode source distributions for
the 3D FDFDsolver. It turns out the coordinate permutations were not being handled
properly which resulted in incorrect mode sources when injecting along x or y (z
worked fine).

5. Major updates to FDFD_3D. 3D solver is pretty much fully functional at this point
with the exception that symmetry boundary conditions have not yet been implemented.
See examples/simple_waveguide/simple_waveguide_3D.py and
examples/simple_waveguide/simple_waveguide_mode.py for examples of how to use it.

Some important notes:
    a. The 3D solver is not quite as straightforward to use as the 2D solver. This is
    a consequence of how the 3D equations are solved (preconditioned iterative solver
    vs direct solver). In particular, the 3D solver uses a multigrid preconditioner
    which solves the problem on coarser grids in order to speed up convergence of the
    simulation. An important parameter is the number of coarser grids to use which
    can be modified using the optional 'mglevels' argument. If not specified, FDFD_3D
    will attempt to automatically choose a good value, but this may not always work.
    If you find that iterations are fast but convergence is slow, try reducing the
    number of mglevels. If the simulation just hangs before getting to the iteration
    step, mglevels should be increased. If confused about this, feel free to drop an
    issue on github!

    b. In theory, the solver should be compatible with adjoint_method.AdjointMethod.
    I haven't really tested this yet, so it may still be buggy. It definitely is not
    compatible with adjoint_method.AdjointMethodFM and AdjointMethodPNF, so stay
    tuned for updates.

    c. Currently the 3D solver is far from the fasts solver out there. Right now, it
    is pretty useable for smaller/coarser simulations. The supplied example which is
    5x3x2.5 um with a 40 nm grid spacing runs in about 50 seconds on two Intel Xeon E5
    v4 using 16 mpi processes. Reducing the mesh size to 30 nm increases the solve
    time to ~2m30s-3m00s.

Released on April 25, 2018
--------------------------

1. Cleaned up 2D FDFD solver build function. Reordered the equations in order to
(potentially) reduce fill-in during LU factorization and also better distribute the
load between processes when updating the matrix. This should improve performance by a
small amount.

2. Added a note about OpenMP threads to readme. Essentially, on many systems, the
default number of threads used for OpenMP will equal the number of cores available.
This is bad when running on a single machine. When running on a single machine, the
number of OpenMP threads should be set such the (# MPI procs)*(# OpenMP threads) <=
(# of cores). Based on experimentation, it looks like OMP_NUM_THREADS=1 with more MPI
proces leads to higher performance on a single machine. On a cluster, however, go
nuts with OMP_NUM_THREADS.

Released on April 24, 2018
--------------------------
1. (in previous commit) Change compile command for OS X in src/Makefile from g++-5 to
g++. This prevented compilation on OS X with up-to-date gcc/g++ installations.

2. Added FieldComponent and SourceComponent "enums" to ensure internal consistency
when retrieving field and source quantities.

3. Changed grid.Material2D.get_values_on to grid.Material2D.get_values_in. Identical
change for grid.Material3D. Furthermore, an optional parameter "squeeze" hass been
added to both Material2D and Material3D which, when set to True, will squeeze out any
length-1 axes from the array. For example, if you want to extract a line from a 2D
Material, then setting squeeze=True will cause the vector which is returned to have
shape (N,). This can be useful for things like plotting.

4. Overhaul of full vector mode solver (now call ModeFullVector). Solver interface has 
been changed a bit and now accepts Material2D/Material3D objects and a DomainCoordinates
to extract the material slice for which the mode is calculated (this change is needed 
for the solver to work properly and it should also make the usage of the mode solve 
less verbose for users). The modes and eigenvalues of the new mode solver have been 
compared against commercial packages and it appears to produce accurate results. At this 
point, the ModeFullVector class should be ready for prime time. 

Furthermore, the ModeFullVector class can now be used to calculate current sources
for exciting waveguide modes in 3D simulations. These waveguide modes can be injected
along any cartesian axis.

5. Renamed Mode_TE to ModeTE and Mode_TM to ModeTM for the sake of internal
consistency.

6. Changed the ModeTE and ModeTM init function to match ModeFullVector. These objects
how accept eps, mu, and a DomainCoordinates and the material distribution is
extracted from those. This, once again, reduces code verbosity a bit and will also
allow us to handle certain internal calculations more effectively (i.e. take
adavantage of grid smoothing better, calculate sources which can inject along either
x or y axis, etc).

7. Fixed some errors in the documentation at the top of modes.py

8. Updated examples according to library changes.

9. Implemented rough 3D solver. Note: This is NOT ready for use, however the
foundation is set. Due to the structure of the 3D equations, a preconditioned
iterative solver is used. The preconditioner combines multigrid with an A^H matrix
preconditioner (which produces makes the system positive semi-definite) and uses
GMRES as a smoother. GCR with restart is used for the outer iteration and MUMPS/LU
factorization is used for the coarse solve. Currently, this works "well enough,"
however it is still very slow compared to FDTD (for example). The 3D solver interface
will be improved in the near future and hopefully be made useable.



Released on March 5, 2018
-------------------------
1. Implemented "general" DomainCoordinates when now replaces LineCoordinates and
PlaneCoordinates. The DomainCoordinates takes the following parameters as input:
xmin, xmax, ymin, ymax, zmin, zmax, dx, dy, dz. All must be specified when creating a
DomainCoordinates object. In order. When all values are unique, the DomainCoordinates
corresponds to a rectangular volume. Planes can be specified by specifying the same
value for two of the min/max parameters. For example, to create a plane normal to the
z-direction located at z=0.5:

    plane = DomainCoordinates(0, W, 0, H, 0.5, 0.5, dx, dy, dz)

In 2D, the exact value for dz is unimportant as long as it is not zero. It is
recommended that 1.0 be used. In order to specify a line, two of the min/max pairs
must have the same values. For example, to create a line along the x direction at
y=0.75, and z=0.5:

    line = DomainCoordinates(0, W, 0.75, 0.75, 0.5, 0.5, dx, dy, dz)

This way of handling slicing sub domains of the simulation allows for a simpler
codebase and also makes slicing domains consistent in 1, 2, and 3 dimensions.

2. Fixed license discrepancy in python documentation

3. Crated new emopt.io module which now houses save_results, load_results, and
plot_iteration

4. Removed emopt.modedata module. Use emopt.modes instead (it is plenty fast).

5. Replaced plot_iterations function with new plot_iteration function (notice that
new function does not end in 's'). The old function was a bit hacky and the results
were not very pleasant to look at. The new function should produce a more informative
plot. Furthermore, it manages plot placement in a more sophisticated way, which
should make it better for structures with different aspect ratios. Finally, plots can
be generated with either a dark or light color scheme (dark is the current default).
This was a personal choice as I generally write code using editors with darker color
schemes, so the new dark scheme is gentler on my eyes :)

6. Refactored code in preparation for 3D code in the future. The following classes
have been renamed:

    emopt.grid.Material -> emopt.grid.Material2D
    emopt.grid.GridMaterial -> emopt.grid.GridMaterial2D
    emopt.grid.StructuredMaterial -> emopt.grid.StructuredMaterial2D

7. The ordering of index ranges in emopt.grid.Material2D.get_values() has been
modified such that it is consistent with get_value. It now takes the form
    
    get_values(k1, k2, j1, j2)

where k1 to k2 is the range of x coordinates (columns in array storage) and j1 to j2
is the range of y coordinates (rows in array storage).

8. A new Material2D structure called ConstantMaterial2D has been added which can
(should) be used for totally uniform structures. This will be faster than using a
StructuredMaterial2D for things like uniform permeability which is often the case in
silicon photonics.

9. All of the examples have been updated and tested.

Released on January 23, 2017 
----------------------------
Note: This is the first entry to the CHANGELOG

1. Fixed boundary condition bugs in mode solver (the problem only manifested if the
permittivity and permeability were uniform across the entire simulation domain and
thus likely did not affect most if any actual simulations.)

2. Changed how PMLs are set. Now, PML is not passed to the constructor of the FDFD
object. Instead, it defaults to a half wavelength everywhere and can be set using the
w_pml property. e.g.

     sim = FDFD_TE(W,H,dx,dy,wlen) 
     sim.w_pml = [1, 1, 1, 1]

3. Introduced ability to change boundary conditions in FDFD_TE and FDFD_TM. This is
accomplished using the "bc" property. The bc property is a 2 character string which
specifies the boundary condition along the x and y directions in the simulation
region. The default value is '00' which corresponds to perfect electrical conductors
on all four edges. The boundary conditions that are currently implemented are

     0 - "perfect electric conductor" : the boundaries of the simulation perfectly
     conduct electrical current.

     M - "perfect Magnetic conductor" : the boundaries of the simulation perfectly
     conduct magnetic current (magnetic monopoles). In most cases, this likely be
     used explicitly by the user.  
     
     E - "Electric field symmetry": the electric field is mirrored across the
     bottom boundary of the simulation domain (y=0).
     
     H - "Magnetic field symmetry": the magnetic field is mirrored across the
     bottom boundary of the simulation domain (y=0). 
     
     P - "Periodic"

Furthermore, in the future, Bloch periodic boundary conditions (B) will be
implemented.

In addition to choosing the desired boundary conditions, the PML widths must be
changed accordingly. For example, if you choose your boundary condition to be '0E'
(PEC on left and right boundaries, E symmetry on bottom boundary, and PEC on top
boundary), then you must be also careful to set the bottom PML width to be 0.0. EMOpt
will throw a warning if you do not do this, so keep an eye out.

Important note: some testing has been done to verify the different boundary
conditions work as intended, however this testing is not exhaustive. 

Examples demonstrating the new boundary conditions can be found in
examples/simple_waveguide/ and examples/periodic_Mie

4. These same boundary conditions have been implemented in the 1D mode solver. The
boundary condition for the mode solver is specified (in most cases) with a length-1
string (i.e. '0', 'M', 'E', 'H', or 'P'). These boundary conditions are designed to
be compatible with the FDFD solver's boundary conditions. In particular, if you
choose E-field symmetry ('E'), for example, the result modes will have an electric
field which is mirrored across the bottom (y=0) boundary. 

5. Fixed bug in source power calculation which occurs when zero-width PMLs are used.