MESTI (Maxwell's Equations Solver with Thousands of Inputs) is an open-source software for full-wave electromagnetic simulations in frequency domain using finite-difference discretization on the Yee lattice.
MESTI implements the augmented partial factorization (APF) method described in this paper. While conventional methods solve Maxwell's equations on every element of the discretization basis set (which contains much more information than is typically needed), APF bypasses such intermediate solution step and directly computes the projected quantities of interest: a generalized scattering matrix given any list of input source profiles and any list of output projection profiles. It can jointly handle thousands of inputs without a loop over them, using fewer computing resources than what a conventional direct method uses to handle a single input. It is exact with no approximation beyond discretization.
MESTI.m here uses MATLAB with double-precision arithmetic and considers 2D systems either in transverse-magnetic (TM) polarization (Hx,Hy,Ez) with
or in transverse-electric (TE) polarization (Ex,Ey,Hz) with
where b(x,y) is the source profile.
We also have a Julia version, MESTI.jl, that supports everything in MESTI.m here plus 3D vectorial capability with anisotropic ε, MPI parallelization, subpixel smoothing for common geometries, and single-precision arithmetic.
MESTI.m is a general-purpose solver with its interface written to provide maximal flexibility. It supports
- TM or TE polarization.
- Any relative permittivity profile ε(x,y), which can be real-valued or complex-valued. The imaginary part of ε(x,y) describes absorption and linear gain.
- Open boundary modeled by a perfectly matched layer (PML) placed on any side(s), with both imaginary-coordinate and real-coordinate stretching (so the PML can accelerate the attenuation of evanescent waves in addition to attenuating the propagating waves).
- Periodic, Bloch periodic, perfect electrical conductor (PEC), and/or perfect magnetic conductor (PMC) boundary conditions.
- Any material dispersion ε(ω), since this is a frequency-domain method.
- Any list of input source profiles (user-specified or automatically built).
- Any list of output projection profiles (or no projection, in which case the complete field profiles are returned).
- Exact outgoing boundaries based on the analytic self-energy in two-sided or one-sided geometries.
- Real-valued or complex-valued frequency ω.
- Automatic or manual choice between APF, a conventional direct solver (e.g., to compute the full field profile), and the recursive Green's function method as the solution method.
- Linear solver using MUMPS (requires installation) or the built-in routines in MATLAB (which uses UMFPACK).
- Shared memory parallelism (with multithreaded BLAS and with OpenMP in MUMPS).
MESTI.m can perform most linear-response computations in 2D and 1D for arbitrary structures, such as
- Scattering problems: transmission, reflection, transport through complex media, waveguide bent, grating coupler, radar cross-section, controlled-source electromagnetic surveys, etc.
- Thermal emission.
- Local density of states.
- Inverse design based on the above quantities.
Since MESTI can use the APF method to handle a large number of input states simultaneously, the computational advantage of MESTI is the most pronounced in multi-input systems.
There are use cases that MESTI can handle but is not necessarily the most efficient, such as
- Broadband response problems involving many frequencies but only a few input states. Time-domain methods like FDTD may be preferred as they can compute a broadband response without looping over the frequencies.
- Problems like plasmonics that require more than an order of magnitude difference in the discretization grid size at different regions of the structure. Finite-element methods may be preferred as they can handle varying spatial resolutions. (Finite-element methods can also adopt APF, but MESTI uses finite difference with a fixed grid size.)
- Homogeneous structures with a small surface-to-volume ratio. Boundary element methods may be preferred as they only discretize the surface.
Problems that MESTI does not handle:
- Nonlinear systems (e.g., χ(2), χ(3), gain media).
- Magnetic systems (e.g., spatially varying permeability μ).
For eigenmode computation, such as waveguide mode solver and photonic band structure computation, one can use mesti_build_fdfd_matrix.m to build the matrix and then compute its eigenmodes. However, we don't currently provide a dedicated function to do so.
No installation is required for MESTI.m itself. To use, simply download it and add the MESTI.m/src folder to the MATLAB search path using the addpath command. The MATLAB version should be R2019b or later. (Using an earlier version is possible but requires minor edits.)
However, to use the APF method, the user needs to install the serial version of MUMPS and its MATLAB interface (note: the serial version of MUMPS already supports multithreading). Without MUMPS, MESTI.m will still run but will only use other methods, which are significantly slower and use much more memory, especially for large systems and systems involving many channels. See this MUMPS installation page for steps to install MUMPS.
The function mesti(syst, B, C, D) provides the most flexibility. Structure syst specifies the polarization to use, permittivity profile, boundary conditions in x and y, which side(s) to put PML with what parameters, the wavelength, and the discretization grid size. Any list of input source profiles can be specified with matrix B, each column of which specifies one source profile b(x,y). Any list of output projection profiles can be specified with matrix C. Matrix D is optional (treated as zero when not specified) and subtracts the baseline contribution; see the APF paper for details.
The function mesti2s(syst, in, out) deals specifically with scattering problems in two-sided or one-sided geometries where ε(x,y) consists of an inhomogeneous scattering region with homogeneous spaces on the left (-x) and right (+x), light is incident from those sides, the boundary condition in x is outgoing, and the boundary condition in y is closed (e.g., periodic or PEC). The user only needs to specify the input and output sides or channel indices or wavefronts through in and out. The function mesti2s() automatically builds the source matrix B, projection matrix C, baseline matrix D, and calls mesti() for the computation.
Flux normalization in x is applied automatically and exactly, so the full scattering matrix is always unitary when ε(x,y) is real-valued.
mesti2s() also offers the additional features of (1) exact outgoing boundaries in x based on the Green's function in free space, and (2) the recursive Green's function method when TM polarization is used; they are efficient for 1D systems and for 2D systems where the width in y is not large.
To compute the complete field profiles, simply omit the argument C or out, or set it to [].
The solution method, the linear solver to use, and other options can be specified with a structure opts as an optional input argument to mesti() or mesti2s(); see documentation for details. They are chosen automatically when not explicitly specified.
The function mesti_build_channels() can be used to build the input and/or output matrices when using mesti() or to determine which channels are of interest when using mesti2s().
Additional functions that build the input/output matrices for different applications and the anisotropic ε(x,y) from subpixel smoothing will be added in the future.
Detailed documentation is given in comments at the beginning of the function files:
For example, typing help mesti in MATLAB brings up the documentation for mesti().
MESTI.m can use shared memory parallelization across the cores within one node/socket through multithreading if MUMPS was compiled with multithreading enabled. The multithreading speed-up comes from using a multithreaded BLAS library and L0-threads feature in MUMPS. Parts of the MUMPS utilize multithreading with L0-threads through OpenMP directives. With APF, most of the computing time is spent on factorization within MUMPS (e.g., see Fig 2d of the APF paper). The factorization and solving stages within MUMPS are parallelized. The building and analyzing stages are not performance-critical and are not parallelized.
L0-threads (see the MUMPS Users' guide for details) enhances the time performance, but marginally increases the memory usage. It is enabled by default. We can change from the default by setting the field opts.use_L0_threads = false in the input argument opts.
The default number of threads is the number of cores available on the machine (either the number of physical cores, or the number of cores requested when running the job with a scheduler like Slurm on a cluster).
We can set the number of threads to be different from the default by setting the environment variable OMP_NUM_THREADS or the field opts.nthreads_OMP of the optional input argument opts.
To check the actual number of threads used in MUMPS, set opts.verbal_solver = true in the input argument and look at the standard output from MUMPS. For example, the following output
executing #MPI = 1 and #OMP = 4
shows that the number of threads used (#OMP) is 4. Since MESTI.m uses the serial version of MUMPS without MPI, #MPI is always 1.
On an Apple Silicon Mac, we typically compile MUMPS by linking to Apple's vecLib for the BLAS library, without multithreading. See details in the MUMPS installation page.
Examples in the examples folder illustrate the usage and the main functionalities of MESTI.m. Each example has its own folder, with its .m script, auxiliary files specific to that example, and a README.md page that shows the example script with its outputs:
- Fabry–Pérot etalon: 1D, using
mesti2s(), with comparison to analytic solution. - Distributed Bragg reflector: 1D, using
mesti2s(), with comparison to analytic solution. - Open channel in a disordered system: 2D, using
mesti2s(), transmission matrix & field profile with customized wavefronts. - Reflection matrix in Gaussian-beam basis: 2D, using
mesti(), reflection matrix in customized basis for a fully open system. - Meta-atom design for metasurfaces: 2D, using
mesti2s()with Bloch periodic boundary. - Angle dependence of a mm-wide metalens: 2D, using
mesti()with compressed input/output matrices (APF-c).
Also see the following repositories:
- APF inverse design: Using MESTI.m to perform multi-channel inverse design.
- Imaging simulations: Using MESTI.m to perform full-wave simulations of imaging inside scattering media.
Here are some animations from the examples above:
- Propagation through a Fabry–Pérot etalon
- Open channel propagating through disorder
- Reflection matrix of a scatterer in Gaussian-beam basis:
- Angle dependence of a mm-wide hyperbolic metalens
For more information on the theory, capability, and benchmarks (e.g., scaling of computing time, memory usage, and accuracy), please see:
- Ho-Chun Lin, Zeyu Wang, and Chia Wei Hsu. Fast multi-source nanophotonic simulations using augmented partial factorization. Nature Computational Science 2, 815–822 (2022).
@article{2022_Lin_NCS,
title = {Fast multi-source nanophotonic simulations using augmented partial factorization},
author = {Lin, Ho-Chun and Wang, Zeyu and Hsu, Chia Wei},
journal = {Nat. Comput. Sci.},
volume = {2},
issue = {12},
pages = {815--822},
year = {2022},
month = {Dec},
doi = {10.1038/s43588-022-00370-6}
}Please cite this paper when you use MESTI.