diff --git a/doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md b/doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md
new file mode 100644
index 000000000..23cf53898
--- /dev/null
+++ b/doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md
@@ -0,0 +1,15 @@
+
+
+
+$$
+ \newcommand{\vb}{\mathbf}
+ \newcommand{\wt}{\widetilde}
+ \newcommand{\mc}{\mathcal}
+ \newcommand{\bmc}[1]{\boldsymbol{\mathcal{#1}}}
+ \newcommand{\sup}[1]{^{\text{#1}}}
+ \newcommand{\sups}[1]{^{\text{#1}}}
+ \newcommand{\sub}[1]{_{\text{#1}}}
+ \newcommand{\subs}[1]{_{\text{#1}}}
+ \newcommand{\pard}[2]{\frac{\partial #1}{\partial #2}}
+ \newcommand{\VMV}[3]{ \Big\langle #1 \Big| #2 \Big| #3 \Big\rangle}
+$$
diff --git a/doc/docs/AdjointSolver/AdjointImplementationNotes.md b/doc/docs/AdjointSolver/AdjointImplementationNotes.md
new file mode 100644
index 000000000..c2f758aaf
--- /dev/null
+++ b/doc/docs/AdjointSolver/AdjointImplementationNotes.md
@@ -0,0 +1,424 @@
+--8<-- "doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md"
+
+# Adjoint-based optimization in meep: Implementation Notes
+
+These notes are intended as something of a companion to the
+[user's manual and tutorial documentation for adjoint-based
+optimization in meep](Overview.md);
+whereas the goal of those pages is to document the user interface
+and explain how to *use* the solver for practical problems,
+our focus here will be what's going on beneath the hood---how
+the solver actually *works.*
+
+Actually, as will be clear to anyone who has ever reviewed the
+fundamentals of
+[adjoint sensitivity analysis](https://en.wikipedia.org/wiki/Adjoint_state_method),
+the theoretical basis of the method and the derivation of its key
+formulas are almost trivially straightforward, with the only
+potential source of difficulty being how to massage the mechanics
+of the procedure into a meep-friendly form.
+
+## Toy problem
+
+Ultimately we will want to use adjoint methods to differentiate
+complex objective functions---involving quantities such as
+Poynting fluxes and mode-expansion coefficients---with respect
+to various different types of parameters describing material geometries.
+However, before tackling the problem in that full generality,
+it's useful to build up to it by starting with a simple toy problem
+and adding complications one at a time. Thus we consider a
+simple waveguide geometry, excited by a point source at
+$\vb{x}\sup{src}$, and define our objective function to be
+simply the frequency-domain electric-field amplitude at
+a point $\vb{x}\sup{obj}$; we will compute the derivative
+of this objective function with respect to the permittivity
+$\epsilon\sup{des}$
+in a small "design region" $\mc{V}\sup{des}$
+centered at a third point $\vb{x}\sup{des}$.
+
+
+
+### Permittivity derivative by finite-differencing
+
+An obvious brute-force way to get at this is simply
+to do two Meep calculations, with $\epsilon\sup{des}$
+augmented by a small finite amount $\Delta\epsilon$ on the
+second run, then compute the difference between the frequency-domain
+electric fields at $\vb{x}\sup{obj}$ and divide
+by $\Delta\epsilon$ to estimate the derivative.
+Here are plots of the unperturbed field and the field perturbation
+(difference between perturbed and unperturbed field):
+
++ $\wt{E_{z0}}(\omega_0, \vb{x})$ (unperturbed):
+
+
+
++ $\Delta {\wt E_z}(\omega_0, \vb{x})$ (perturbed-unperturbed):
+
+
+(Here and below we use the tilde sign ($\sim$) to indicate frequency-domain
+fields and sources.)
+
+### Permittivity derivative from effective sources
+
+One way to think about the effect of a localized permittivity
+bump goes like this: Increasing the permittivity in some localized
+region of a material body corresponds to increasing the
+polarizability in that region---that is, the ease with which
+positive and negative charges in the material, ordinarily bound
+so tightly together that they neutralize each other as sources,
+can be induced by applied electric fields to separate ("polarize"),
+whereupon they cease to cancel each other and act as effective
+sources contributing to electromagnetic fields.
+Of course, if there were no electric field in the material,
+then we could increase its polarizability as much as we pleased
+without producing any sources---zero times a bigger
+coefficient being still zero---but here there *is* a nonzero
+electric field throughout our geometry, due to the point source
+in the unperturbed problem, which means that the effect of bumping the
+permittivity of the design region may be approximated by
+adding new *sources* in that region, with strength
+proportional to $\Delta\epsilon$ and to the unperturbed electric field.
+More specifically, in a frequency-domain problem involving time-harmonic
+fields and sources with angular frequency $\omega$ (time dependence
+$\propto e^{-i\omega t}$), the following perturbations are
+equivalent:
+$$
+\begin{array}{ccc}
+ \left(\begin{array}{c}
+ \text{a permittivity shift of } \Delta\epsilon \\
+ \text{over a small region } \mc{V} \text{ in which} \\
+ \text{the electric field is } \wt{\vb{E}}
+ \end{array}\right)
+ &\Longleftrightarrow &
+ \left(\begin{array}{c}
+ \text{a localized electric current } \\
+ \text{flowing in }\mc{V} \text{ with amplitude } \\
+ \Delta\wt{\vb J}=-i\omega\Delta\epsilon \wt{\vb{E}}
+ \end{array}\right)
+\end{array}
+$$
+
+
+
+Superposing this effective source with the original point source
+at $\vb{x}\sup{src}$ yields a source configuration that,
+acting on the *unperturbed* geometry, produces the same fields
+as the point source alone acting on the *perturbed* geometry.
+
+Alternatively, by exploiting the linearity of Maxwell's equations
+(and assuming we have linear media!) we could just as easily
+remove the original point source and compute the fields of
+$\wt{\Delta \vb{J}}$ *alone*, which, upon dividing through by
+$\Delta\epsilon$, give just the derivatives of field components
+with respect to $\epsilon$. In other words,
+$$ \pard{ \wt{\vb E} (\vb x\sup{obj}) }{\epsilon\sup{des}}
+ \equiv
+ \left(\begin{array}{cc}
+ \text{electric field at }\vb{x}\sup{obj} \text{ due to }\\
+ \text{current at }\vb{x}\sup{des}\text{ with amplitude}\\
+ \wt{\Delta \vb J}=-i\omega\wt{\vb{E}}(\vb{x}\sup{obj})
+ \end{array}\right)
+ \tag{1a}
+$$
+Analogous reasoning yields a prescription for magnetic-field derivatives:
+$$ \pard{\wt{\vb{H}}(\vb x\sup{obj})}{\epsilon\sup{des}}
+ \equiv
+ \left(
+ \begin{array}{cc}
+ \text{magnetic field at }\vb{x}\sup{obj} \text{ due to }\\
+ \text{current at }\vb{x}\sup{des}\text{ with amplitude}\\
+ \wt{\Delta \vb J}=-i\omega\wt{\vb{E}}(\vb{x}\sup{obj})
+ \end{array}
+ \right)
+ \tag{1b}
+$$
+
+### Digression: Configuring time-domain sources for desired frequency-domain fields in Meep
+
+In frequency-domain electromagnetism we usually consider
+a time-harmonic source distribution of the form
+$$
+ \vb{J}\sup{monochromatic}(t,\vb{x})\equiv
+ \wt{\vb{J}}(\vb x)e^{-i\omega t}
+$$
+and we ask for the time-harmonic electric field distribution
+radiated by this distribution:
+$$
+ \vb{E}\sup{monochromatic}(t,\vb{x})\equiv
+ \wt{\vb{E}}(\vb x)e^{-i\omega t}
+$$
+where $\sim$ indicates frequency-domain amplitudes.
+A typical frequency-domain solver might input
+$\wt{\vb J}(\vb x)$ and output
+$\wt{\vb E}(\vb x)$:
+$$ \wt{\vb J}(\vb x)
+ \quad \Longrightarrow \quad
+ \begin{array}{|c|}\hline\\
+ \text{ frequency-domain solver }\\
+ \\\hline\end{array}
+ \quad \Longrightarrow \quad
+ \wt{\vb E}(\vb x)
+$$
+
+On the other hand, when using Meep to compute
+the fields produced by a given spatial source distribution,
+we typically construct a time-domain source of the form
+$$\vb{J}\sup{meep}(t,\vb{x})=G(t)\wt{\vb{J}}(\vb x)$$
+where $G(t)$ is a Gaussian temporal envelope.
+More specifically, for Meep's `GaussianSrc` with
+center frequency $\omega_0=2\pi f_0$,
+frequency width $\Delta \omega =2\pi \Delta f$, and
+peak time $t_0$, we have
+$$ G(t) = e^{-i\omega_0(t-t_0) - \frac{1}{2}[\Delta f(t-t_0)]^2}.$$
+The Fourier transform of this is
+$$
+ \wt G(\omega) \equiv \frac{1}{\sqrt{2\pi}}
+ \int e^{i\omega t}G(t)\,dt =
+ \frac{1}{\Delta f}
+ e^{i\omega t_0 -\frac{(\omega-\omega_0)^2}{2\Delta f^2}}.
+$$
+The Meep version of the above
+input/output diagram looks like
+$$ G(t)\wt{\vb J}(\vb x)
+ \quad \Longrightarrow \quad
+ \begin{array}{|c|}\hline\\
+ \text{ MEEP }\\
+ \text{ (timestepping + DFT) } \\
+ \\\hline\end{array}
+ \quad \Longrightarrow \quad
+ \wt{G}(\omega)\wt{\vb E}(\vb x)
+$$
+
+The upshot is that the frequency-domain fields obtained from a
+Meep run with a Gaussian source
+come out multiplied by a factor of $\wt{G}(\omega)$ that should
+be divided out to yield the desired frequency-domain quantities.
+
+## Invoking reciprocity
+
+It is convenient to describe the process described above
+in the language of frequency-domain Green's functions, which
+express the fields radiated by monochromatic source distributions
+as spatial convolutions:
+$$ \wt{E_i}(\omega, \vb{x}\sup{dest}) =
+ \int
+ \mc{G}\sup{EE}_{ij}(\omega, \vb{x}\sup{dest}, \vb{x}\sup{src})
+ \wt{J_j}(\omega, \vb{x}\sup{src})
+ \,d\vb{x}\sup{src}
+$$
+with $\bmc{G}\sup{EE}$ the
+electric-electric dyadic Green's function of the material geometry
+(giving the electric field produced by a unit-strength electric
+current).
+In this language, the effective-source representation
+of the permittivity derivative reads
+$$ \pard{\wt{E}_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ =
+ \int \mc{G}\sup{EE}_{ij}(\vb{x}\sup{obj}, \vb{x}\sup{des})
+ \left[-i\omega \wt{E}_j(\vb{x}\sup{des})\right]
+ \,d\vb{x}\sup{des}
+$$
+It is convenient to think of the RHS here as a double convolution
+of two vector-valued functions with the $\bmc{G}\sup{EE}$ kernel:
+$$ \pard{\wt{E}_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ =
+ \left[ \vphantom{\wt{\vb E}\sup{des}}
+ \,\, \boldsymbol{\delta}_i\sup{obj}\,\,
+ \right]
+ \star \bmc{G}\sup{EE} \star
+ \left[-i\omega \wt{\vb E}\sup{des}\right]
+$$
+or
+$$ \pard{\wt E_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ =-i\omega\VMV{ \boldsymbol{\delta}_i\sup{obj} }
+ { \bmc G\sup{EE} }
+ { \wt{\vb E}\sup{des} }
+$$
+where $\star$ denotes convolution,
+$\boldsymbol{\delta}_i\sup{obj}$ is short for
+$\delta_{ij} \delta(\vb{x}-\vb{x}\sup{obj}),$
+and the bra-ket notation describes a machine that inputs two
+vector-valued functions $\vb{f},\vb{g}$ and a kernel $\mc{K}$
+and outputs a scalar quantity:
+$$ \VMV{\vb f}{\bmc K}{\vb g}
+ \equiv \sum_{ij} \iint f_i(\vb{x})
+ \mc{K}_{ij}(\vb{x},\vb{x}^\prime)
+ g_j(\vb{x}^\prime)
+ \,d\vb{x} \,d\vb{x}^\prime
+$$
+(Note that this is not a Hermitian inner product, i.e. the first
+factor is not conjugated.)
+
+For the magnetic-field derivative we have similarly
+$$ \pard{\wt{H}_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ =-i\omega\VMV{ \boldsymbol{\delta}_i\sup{obj} }
+ { \boldsymbol{\mc G}\sup{ME} }
+ { \wt{\vb E}\sup{des} }
+$$
+where $\bmc{G}\sup{ME}$ is the magnetic-electric
+Green's function, giving the magnetic field produced
+by an electric current.
+
+Computationally, inner products like
+$\VMV{\vb f}{\bmc G\sup{EE}}{\vb g}$
+for arbitrary functions $\vb{f}(\vb x), \vb{g}(\vb x)$
+may be evaluated in Meep
+as follows:
+
+1. Create an electric current source with
+ spatially-varying amplitude $\vb{g}(\vb x)$
+ and Gaussian temporal envelope $G(t)$.
+
+2. Timestep and DFT to compute the frequency-domain electric field
+ $\wt{\vb E}(\omega; \vb{x})$ produced by this source.
+
+3. Compute the inner product
+ $[\wt{G}(\omega)]^{-1}
+ \int \vb{f}\cdot \wt{\vb E}\,dV.$
+ (The normalization prefactor was discussed above.)
+
+The virtue of writing things this way is that it allows the physical
+property of reciprocity to be expressed as the mathematical property
+that the aforementioned inner-product machine is insensitive to the
+order of its arguments, i.e. we can flip the $\vb f$ and $\vb g$
+inputs and still get the same scalar output:
+$$ \VMV{\vb f}{\bmc K}{\vb g} = \VMV{\vb g}{\bmc K}{\vb f}
+ \quad\text{ for }\quad \bmc K= \bmc{G}\sup{EE}, \bmc{G}\sup{ME}.
+$$
+
+Applying reciprocity to the above expressions for field derivatives yields
+\begin{align*}
+\pard{\wt E_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ &=-i\omega\VMV{ \wt{\vb E }\sup{des} }
+ { \bmc{G}\sup{EE} }
+ { \boldsymbol{\delta}_i\sup{obj} }
+ \tag{2}
+\\[5pt]
+\pard{\wt H_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}
+ &=-i\omega\VMV{ \wt{\vb E }\sup{des} }
+ { \bmc{G}\sup{ME} }
+ { \boldsymbol{\delta}_i\sup{obj} }
+ \tag{3a}
+\\
+ \hphantom{\pard{\wt H_i(\vb{x}\sup{obj})}{\epsilon\sup{des}}}
+ &=+i\omega\VMV{ \wt{\vb E}\sup{des} }
+ { \bmc{G}\sup{EM} }
+ { \boldsymbol{\delta}_i\sup{obj} }
+ \tag{3b}
+\end{align*}
+where in going to the last line we invoked the identity
+$\bmc{G}\sup{EM}=-\bmc{G}\sup{ME}.$
+
+Note that equations (3a) and (3b), notwithstanding their nearly
+identical appearance, describe two rather different
+Meep calculations: In the former case
+we place an electric source at $\vb x\sup{obj}$ and timestep to
+compute the resulting magnetic field, while in the latter
+case we place a magnetic source and timestep
+to compute the resulting electric field. (In both cases,
+upon computing the field in question we proceed to compute its
+overlap with the unperturbed $\vb E$ field in the design region.)
+
+## Differentiating more complicated functions of field components
+
+Thus far we have only considered derivatives of individual
+field components, and then only at a single point $\vb{x}\sup{obj}$;
+more generally, we will want to differentiate functions of
+multiple field components over a subregion of the grid,
+which we will call the *objective region* $\mc{V}\sup{obj}$.
+
+### **E**-field energy in region
+
+As one example, the electric field energy in the objective
+region is defined by an integral over that region, which Meep
+approximate by a weighted sum over grid points:
+$$ \mc{E}=
+ \frac{1}{2}\int_{\mc{V}\sup{obj}} \epsilon |\wt{\vb E}|^2 \,d\mc{V}
+ \approx
+ \frac{1}{2}\sum_{i,\vb{n}\in\mc{V}\sup{obj}} w_{\vb{n}}
+ \epsilon_\vb {n} \wt E_{i\vb n}^* \wt E_{i\vb n}
+$$
+Here the sum is over all field components $i=\{x,y,z\}$ and
+all grid points $\vb{n}$ lying in $\mc{V}\sup{obj}$,
+and $w_{\vb{n}}$ is a cubature weight associated with point $\vb{n}$
+(as returned by [`get_dft_array_metadata`](/Python_User_Interface#metadata)).
+Differentiating, we have
+\begin{align*}
+ \pard{\mc E}{\epsilon\sup{des}}
+ &=\text{Re }\sum_{i\vb{n}\in\mc V\sup{obj}} w_{\vb n}\epsilon_{\vb n}
+ \wt{E}^*_{i\vb n} \pard{\wt{E}_{i\vb n}} {\epsilon\sup{des}}
+\\
+&\hspace{-1.5in}\text{Insert equation (1a):}
+\\
+ &=\text{Re }\left\{ -i\omega \VMV{\epsilon \wt{\vb E}\sup{obj*}}
+ {{\bmc G}\sup{EE}}
+ {\wt{\vb E}\sup{des}}
+ \right\}
+\\
+&\hspace{-1.5in}\text{Invoke reciprocity:}
+\\
+ &=\text{Re }\left\{ -i\omega \VMV{\wt{\vb E}\sup{des}}
+ {{\bmc G}\sup{EE}}
+ {\epsilon \wt{\vb E}\sup{obj*}}
+ \right\}
+\end{align*}
+
+### Poynting flux
+
+A case that arises frequently is that in which the objective region
+is a cross-sectional surface $\mc S\sup{obj}$ cutting normally through a
+waveguide or similar structure and the objective function is the
+normal Poynting flux through $\mc S$.
+For example, the $x$-directed Poynting flux is given by
+$$
+ S_x
+=
+\frac{1}{2}\text{Re }\left\{ \int_{\mc S\sup{obj}} \Big(E^*_y H_z + \cdots\Big) \, d\mathbf{x} \right\}
+\approx \frac{1}{2}\text{Re }\sum_{\vb n\in \mc S\sup{obj}} w_{\vb n}
+ \Big(E^*_{y\vb n} H_{z\vb n} + \cdots\Big)
+$$
+where $\cdots$ refers to three other terms of the form $\pm E^*_i H_j$.
+Differentiating and rearranging slightly, we have
+\begin{align*}
+ \pard{S_x}{\epsilon\sup{des}}
+ &=\text{Re }\sum_{\vb n\in \mc S\sup{obj}} w
+ \left\{ \wt{E}^*_{y} \pard{\wt{H}_z}{\epsilon\sup{des}}
+ +\wt{H}^*_{z} \pard{\wt{E}_y}{\epsilon\sup{des}}
+ +\cdots
+ \right\}
+\\[5pt]
+&\hspace{-0.5in}\text{Use (1a) and (1b):}
+\\[5pt]
+ &=\text{Re }\left\{ -i\omega \VMV{\wt{\vb E}_{y}\sup{obj*}}
+ {\bmc G\sup{ME}}
+ {\wt{\vb E}\sup{des}}
+ -i\omega \VMV{\wt{\vb H}_z\sup{obj*}}
+ {\bmc G\sup{EE}}
+ {\wt{\vb E}\sup{des}}
+ +\cdots
+ \right\}
+\\[5pt]
+&\hspace{-0.5in}\text{Use reciprocity:}
+\\[5pt]
+&=\text{Re }\left\{ -i\omega \VMV{\wt{\vb E}\sup{des}}
+ {\bmc G\sup{ME}}
+ {\wt{\vb E}_y\sup{obj*}}
+ -i\omega \VMV{\wt{\vb E}\sup{des}}
+ {\bmc G\sup{EE}}
+ {\wt{\vb H}_z\sup{obj*}}
+ +\cdots
+ \right\}
+\end{align*}
+
+### Mode coefficient
+
+$$ \alpha_m^\pm = C_1 \pm C_2 $$
+\begin{align*}
+ C_1 &=\frac{1}{\mathcal{N}}
+ \int_{\mc S\sup{obj}} \Big(e^*_y H_z - e^*_z H_y\Big)d\mathbf{x}\\
+ C_2 &=
+ \frac{1}{\mathcal{N}}\int_{\mc S\sup{obj}} \Big(h^*_z E_y - h^*_y E_z\Big)d\mathbf{x}\\
+\end{align*}
+
+--8<-- "doc/docs/AdjointSolver/AdjointLinks.md"
diff --git a/doc/docs/AdjointSolver/ExampleGallery.md b/doc/docs/AdjointSolver/ExampleGallery.md
new file mode 100644
index 000000000..0641ec132
--- /dev/null
+++ b/doc/docs/AdjointSolver/ExampleGallery.md
@@ -0,0 +1,248 @@
+--8<-- "doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md"
+
+---
+# Design optimization with the Meep adjoint solver: Gallery of worked examples
+---
+
+> :bookmark:{.summary .quirklist} **`table of contents`**
+>
+>
+> - [**Full automated optimization of a cross-router device**](#full-automated-optimization-of-a-cross-router-device)
+>
+> An full example of a successful automated design automation in which `mp.adjoint`
+> automatically redesigns a photonic component to improve its performance by
+> several orders of magnitude with no human guidance whatsoever.
+>
+>
+>
+> - [**Numerical validation of adjoint gradients**](#numerical-validation-of-adjoint-gradients)
+>
+> A warm-up validation in which we test the accuracy of objective-function derivatives
+> computed by `mp.adjoint` against known results and finite-difference derivatives
+> for a simple toy model.
+
+
+## **Full automated optimization of a cross-router device**
+
+In this example, `meep.adjoint` automatically designs the
+central section of a [four-way router][CrossRouterExample]
+to achieve near-perfect routing of power from the `west` input
+port to the `north` output port.
+
+The example may be executed as follows:
+
+```bash
+ cd ${MEEP_INSTALLATION}/python/examples/adjoint_optimization
+ python CrossRouter.py --verbose --visualize --optimize
+```
+
+The figures below show results for the 0th (initial), 1st, and final (50th)
+iterations. For each iteration, we show **(1)** the geometry, including the
+central region that is subject to design optimization, **(2)** the frequency-domain
+fields obtained by the forward calculation for the geometry on the given iteration,
+**(3)** the adjoint derivative $\partial f/\partial \epsilon$ indicating
+how the permittivity in the design region should be updated to improve the objective
+function.
+
+Notice in particular that the forward-field results (center diagram) indicate
+that, in the initial geometry, power flows almost directly from the west (input)
+port to the east port; in contrast, by the end of the optimization
+the power flow exhibits a near-perfect 90° bend with almost all power
+flowing out through the north port.
+
+>:bookmark: **Iteration Zero: Initial objective-function value $f_0 = 8.5\times 10^{-4}$**
+
+
+
+
+| Geometry | Forward fields | Adjoint derivative $\partial f /\partial\epsilon$ |
+|:-------------------------------------------:|: -----------------------------------------:|:---------------------------------------------------:|
+|  |  |  |
+
+
+
+>:bookmark: **Iteration One: Objective-function value $f_1 = 3.8\times 10^{-1}$**
+
+| Geometry | Forward fields | Adjoint derivative $\partial f /\partial\epsilon$ |
+|:-------------------------------------------:|: -----------------------------------------:|:---------------------------------------------------:|
+|  |  |  |
+
+
+>:bookmark: **After 50 Iterations: Objective-function value $f_1 = 0.9605$**
+
+
+| Geometry | Forward fields | Adjoint derivative $\partial f /\partial\epsilon$ |
+|:-------------------------------------------:|: -----------------------------------------:|:---------------------------------------------------:|
+|  |  |  |
+
+
+>:bookmark: **Objective-function evolution by iteration**
+ 
+
+
+--------------------------------------------------
+
+## Numerical validation of adjoint gradients
+
+As discussed in the [Reference Manual](ReferenceManual.md),
+the `OptimizationProblem` base class offers command-line options
+for estimating objective-function derivatives with respect
+to individual design variables---that is, individual components
+of the objective-function gradient vector---by numerical finite-differencing.
+This provides a useful check on the accuracy of the
+gradient computed via the adjoint method. In this example
+we'll apply this technique to the simple
+[holey waveguide][HoleyWaveguide] geometry discussed
+in the [Overview][Overview], involving a circular hole
+in an otherwise perfect section of a slab waveguide of
+permittivity $\epsilon\sups{wvg}$.
+For this problem, the design region is the area of the hole,
+and we will consider a particularly simple basis consisting of
+the single basis function $\{b_0(\vb x)\equiv 1\}$, so that
+our sole design variable $\beta_0$ is just the constant permittivity$
+$\epsilon\sups{hole}$ of the hole region; the design problem is
+to tweak $\epsilon\sups{hole}$
+to maximize output power flux (flux through the port labeled 'east')
+for fixed input flux from by the eigenmode source2.
+
+The advantage of this as a validating sanity-check for
+adjoint calculations is the simple dependence
+of the objective on the design variable:
+clearly, transmission through the waveguides can't
+be better than perfect, which is what it is
+when $\epsilon\sups{hole}\equiv \epsilon\sups{wvg}
+and there \textit{is} no hole, so
+the function $f(\epsilon\sup{hole})$ must
+be peaked at $\epsilon\sup{wvg}$ and its
+derivative vanishes there. This gives one
+test on the correctness of an adjoint implementation,
+and we get others by looking at the derivative
+at other points on the curve, with reference values
+be estimated via finite-differencing.
+
+Thus, for a set of 24 $\epsilon\sub{hole}$ values
+ranging from 1 to $2\epsilon\sups{wvg}$ we will
+compute the objective function and its adjoint-based
+and finite-difference derivatives with respect to
+$\epsilon\sups{hole}$.
+
+### Running the calculation serially from a shell script
+
+The individual calculations may be performed by executing the
+[`HoleyWaveguide.py`][`HoleyWaveyguide.py`] script from the
+shell with command-line options, i.e.
+
+```bash
+#!/bin/bash
+
+python HoleyWaveguide.py ${COMMON} --beta 0 1.0 --filebase B1P0
+python HoleyWaveguide.py ${COMMON} --beta 0 1.5 --filebase B1P5
+python HoleyWaveguide.py ${COMMON} --beta 0 2.0 --filebase B120
+...
+python HoleyWaveguide.py ${COMMON} --beta 0 12.0 --filebase BW12P0
+```
+
+Note that the only options that vary from case to case
+are `--beta 0 xx`, setting distinct values for the design
+variable $\beta_0\equiv \epsilon\sups{hole}$,
+and `--filebase,` giving different names to output
+files so they don't overwrite one another.
+The options in common to all cases are:
+
+```bash
+COMMON="--verbose --eval_gradient --fd_index 0 --fd_order 2 --res 40 --df 0.25"
+```
+
+### Running the calculation serially from Python
+
+The effect of running the shell command
+`% python HoleyWaveguide.py OPTIONS`
+can be reproduced from a python script or console
+by creating an instance of `HoleyWaveguide`
+and calling its `run()` method.
+In this case, the command-line options OPTIONS
+should be passed (as a single string separated by spaces)
+to the optional `cmdline` parameter of the `HoleyWaveguide`
+constructor, while `run()` takes no arguments:
+
+```py
+import numpy as np
+
+# generate command lines
+bs = np.linspace(1.0,12.0,23)
+args = " --eval_gradient --fd_order 2 --fd_index 0 "
+cmdlines = [ args + '--beta 0 {} --filebase b{}'.format(b,b) for b in bs ]
+
+# run calculations
+for cmdline in cmdlines:
+ HW=HoleyWaveguide(cmdline=cmdline)
+ HW.run()
+```
+
+### Running the calculation in parallel using `ParallelDesignTester`
+
+Although these calculations aren't particularly time-consuming,
+in a multiprocessor environment
+we can make them go *really* fast by using the `ParallelDesignTester`
+utility distributed with `meep.adjoint`. As described in the
+[ReferenceManual](ReferenceManual#2d-running-many-single-point-calculations-in-parallel-paralleldesigntester),
+this is a simple
+tool provided by `meep.adjoint` for parallel batch processing of jobs
+like the one considered here, which has the advantage of
+presenting as a completely transparent drop-in replacement
+for your existing serial loop:
+the `entirety` of the rejiggering required to parallelize
+your batch is to replace the `for cmdline in cmdlines` loop
+above with the single line
+
+```py
+
+ParallelDesignTester('HoleyWaveguide',cmdlines).launch()
+
+```
+
+This will launch a pool of $N$ server processes to whittle
+its way in parallel through your batch; for each string
+in the `cmdlines` list, one or another server will
+instantiate a `HoleyWaveguide` (or any other class you
+specify) with those command lines, then call its `run()`
+method.
+
+### Results
+
+Whether executed serially or in parallel, each of the 23 jobs
+produces text-based output files named e.g. `b13.out`, `b13.digest,`
+and `b13.legend`, where `b13` is the job-dependent
+label assigned the `filebase` option. The `out`
+file uses a compressed format in which all data
+(objective-function values, adjoint gradient values,
+and finite-difference derivatives) for a specific design
+point (in this case, that means a specific value of
+$\beta_0\equiv \epsilon\sups{hole}$) appear on a single
+line of the file with labels to indicate the design point;
+this is indended for plotting vs. $\beta_0$ or other
+post-processing. The `digest` file reports the same
+data in a more human-readable format. The `legend`
+file contains information helping to interpret the content
+of the other files.
+
+???+ warning "Working directories in `ParallelDesignTester` runs"
+ `ParallelDesignTester` creates and changes to a new
+ timestamped working directory for each pool of parallel
+ batch jobs; look for your output files in
+ subdirectory called something like `HoleyWaveguide_0323.022256,`
+
+For the case at hand, we want to **(a)** concatenate all of the
+`.out` files into an omnibus data file, **(b)** sort by the
+design-variable values reported on specific columns, then
+**(****c)** plot various data quantities (e.g. objective-function
+value, adjoint gradient components, and finite-difference derivatives,
+each of which appear in their own specific columns)
+versus design variable.
+
+
+
+We see that $\epsilon\sup{hole}$ derivatives computed by
+adjoints agree well with numerical finite-difference data.
+
+--8<-- "doc/docs/AdjointSolver/AdjointLinks.md"
diff --git a/doc/docs/AdjointSolver/Overview.md b/doc/docs/AdjointSolver/Overview.md
new file mode 100644
index 000000000..5a1960e51
--- /dev/null
+++ b/doc/docs/AdjointSolver/Overview.md
@@ -0,0 +1,1104 @@
+--8<-- "doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md"
+
+---
+# Adjoint sensitivity analysis for automated design optimization
+---
+
+This section of the Meep documentation
+covers `meep.adjoint`, a submodule of the Meep python module
+that implements an [*adjoint-based sensitivity solver*](https://en.wikipedia.org/wiki/Adjoint_state_method)
+to facilitate automated design optimization via derivative-based numerical optimizers.
+
+> :bookmark:{.tocfloat .summary} **`table of contents`**
+>
+> The `meep.adjoint` documentation is divided into a number of subsections:
+>
+> + This **Overview** page reviews some basic facts about adjoints and optimizers,
+> outlines the steps needed to prepare a Meep
+> geometry for optimization, and sketches the mechanics of
+> the `meep.adjoint` design process.
+> (This page is designed to be a gentle introduction for the
+> adjoint neophyte; experts may want only to skim it before
+> skipping to the next section.)
+>
+> + The [**Reference Manual**](ReferenceManual.md) fills in the details of
+> the topics outlined on this page, spelling out exactly how to
+> write the python script that will drive your optimization
+> session.
+>
+> + The [**Example Gallery**](ExampleGallery.md) presents a number
+> of worked examples that illustrate how `meep.adjoint` tackles
+> practical problems in various settings.
+>
+> + The [**Implementation Notes**](AdjointImplementationNotes.md) page
+> offers a glimpse of what's behind the hood---the physical
+> and mathematical basis of the adjoint method and how they
+> are implemented by `meep.adjoint.` An understanding of this
+> content is not strictly necessary to use the solver, but may
+> help you get more out of the process.
+>
+> + Although logically independent of the adjoint solver,
+> the [**Visualization**](Visualization.md) package bundled
+> with the `meep.adjoint` module offers several general-purpose
+> utilities for convenient visualization of various aspects
+> of Meep calculations, which are
+> useful in *any* meep calculation whether adjoint-related
+> or not.
+
+
+## Overview: Adjoint-based optimization
+
+A common task in electromagnetic engineering is to custom-tune the design
+of some component of a system---a waveguide taper, a power splitter,
+an input coupler, an antenna, etc.---to optimize the performance of the system
+as defined by some problem-specific metric. For our purposes,
+a "design" will consist of a specification of the spatially-varying
+scalar permittivity $\epsilon(\mathbf x)$ in some subregion
+of a Meep geometry, and the performance metric
+will be a physical quantity computed from frequency-domain
+fields---a [power flux][GetFluxes],
+an [energy density][DFTEnergy],
+an
+[eigenmode expansion coefficient][EigenCoefficients],
+or perhaps some mathematical function of one or more of these
+quantities. We will shortly present a smorgasbord of examples; for now,
+perhaps a good one to have in mind is the
+[hole cloak ](#HoleCloak) discussed below, in which a
+chunk of material has been removed from an otherwise perfect waveguide
+section, ruining the otherwise perfectly unidirectional (no scattering or reflection)
+flow of power from a source at one end of the guide to a sink at the other;
+our task is to tweak the permittivity in an annular region
+surrounding the defect (the *cloak*) so as to restore
+as much as possible the reflectionless transfer of power
+across the waveguide---thus hiding or "cloaking"
+the presence of defect from external detection.
+
+Now, given a candidate design
+$\epsilon\sup{trial}(\mathbf{x})$, it's easy enough to see
+how we can use Meep to evaluate
+its performance---just create a
+Meep geometry with $\epsilon\sup{trial}$ as a
+[spatially-varying permittivity function][EpsFunc],
+in the design region,
+add [DFT cells][FluxSpectra]
+to tabulate the frequency-domain Poynting flux entering and departing
+the cloak region,
+[timestep][RunStepFunctions] until
+the DFTs converge, then use post-processing routines like
+[`get_fluxes()`][GetFluxes]
+or perhaps
+[`get_eigenmode_coefficients()`][EigenCoefficients]
+to get the quantities needed to evaluate the performance of the device.
+Thus, for the cost of one full Meep timestepping
+run we obtain the value of our objective function at one point
+in the parameter space of possible inputs.
+
+But *now* what do we do?! The difficulty is that the computation
+izinust described furnishes only the *value* of the objective function
+for a given input, not its *derivatives* with respect to the
+design variables---and thus yields zero insight into how we should
+tweak the design to improve performance.
+In simple cases we might hope to proceed on the basis of physical
+intuition, while
+for small problems with just a few parameters we might try our luck with a
+[derivative-free optimization algorithm](https://en.wikipedia.org/wiki/Derivative-free_optimization);
+however, both of these approaches will run out of steam long before
+we scale up to
+the full complexity of a practical problem with thousands
+of degrees of freedom.
+Alternatively, we could get approximate derivative information by brute-force
+finite-differencing---slightly tweaking one design variable, repeating
+the full timestepping run, and asking how the results changed---but
+proceeding this way to compute derivatives with respect to all $D$
+design variables would require fully $D$ separate timestepping runs;
+for the problem sizes we have in mind, this would make calculating the
+objective-function gradient
+*several thousand times* more costly than calculating its value.
+So we face a dilemma: How can we obtain the derivative information
+necessary for effective optimization in a reasonable amount of time?
+This is where adjoints come to the rescue.
+
+The *adjoint method* of sensitivity analysis is a technique in which
+we exploit certain facts about the physics of a problem and the
+consequent mathematical structure---specifically, in this case, the
+linearity and reciprocity of Maxwell's equations---to rearrange the
+calculation of derivatives in a way that yields an *enormous* speedup
+over the brute-force finite-difference approach. More specifically,
+after we have computed the objective-function value by doing
+the full Meep timestepping run mentioned
+above---the "forward" run in adjoint-method parlance---we can magically
+compute its derivatives with respect to *all* design variables by doing
+just *one* additional timestepping run with a funny-looking choice
+of sources and outputs (the "adjoint" run).
+Thus, whereas gradient computation via finite-differencing is at least $D$
+times more expensive than computing the objective function value,
+with adjoints we get both value and gradient for roughly just *twice* the
+cost of the value alone. Such a bargain! At this modest cost, derivative-based
+optimization becomes entirely feasible.
+
+> :memo:{.center .pct80} **More general materials**
+> Although for simplicity we focus here on case of isotropic,
+> non-magnetic materials, the adjoint solver is also capable
+> of optimizing geometries involving permeable ($\mu\ne 1$) and
+> anisotropic (tensor-valued $\boldsymbol{\epsilon},\boldsymbol{\mu}$) media.
+
+
+---
+
+## Examples of optimization problems
+
+Throughout the `meep.adjoint` documentation we will refer to a running collection of
+simple optimization problems to illustrate the mechanics of optimization,
+among which are the following; click the geometry images to view
+in higher resolution.
+
+
+--------------------------------------------------
+
+
+### The Holey Waveguide
+
+By way of warm-up, a useful toy version of an optimization problem
+is an otherwise pristine length of dielectric slab waveguide in
+which a careless technician has torn a circular `hole` of variable
+permittivity $\epsilon\sup{hole}$.
+
+
+
+
+
+> :bookmark:{.center}
+>
+> 
+
+
+
+
+Incident power from an
+[eigenmode source][EigenModeSource] (cyan line in figure)
+travels leftward through the waveguide, but is partially
+reflected by the hole, resulting in less than 100% power
+the waveguide output (as may be
+characterized in Meep
+by observing power flux and/or
+eigenmode expansion coefficients at the two
+flux monitors, labeled `east` and `west`).
+Our objective is to tweak the value of
+$\epsilon\sup{hole}$ to maximize transmission
+as assessed by one of these metrics.
+The simplicity of this model makes it a useful
+initial warm-up and sanity check for making sure we know
+what we are doing in design optimization; for example,
+[in this worked example][AdjointVsFDTest]
+we use it to confirm the numerical accuracy of
+adjoint-based gradients computed by `mp.adjoint`
+
+
+--------------------------------------------------
+
+### The Hole Cloak
+
+We obtain a more challenging variant of the holey-waveguide problem
+be supposing that the material in the hole region is *not* a
+tunable design parameter---it is fixed at vacuum, say, or
+perfect metal---but that we *are* allowed to vary the permittivity
+in an annular region surrounding the hole in such a way
+as to mimic the effect of filling in the hole, i.e. of hiding
+or "cloaking" the hole as much as possible from external
+ detection.
+
+> :bookmark:{.center}
+>
+> 
+
+For the hole-cloak optimization, the objective function
+will most likely the same as that considered above---namely,
+to maximize the Poynting flux through the flux monitor
+labeled `east` (a quantity we label $S\subs{east}$)
+ or perhaps to maximize the overlap coefficient
+between the actual fields passing through monitor
+``east`` and the fields of (say)
+the $n$th forward- or backward-traveling eigenmode
+of the waveguide (which we label $\{P,M\}_{n,\text{east}}$
+with $P,M$ standing for "plus and minus.")
+On the other hand, the design space here is more
+complicated than for the simple hole, consisting
+of all possible scalar functions $\epsilon(r,\theta)$
+defined on the annular cloak region.
+
+
+--------------------------------------------------
+
+### The cross-router
+
+A different flavor of waveguide-optimization problem arises when we
+consider the *routing* of signals from given inputs to
+given destinations. One example is the *cross-router*, involving
+an intersection between $x-$directed and $y-$directed waveguides,
+with center region of variable permittivity that we may
+tweak to control the routing of power through it.
+
+> :bookmark:{.center}
+>
+> 
+
+Whereas in the previous examples there was more or less
+only one reasonable design objective one might realistically
+want to optimize,
+for a problem like this there are many possibilities.
+For example, given fixed input power supplied by an eigenmode
+source on the "western" branch (cyan line),
+we might be less interested in the absolute output
+power at any port and more concerned with
+achieving maximal *equality* of output
+power among the north, south, and east outputs,
+whereupon we might minimize an objective function of
+the form
+$$f\sub{obj} =
+ \Big( S\sub{north} - S\sub{south}\Big)^2
+ +\Big( S\sub{north} - S\sub{east}\Big)^2
+ + \Big( S\sub{east} - S\sub{south}\Big)^2
+$$
+(or a similar functional form involving eigenmode
+coefficients).
+Alternatively, perhaps we don't care what happens in
+the southern branch, but we really want the fields
+traveling past the `north` monitor
+to have twice as much
+overlap with the forward-traveling 3rd eigenmode of that
+waveguide
+as the `east` fields have with their backward-traveling
+2nd eigenmode:
+
+$$ f\sub{obj} \equiv \Big( P\sub{3,north} - 2M\sub{2,east}\Big)^2$$
+
+The point is that the definition of an optimization problem
+involves not only a set of physical quantities (power fluxes, eigenmode coefficients,
+etc.) that we compute from Meep calculations,
+but also a rule (the objective function $f$) for crunching those
+numbers in some specific way to define a single scalar figure of merit.
+
+In `mp.adjoint` we use the collective term *objective quantities*
+for the power fluxes, eigenmode coefficients, and other physical quantities
+needed to compute the objective function.
+Similarly, the special geometric subregions of
+Meep geometries with
+which objective quantities are associated---the
+cross-sectional flux planes of `DFTFlux` cells or
+field-energy boxes of `DFTField` cells----are known as *objective regions.*
+
+The [Example Gallery][ExampleGallery.md] includes a worked example
+of a full successful iterative optimization in which
+`mp.adjoint` begins with the design shown above and thoroughly rejiggers
+it over the course of 50 iterations to yield a device
+that efficiently routs power around a 90°ree; bend
+from the eigenmode source (cyan line above)
+to the 'north' output port.
+
+--------------------------------------------------
+
+
+
+### The asymmetric splitter
+
+A `splitter` seeks to divide incoming power from one source
+in some specific way among two or more destinations.,
+We will consider an asymmetric splitter in which power
+arriving from a single incoming waveguide is to be routed
+into two outgoing waveguides by varying the design of the
+central coupler region:
+
+> :bookmark:{.center}
+>
+> 
+
+
+--------------------------------------------------
+
+
+## Common elements of optimization geometries: Objective regions, objective functions, design regions, basis sets
+
+The examples above, distinct though they all are, illustrate
+the common defining features that are present in every
+Meep optimization problem:
+
++ **Objective regions:** One or more [regions over which to tabulate frequency-domain fields (DFT cells)][DFTObj]
+ for use in computing power fluxes, mode-expansion coefficients, and other frequency-domain
+ quantities used in characterizing device performance. Because these regions are used to evaluate
+ objective functions, we refer to them as *objective regions.*
+
+>:memo:{.center .pct80} **Objective regions may or may not have zero thickness**
+> In the examples above, it happens that all objective regions are one-dimensional
+> (zero-thickness) flux monitors, indicated by magenta lines; in a 3D geometry they
+> would be two-dimensional flux planes, still of zero thickness in the normal
+> direction. However, objective regions may also be of nonzero thickness, as for
+> instance if the objective function involves the [field energy in a box-shaped
+> subregion of a geometry.][Energy]
+
+ + **Objective quantities and the objective function:**
+ A specification of which quantities (power fluxes, mode coefficients,
+ energies, etc.) are to be computed for each objective region, and how
+ those quantities are to be crunched mathematically to yield a single number
+ measuring device performance. We refer to the individual quantities as
+ *objective quantities*, while the overall function that inputs one more more
+ objective quantities and outputs a single numerical score is the
+ *objective function.*
+
++ **Design region:** A specification of the region over which the material design is to be
+ optimized, i.e. the region in which the permittivity is given by the
+ design quantity $\epsilon\sup{des}(\mathbf x)$.
+ We refer to this as the *design region* $\mathcal{V}\sup{des}$.
+
++ **Basis:** Because the design variable $\epsilon\sup{des}(\mathbf x)$
+ is a continuous function defined throughout a finite volume of space,
+ technically it involves infinitely many degrees of freedom.
+ To yield a finite-dimensional optimization problem, it is convenient
+ to approximate $\epsilon\sup{des}$ as a finite expansion in some
+ convenient set of basis functions, i.e.
+ $$ \epsilon(\mathbf x) \equiv \sum_{d=1}^N \beta_d \mathcal{b}_d(\mathbf x),
+ \qquad \mathbf x\in \mathcal{V}\sup{des},
+ $$
+ where $\{\mathcal{b}_n(\mathbf x)\}$ is a set of $D$ scalar-valued
+ basis functions defined for $\mathbf x\in\mathcal{V}\sup{des}$.
+ The task of the optimizer then becomes to determine
+ numerical values for the $N$-vector of coefficients
+ $\boldsymbol{\beta}=\{\beta_n\},n=1,\cdots,N.$
+
+ For adjoint optimization in Meep, the
+ basis set is chosen by the user, either from among a predefined collection of
+ common basis sets, or as an arbitrary user-defined basis set specified by
+ subclassing an abstract base class in `mp.adjoint.`
+
+## Mechanics of Meep design optimization
+
+With all that by way of background, here's a quick rundown of the
+process you'll follow to optimize a geometry in `meep.adjoint.`
+
+1. You write a python script that implements a subclass of
+ `OptimizationProblem` (an abstract base class in `meep.adjoint`)
+ to describe your specific problem. In particular, your
+ class must override the following two pure virtual methods
+ in `OptimizationProblem:`
+
+??? summary "`init_problem`: One-time initialization"
+ Inputs an `args` structure describing command-line options
+ and returns a 5-tuple
+ ```py
+ fstr, objective_regions, extra_regions, design_region, basis
+ ```
+ defining your objective function, the objective regions on which
+ its inputs (the objective variables) are defined, the design region,
+ and an expansion basis.
+
+
+??? summary "`create_sim`: Instantiation of design-dependent geometries"
+ Inputs a vector of expansion coefficients `beta_vector` and
+ returns a `meep.simulation` describing a geometry with the
+ corresponding spatially-varying permittivity.
+
+
+2. You run computations on your geometry either by executing your
+ script from the shell with command-line options:
+
+```bash
+ % python HoleyWaveguide.py --beta 0 2.3 --eval_gradient
+```
+
+or equivalently from a python script or console by
+calling its `run()` method:
+
+```py
+ from HoleyWaveguide import HoleyWaveguide
+
+ HW=HoleyWaveguide(cmdline='--beta 0 2.3 --eval_gradient')
+ HW.run()
+```
+
+The actual calculations that may be run in this way
+range from a single non-iterative computation of the objective
+function and (optionally) its gradient at a given set of design-variable
+values, to [full-blown iterative design optimization][CrossRouterExample].
+
+Here, in their entirety, are the python scripts implementing the 4 examples
+described above. (These may also be found in the `python/examples/adjoint_optimization`
+subdirectory of your Meep installation.)
+
+
+??? example "`HoleyWaveguide.py`"
+ ```py
+
+ import sys
+ import argparse
+ import numpy as np
+ import meep as mp
+
+ from meep.adjoint import (OptimizationProblem, DFTCell, adjoint_options,
+ xHat, yHat, zHat, origin, FluxLine,
+ ParameterizedDielectric, FourierLegendreBasis)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ class HoleyWaveguide(OptimizationProblem):
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def add_args(self, parser):
+
+ # add new problem-specific arguments
+ parser.add_argument('--dair', type=float, default=-1.0, help='')
+ parser.add_argument('--w_wvg', type=float, default=3.0, help='')
+ parser.add_argument('--eps_wvg', type=float, default=6.0, help='')
+ parser.add_argument('--r_disc', type=float, default=0.5, help='')
+ parser.add_argument('--nr_max', type=int, default=0, help='')
+ parser.add_argument('--kphi_max', type=int, default=0, help='')
+
+ # set problem-specific defaults for existing (general) arguments
+ parser.set_defaults(fcen=0.5)
+ parser.set_defaults(df=0.2)
+ parser.set_defaults(dpml=1.0)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def init_problem(self, args):
+
+ #----------------------------------------
+ # size of computational cell
+ #----------------------------------------
+ lcen = 1.0/args.fcen
+ dpml = 0.5*lcen if args.dpml==-1.0 else args.dpml
+ dair = 0.5*args.w_wvg if args.dair==-1.0 else args.dair
+ L = 3.0*lcen
+ Lmin = 6.0*dpml + 2.0*args.r_disc
+ L = max(L,Lmin)
+ sx = dpml+L+dpml
+ sy = dpml+dair+args.w_wvg+dair+dpml
+ cell_size = mp.Vector3(sx,sy)
+
+ #----------------------------------------
+ #- design region
+ #----------------------------------------
+ design_center = origin
+ design_size = mp.Vector3(2.0*args.r_disc, 2.0*args.r_disc)
+ design_region = mp.Volume(center=design_center, size=design_size)
+
+ #----------------------------------------
+ #- objective regions
+ #----------------------------------------
+ x0_east = args.r_disc + dpml
+ x0_west = -args.r_disc - dpml
+ y0 = 0.0
+ flux_length = 2.0*args.w_wvg
+ east = FluxLine(x0_east,y0,flux_length,mp.X,'east')
+ west = FluxLine(x0_west,y0,flux_length,mp.X,'west')
+
+ objective_regions = [east, west]
+
+ #----------------------------------------
+ #- optional extra regions for visualization
+ #----------------------------------------
+ extra_regions = [mp.Volume(center=origin, size=cell_size)] if args.full_dfts else []
+
+ #----------------------------------------
+ # basis set
+ #----------------------------------------
+ basis = FourierLegendreBasis(radius=args.r_disc, nr_max=args.nr_max, kphi_max=args.kphi_max)
+
+ #----------------------------------------
+ #- source location
+ #----------------------------------------
+ source_center = (x0_west - dpml)*xHat
+ source_size = flux_length*yHat
+
+ #----------------------------------------
+ #- objective function
+ #----------------------------------------
+ fstr='Abs(P1_east)**2+0.0*(P2_east+P1_west+P2_west+M1_east+M2_east+M1_west+M2_west+S_east+S_west)'
+
+ #----------------------------------------
+ #- internal storage for variables needed later
+ #----------------------------------------
+ self.args = args
+ self.dpml = dpml
+ self.cell_size = cell_size
+ self.basis = basis
+ self.design_center = design_center
+ self.source_center = source_center
+ self.source_size = source_size
+
+ return fstr, objective_regions, extra_regions, design_region, basis
+
+ ##############################################################
+ ##############################################################
+ ##############################################################
+ def create_sim(self, beta_vector, vacuum=False):
+
+ args=self.args
+ sx=self.cell_size.x
+
+ wvg=mp.Block(center=origin, material=mp.Medium(epsilon=args.eps_wvg),
+ size=mp.Vector3(self.cell_size.x,args.w_wvg))
+ disc=mp.Cylinder(center=self.design_center, radius=args.r_disc,
+ epsilon_func=ParameterizedDielectric(self.design_center,
+ self.basis,
+ beta_vector))
+
+ geometry=[wvg] if vacuum else [wvg, disc]
+
+ envelope = mp.GaussianSource(args.fcen,fwidth=args.df)
+ amp=1.0
+ if callable(getattr(envelope, "fourier_transform", None)):
+ amp /= envelope.fourier_transform(args.fcen)
+ sources=[mp.EigenModeSource(src=envelope,
+ center=self.source_center,
+ size=self.source_size,
+ eig_band=self.args.source_mode,
+ amplitude=amp
+ )
+ ]
+
+ sim=mp.Simulation(resolution=args.res, cell_size=self.cell_size,
+ boundary_layers=[mp.PML(args.dpml)], geometry=geometry,
+ sources=sources)
+
+ if args.complex_fields:
+ sim.force_complex_fields=True
+
+ return sim
+
+ ######################################################################
+ # if executed as a script, we look at our own filename to figure out
+ # the name of the class above, create an instance of this class called
+ # opt_prob, and call its run() method.
+ ######################################################################
+ if __name__ == '__main__':
+ opt_prob=globals()[__file__.split('/')[-1].split('.')[0]]()
+ opt_prob.run()
+
+ ```
+
+--------------------------------------------------
+
+
+??? example "`HoleCloak.py`"
+ ```py
+
+ import sys
+ import argparse
+ import numpy as np
+ import meep as mp
+
+ from meep.adjoint import (OptimizationProblem, DFTCell, adjoint_options,
+ xHat, yHat, zHat, origin, FluxLine,
+ ParameterizedDielectric, FourierLegendreBasis)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ class HoleCloak(OptimizationProblem):
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def add_args(self, parser):
+
+ # add new problem-specific arguments
+ parser.add_argument('--dair', type=float, default=-1.0, help='')
+ parser.add_argument('--w_wvg', type=float, default=4.0, help='')
+ parser.add_argument('--eps_wvg', type=float, default=6.0, help='')
+ parser.add_argument('--r_disc', type=float, default=0.5, help='')
+ parser.add_argument('--r_cloak', type=float, default=1.5, help='')
+ parser.add_argument('--nr_max', type=int, default=3, help='')
+ parser.add_argument('--kphi_max', type=int, default=2, help='')
+ parser.add_argument('--eps_disc', type=float, default=1.0, help='permittivity in hole region (0.0 for PEC)')
+
+ # set problem-specific defaults for existing (general) arguments
+ parser.set_defaults(fcen=0.5)
+ parser.set_defaults(df=0.2)
+ parser.set_defaults(dpml=1.0)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def init_problem(self, args):
+
+ #----------------------------------------
+ # size of computational cell
+ #----------------------------------------
+ lcen = 1.0/args.fcen
+ dpml = 0.5*lcen if args.dpml==-1.0 else args.dpml
+ dair = 0.5*args.w_wvg if args.dair==-1.0 else args.dair
+ L = 3.0*lcen
+ Lmin = 6.0*dpml + 2.0*args.r_cloak
+ L = max(L,Lmin)
+ sx = dpml+L+dpml
+ sy = dpml+dair+args.w_wvg+dair+dpml
+ cell_size = mp.Vector3(sx, sy, 0.0)
+
+ #----------------------------------------
+ #- design region
+ #----------------------------------------
+ design_center = origin
+ design_size = mp.Vector3(2.0*args.r_cloak, 2.0*args.r_cloak)
+ design_region = mp.Volume(center=design_center, size=design_size)
+
+ #----------------------------------------
+ #- objective regions
+ #----------------------------------------
+ fluxW_center = (+args.r_cloak+ dpml)*xHat
+ fluxE_center = (-args.r_cloak- dpml)*xHat
+ flux_size = 2.0*args.w_wvg*yHat
+
+ #fluxW_region = mp.FluxRegion(center=fluxW_center, size=flux_size, direction=mp.X)
+ #fluxE_region = mp.FluxRegion(center=fluxE_center, size=flux_size, direction=mp.X)
+ x0_east = args.r_cloak + dpml
+ x0_west = -args.r_cloak - dpml
+ y0 = 0.0
+ flux_length = 2.0*args.w_wvg
+ east = FluxLine(x0_east,y0,flux_length,mp.X,'east')
+ west = FluxLine(x0_west,y0,flux_length,mp.X,'west')
+
+ objective_regions = [east, west]
+
+ #----------------------------------------
+ #- optional extra regions for visualization
+ #----------------------------------------
+ extra_regions = [mp.Volume(center=origin, size=cell_size)] if args.full_dfts else []
+
+ #----------------------------------------
+ # basis set
+ #----------------------------------------
+ basis = FourierLegendreBasis(outer_radius=args.r_cloak, inner_radius=args.r_disc,
+ nr_max=args.nr_max, kphi_max=args.kphi_max)
+
+ #----------------------------------------
+ #- source location
+ #----------------------------------------
+ source_center = (x0_west-dpml)*xHat
+ source_size = flux_length*yHat
+
+ #----------------------------------------
+ #- objective function
+ #----------------------------------------
+ fstr='Abs(P1_east)**2+0.0*(P1_west + M1_east + M1_west + S_west + S_east)'
+
+ #----------------------------------------
+ #- internal storage for variables needed later
+ #----------------------------------------
+ self.args = args
+ self.dpml = dpml
+ self.cell_size = cell_size
+ self.basis = basis
+ self.design_center = design_center
+ self.source_center = source_center
+ self.source_size = source_size
+
+ return fstr, objective_regions, extra_regions, design_region, basis
+
+ ##############################################################
+ ##############################################################
+ ##############################################################
+ def create_sim(self, beta_vector, vacuum=False):
+
+ args=self.args
+ sx=self.cell_size.x
+
+ wvg=mp.Block(center=origin, material=mp.Medium(epsilon=args.eps_wvg),
+ size=mp.Vector3(self.cell_size.x,args.w_wvg))
+ cloak=mp.Cylinder(center=self.design_center, radius=args.r_cloak,
+ epsilon_func=ParameterizedDielectric(self.design_center,
+ self.basis,
+ beta_vector))
+ disc=mp.Cylinder(center=self.design_center, radius=args.r_disc,
+ material=(mp.metal if args.eps_disc==0 else
+ mp.Medium(epsilon=args.eps_disc)))
+
+ geometry=[wvg] if vacuum else [wvg, cloak, disc]
+
+ envelope = mp.GaussianSource(args.fcen,fwidth=args.df)
+ amp=1.0
+ if callable(getattr(envelope, "fourier_transform", None)):
+ amp /= envelope.fourier_transform(args.fcen)
+ sources=[mp.EigenModeSource(src=envelope,
+ center=self.source_center,
+ size=self.source_size,
+ eig_band=self.args.source_mode,
+ amplitude=amp
+ )
+ ]
+
+ sim=mp.Simulation(resolution=args.res, cell_size=self.cell_size,
+ boundary_layers=[mp.PML(args.dpml)], geometry=geometry,
+ sources=sources)
+
+ if args.complex_fields:
+ sim.force_complex_fields=True
+
+ return sim
+
+ ######################################################################
+ # if executed as a script, we look at our own filename to figure out
+ # the name of the class above, create an instance of this class called
+ # opt_prob, and call its run() method.
+ ######################################################################
+ if __name__ == '__main__':
+ opt_prob=globals()[__file__.split('/')[-1].split('.')[0]]()
+ opt_prob.run()
+ ```
+
+??? example "`CrossRouter.py`"
+
+ ```py
+
+ import numpy as np
+ import meep as mp
+
+ from meep.adjoint import (OptimizationProblem, FluxLine,
+ xHat, yHat, zHat, origin,
+ ParameterizedDielectric, FiniteElementBasis)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ class CrossRouter(OptimizationProblem):
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def add_args(self, parser):
+
+ # add new problem-specific arguments
+ parser.add_argument('--wh', type=float, default=1.5, help='width of horizontal waveguide')
+ parser.add_argument('--wv', type=float, default=1.5, help='width of vertical waveguide')
+ parser.add_argument('--l_stub', type=float, default=3.0, help='waveguide input/output stub length')
+ parser.add_argument('--eps', type=float, default=6.0, help='waveguide permittivity')
+ parser.add_argument('--r_design', type=float, default=0.0, help='design region radius')
+ parser.add_argument('--l_design', type=float, default=4.0, help='design region side length')
+ parser.add_argument('--nfe', type=int, default=2, help='number of finite elements per unit length')
+ parser.add_argument('--n_weight', type=float, default=1.00, help='')
+ parser.add_argument('--s_weight', type=float, default=0.00, help='')
+ parser.add_argument('--e_weight', type=float, default=0.00, help='')
+
+ # set problem-specific defaults for existing (general) arguments
+ parser.set_defaults(fcen=0.5)
+ parser.set_defaults(df=0.2)
+ parser.set_defaults(dpml=1.0)
+ parser.set_defaults(epsilon_design=6.0)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def init_problem(self, args):
+
+ #----------------------------------------
+ # size of computational cell
+ #----------------------------------------
+ lcen = 1.0/args.fcen
+ dpml = 0.5*lcen if args.dpml == -1.0 else args.dpml
+ design_length = 2.0*args.r_design if args.r_design > 0.0 else args.l_design
+ sx = sy = dpml + args.l_stub + design_length + args.l_stub + dpml
+ cell_size = mp.Vector3(sx, sy, 0.0)
+
+ #----------------------------------------
+ #- design region bounding box
+ #----------------------------------------
+ design_center = origin
+ design_size = mp.Vector3(design_length, design_length)
+ design_region = mp.Volume(center=design_center, size=design_size)
+
+ #----------------------------------------
+ #- objective and source regions
+ #----------------------------------------
+ gap = args.l_stub/6.0 # gap between source region and flux monitor
+ d_flux = 0.5*(design_length + args.l_stub) # distance from origin to NSEW flux monitors
+ d_source = d_flux + gap # distance from origin to source
+ d_flx2 = d_flux + 2.0*gap
+ l_flux_NS = 2.0*args.wv
+ l_flux_EW = 2.0*args.wh
+ north = FluxLine(0.0, +d_flux, l_flux_NS, mp.Y, 'north')
+ south = FluxLine(0.0, -d_flux, l_flux_NS, mp.Y, 'south')
+ east = FluxLine(+d_flux, 0.0, l_flux_EW, mp.X, 'east')
+ west1 = FluxLine(-d_flux, 0.0, l_flux_EW, mp.X, 'west1')
+ west2 = FluxLine(-d_flx2, 0.0, l_flux_EW, mp.X, 'west2')
+
+ objective_regions = [north, south, east, west1, west2]
+
+ source_center = mp.Vector3(-d_source, 0.0)
+ source_size = mp.Vector3(0.0,l_flux_EW)
+
+ #----------------------------------------
+ #- optional extra regions for visualization
+ #----------------------------------------
+ extra_regions = [mp.Volume(center=origin, size=cell_size)] if args.full_dfts else []
+
+ #----------------------------------------
+ # basis set
+ #----------------------------------------
+ basis = FiniteElementBasis(lx=args.l_design, ly=args.l_design, density=args.nfe)
+
+ #----------------------------------------
+ #- objective function
+ #----------------------------------------
+ fstr=( ' {:s}*Abs(P1_north)**2'.format('0.0' if args.n_weight==0.0 else '{}'.format(args.n_weight))
+ + ' + {:s}*Abs(M1_south)**2'.format('0.0' if args.s_weight==0.0 else '{}'.format(args.s_weight))
+ + ' + {:s}*Abs(P1_east)**2'.format('0.0' if args.e_weight==0.0 else '{}'.format(args.e_weight))
+ + ' + 0.0*(P1_north + M1_south + P1_east + P1_west1 + P1_west2)'
+ + ' + 0.0*(M1_north + M1_south + M1_east + M1_west1 + M1_west2)'
+ + ' + 0.0*(S_north + S_south + S_east + S_west1 + S_west2)'
+ )
+
+ #----------------------------------------
+ #- internal storage for variables needed later
+ #----------------------------------------
+ self.args = args
+ self.dpml = dpml
+ self.cell_size = cell_size
+ self.basis = basis
+ self.design_center = design_center
+ self.design_size = design_size
+ self.source_center = source_center
+ self.source_size = source_size
+
+ if args.eps_design is None:
+ args.eps_design = args.eps
+
+ return fstr, objective_regions, extra_regions, design_region, basis
+
+ ##############################################################
+ ##############################################################
+ ##############################################################
+ def create_sim(self, beta_vector, vacuum=False):
+
+ args=self.args
+
+ hwvg=mp.Block(center=origin, material=mp.Medium(epsilon=args.eps),
+ size=mp.Vector3(self.cell_size.x,args.wh))
+ vwvg=mp.Block(center=origin, material=mp.Medium(epsilon=args.eps),
+ size=mp.Vector3(args.wv,self.cell_size.y))
+
+ if args.r_design>0.0:
+ router=mp.Cylinder(center=self.design_center, radius=args.r_design,
+ epsilon_func=ParameterizedDielectric(self.design_center,
+ self.basis,
+ beta_vector))
+ else:
+ router=mp.Block(center=self.design_center, size=self.design_size,
+ epsilon_func=ParameterizedDielectric(self.design_center,
+ self.basis,
+ beta_vector))
+ geometry=[hwvg, vwvg, router]
+
+ envelope = mp.GaussianSource(args.fcen,fwidth=args.df)
+ amp=1.0
+ if callable(getattr(envelope, "fourier_transform", None)):
+ amp /= envelope.fourier_transform(args.fcen)
+ sources=[mp.EigenModeSource(src=envelope,
+ center=self.source_center,
+ size=self.source_size,
+ eig_band=args.source_mode,
+ amplitude=amp
+ )
+ ]
+
+ sim=mp.Simulation(resolution=args.res, cell_size=self.cell_size,
+ boundary_layers=[mp.PML(self.dpml)], geometry=geometry,
+ sources=sources)
+
+ if args.complex_fields:
+ sim.force_complex_fields=True
+
+ return sim
+
+ ######################################################################
+ # if executed as a script, we look at our own filename to figure out
+ # the name of the class above, create an instance of this class called
+ # op, and call its run() method.
+ ######################################################################
+ if __name__ == '__main__':
+ op=globals()[__file__.split('/')[-1].split('.')[0]]()
+ op.run()
+
+ ```
+
+??? example "`AsymmetricSplitter.py`"
+ ```py
+ import sys
+ import argparse
+ import numpy as np
+ import meep as mp
+
+ from meep.adjoint import (OptimizationProblem, DFTCell, adjoint_options,
+ xHat, yHat, zHat, origin, FluxLine,
+ ParameterizedDielectric, FiniteElementBasis)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ class AsymmetricSplitter(OptimizationProblem):
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def add_args(self, parser):
+
+ # add new problem-specific arguments
+ parser.add_argument('--dair', type=float, default=-1.0, help='')
+ parser.add_argument('--w_in', type=float, default=1.0, help='width of input waveguide')
+ parser.add_argument('--w_out1', type=float, default=0.5, help='width of output waveguide 1')
+ parser.add_argument('--w_out2', type=float, default=0.5, help='width of output waveguide 2')
+ parser.add_argument('--l_stub', type=float, default=3.0, help='length of waveguide input/output stub')
+ parser.add_argument('--l_design', type=float, default=2.0, help='length of design region')
+ parser.add_argument('--h_design', type=float, default=6.0, help='height of design region')
+ parser.add_argument('--eps_in', type=float, default=6.0, help='input waveguide permittivity')
+ parser.add_argument('--eps_out1', type=float, default=2.0, help='output waveguide 1 permittivity')
+ parser.add_argument('--eps_out2', type=float, default=12.0, help='output waveguide 2 permittivity')
+ parser.add_argument('--nfe', type=int, default=2, help='number of finite elements per unit length')
+
+ # set problem-specific defaults for existing (general) arguments
+ parser.set_defaults(fcen=0.5)
+ parser.set_defaults(df=0.2)
+ parser.set_defaults(dpml=1.0)
+
+ ##################################################
+ ##################################################
+ ##################################################
+ def init_problem(self, args):
+
+ #----------------------------------------
+ # size of computational cell
+ #----------------------------------------
+ lcen = 1.0/args.fcen
+ dpml = 0.5*lcen if args.dpml==-1.0 else args.dpml
+ dair = 0.5*args.w_in if args.dair==-1.0 else args.dair
+ sx = dpml + args.l_stub + args.l_design + args.l_stub + dpml
+ sy = dpml + dair + args.h_design + dair + dpml
+ cell_size = mp.Vector3(sx, sy, 0.0)
+
+ #----------------------------------------
+ #- design region
+ #----------------------------------------
+ design_center = origin
+ design_size = mp.Vector3(args.l_design, args.h_design, 0.0)
+ design_region = mp.Volume(center=design_center, size=design_size)
+
+ #----------------------------------------
+ #- objective regions
+ #----------------------------------------
+ x_in = -0.5*(args.l_design + args.l_stub)
+ x_out = +0.5*(args.l_design + args.l_stub)
+ y_out1 = +0.25*args.h_design
+ y_out2 = -0.25*args.h_design
+
+ flux_in = FluxLine(x_in, 0.0, 2.0*args.w_in, mp.X, 'in')
+ flux_out1 = FluxLine(x_out, y_out1, 2.0*args.w_out1, mp.X, 'out1')
+ flux_out2 = FluxLine(x_out, y_out2, 2.0*args.w_out2, mp.X, 'out2')
+
+ objective_regions = [flux_in, flux_out1, flux_out2]
+
+ #----------------------------------------
+ #- optional extra regions for visualization if the --full-dfts options is present.
+ #----------------------------------------
+ extra_regions = [mp.Volume(center=origin, size=cell_size)] if args.full_dfts else []
+
+ #----------------------------------------
+ # forward source region
+ #----------------------------------------
+ source_center = (x_in - 0.25*args.l_stub)*xHat
+ source_size = 2.0*args.w_in*yHat
+
+ #----------------------------------------
+ # basis set
+ #----------------------------------------
+ basis = FiniteElementBasis(args.l_design, args.h_design, args.nfe)
+
+ #----------------------------------------
+ #- objective function
+ #----------------------------------------
+ fstr = ( 'Abs(P1_out1)**2'
+ + '+0.0*(P1_out1 + M1_out1)'
+ + '+0.0*(P1_out2 + M1_out2)'
+ + '+0.0*(P1_in + M1_in + S_out1 + S_out2 + S_in)'
+ )
+
+ #----------------------------------------
+ #- internal storage for variables needed later
+ #----------------------------------------
+ self.args = args
+ self.dpml = dpml
+ self.cell_size = cell_size
+ self.basis = basis
+ self.design_center = origin
+ self.design_size = design_size
+ self.source_center = source_center
+ self.source_size = source_size
+
+ return fstr, objective_regions, extra_regions, design_region, basis
+
+ ##############################################################
+ ##############################################################
+ ##############################################################
+ def create_sim(self, beta_vector, vacuum=False):
+
+ args=self.args
+ sx=self.cell_size.x
+
+ x_in = -0.5*(args.l_design + args.l_stub)
+ x_out = +0.5*(args.l_design + args.l_stub)
+ y_out1 = +0.25*args.h_design
+ y_out2 = -0.25*args.h_design
+
+ wvg_in = mp.Block( center=mp.Vector3(x_in,0.0),
+ size=mp.Vector3(args.l_stub,args.w_in),
+ material=mp.Medium(epsilon=args.eps_in))
+ wvg_out1 = mp.Block( center=mp.Vector3(x_out,y_out1),
+ size=mp.Vector3(args.l_stub,args.w_out1),
+ material=mp.Medium(epsilon=args.eps_out1))
+ wvg_out2 = mp.Block( center=mp.Vector3(x_out,y_out2),
+ size=mp.Vector3(args.l_stub,args.w_out2),
+ material=mp.Medium(epsilon=args.eps_out2))
+ design = mp.Block( center=origin,
+ size=mp.Vector3(args.l_design,args.h_design),
+ epsilon_func=ParameterizedDielectric(self.design_center,
+ self.basis,
+ beta_vector)
+ )
+
+ geometry=[wvg_in, wvg_out1, wvg_out2, design]
+
+ envelope = mp.GaussianSource(args.fcen,fwidth=args.df)
+ amp=1.0
+ if callable(getattr(envelope, "fourier_transform", None)):
+ amp /= envelope.fourier_transform(args.fcen)
+ sources=[mp.EigenModeSource(src=envelope,
+ center=self.source_center,
+ size=self.source_size,
+ eig_band=self.args.source_mode,
+ amplitude=amp
+ )
+ ]
+
+ sim=mp.Simulation(resolution=args.res, cell_size=self.cell_size,
+ boundary_layers=[mp.PML(args.dpml)], geometry=geometry,
+ sources=sources)
+
+ if args.complex_fields:
+ sim.force_complex_fields=True
+
+ return sim
+
+ ######################################################################
+ # if executed as a script, we look at our own filename to figure out
+ # the name of the class above, create an instance of this class called
+ # opt_prob, and call its run() method.
+ ######################################################################
+ if __name__ == '__main__':
+ opt_prob=globals()[__file__.split('/')[-1].split('.')[0]]()
+ opt_prob.run()
+
+ ```
+
+
+
+--8<-- "doc/docs/AdjointSolver/AdjointLinks.md"
diff --git a/doc/docs/AdjointSolver/ReferenceManual.md b/doc/docs/AdjointSolver/ReferenceManual.md
new file mode 100644
index 000000000..267419343
--- /dev/null
+++ b/doc/docs/AdjointSolver/ReferenceManual.md
@@ -0,0 +1,279 @@
+--8<-- "doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md"
+
+---
+# Design optimization with the Meep adjoint solver: A reference manual
+---
+
+As described in the [Overview](Overview.md), the first step in
+using `mp.adjoint` is to write a python script implementing
+a subclass of the `OptimizationProblem` abstract base class
+defined by `mp.adjoint`. Once that's ready, you can run
+your script with command-line options instructing `mp.adjoint`
+to carry out various calculations; these may be *single-point*
+calculations, in which the geometry is fixed at given set
+of values you specify for the design variables (thus defining
+a single point in the space of possible input), or full-blown
+*iterative optimizations* in which `mp.adjoint` automatically
+evolves the design variables toward the values that optimize
+your objective.
+
+> :bookmark:{.summary} **`table of contents`**
+>
+> 1. **Defining your problem**: Writing a python script for `mp.adjoint`
+>
+> a. The `OptimizationProblem` abstract base class
+>
+> b. Mandatory class-method overrides: `init_problem` and `create_sim`
+>
+> c. Optional class-method override: `add_args`
+>
+> 2. **Running your problem**: Visualization, single-point calculations, and full iterative optimization
+>
+> a. Visualizing your geometry
+>
+> b. Evaluating the objective function value and gradient at a single design point
+>
+> c. Testing gradient components by finite-differencing
+>
+> d. Running many single-point calculations in parallel:`ParallelDesignTester`
+
+> e. Full iterative optimization
+>
+> [4. Built-in command-line options](#4-built-in-command-line-options)
+
+## 1. Defining your problem: Writing a python script for `mp.adjoint`
+
+### 1a. The `OptimizationProblem` abstract base class
+
+As described in the [Overview](Overview.md), the python script
+that drives your `mp.adjoint` session implements a subclass
+of the `OptimizationProblem` base class defined by `mp.adjoint.`
+This is a high-level abstraction of the design-automation
+process; the base class knows how to do various general things
+involving Meep geometries and objective
+functions, but is lacking crucial information from you,
+without which it can't do anything on its own.
+That is to say, `OptimizationProblem` is an
+[abstract base class](https://docs.python.org/3/library/abc.html)
+with two pure virtual methods that your derived class
+must override to describe the specifics of your design problem.
+
+### 1b. Mandatory class-method overrides
+
+More specifically, your subclass of `OptimizationProblem`
+must furnish implementations of the following two pure virtual
+methods left unimplemented in the base class. (Click the header
+bars to unfold the description of each item.)
+
+
+??? summary "`init_problem`: One-time initialization"
+ Your `init_problem` routine will be called once, at the beginning of a `mp.adjoint`
+ session; think of it as the class constructor. (Indeed, it is called from the
+ parent class constructor.)
+ It has two purposes: **(a)** to give you a chance to complete any one-time initialization
+ tasks you need to do, and **(b)** to communication to the base class the
+ [CommonElementsOfOptimizationGeometries](Overview.md#common-elements-of-optimization-geometries-objective-regions-objective-functions-design-regions-basis-sets)
+ needed to turn a generic Meep simulation into an optimization problem.
+
+ - **Calling convention:** `def init_problem(self,args)`
+
+ - `args`: Structure containing values of all command-line options.
+
+ - **Return values:**
+
+ The routine should return a 5-tuple
+
+ ```py
+ fstr, objective_regions, extra_regions, design_region, basis
+ ```
+
+ where
+
+ - `fstr` is a string specifying your objective function
+
+ - `objective_regions` is a list of regions over which to compute
+ frequency-domain (DFT) fields needed to evaluate the quantities
+ on which your objective function depends
+
+ - `extra_regions` is a list of additional regions over which to compute
+ DFT fields for post-processing or visualization; it will often be
+ just the empty list `[]`
+
+ - `design_region` is a region encompassing the variable-permittivity
+ region of your geometry
+
+ - `basis` is a specification of the set of basis functions used to
+ expand the permittivity. This is a subclass of the `Basis` base class
+ defined by `meep.adjoint`; you may implement your own arbitrary
+ basis, or use one of several predefined bases provided by `meep.adjoint.`
+
+---
+
+??? summary "`create_sim`: Instantiating a geometry with given design variables"
+ Your `create_sim` routine will be called each time `mp.adjoint` needs to
+ compute your objective function for a particular set of design-variable
+ values.
+
+--------------------------------------------------
+
+### 1c. Optional class-method override
+
+??? summary "`add_args`: Configure problem-specific command-line arguments"
+ The `OptimizationProblem` base class defines a
+ [large number of general-purpose command-line options](#4-built-in-command-line-options),
+ with judiciously chosen default values, to control `mp.adjoint` calculations.
+ In many cases you will want **(1)** to reconfigure the default values as
+ appropriate for your problem, **(2)** to add additional options relevant to
+ your specific geometry. Your subclass can do both of these things by overriding
+ the `add_args` class method, which will be called just before the actual
+ command-line arguments are parsed.
+
+ **Prototype:** `init_args(self,parser)`
+
+ - `parser`: [`argparse`](https://docs.python.org/3/library/argparse.html)
+
+ A structure that has been initialized with
+ [all built-in `mp.adjoint` options](#4-built-in-command-line-options)
+ and their default values.
+ You may call `parser.set_defaults` to change the default values of
+ built-in options and `parser.add_argument` to add new options.
+ (The actual values specified for all arguments are made available
+ to you via the `args` parameter passed to your `init_problem`
+ routine.)
+
+ **Return values:** None.
+
+--------------------------------------------------
+
+### Objective-function specifications
+
+Your objective function is specified by the string `fstr` that you return as
+the first element in the 5-tuple return value of `init_problem.`
+
+Your objective function will depend on one or more objective quantities,
+such as power fluxes or eigenmode expansion coefficients, associated
+with specific objective regions in your geometry. `mp.adjoint` defines
+a convention for assigning a unique character string to each
+objective function in your geometry. Your `fstr` should use
+these labels to refer to the objective quantities on which it
+depends.
+
+>:bookmark: **Naming convention for objective quantities**
+>
+> | Quantity | Label |
+ |:--------------------------------------------------------------------------------:|:------------|
+ | Power flux through flux region `r` | `S_r` |
+ | Expansion coefficient for forward-travelling eigenmode *n* at flux region `r` | `Pn_r` |
+ | Expansion coefficient for backward-travelling eigenmode *n* at flux region `r` | `Mn_r` |
+>
+> Note that the flux-region label `r` may be a character string like `east` if your `init_problem`
+ method assigned a name to the `DFTCell` for that region. Otherwise, `r` is an integer-valued
+ index corresponding to the zero-based index of the DFT cell in the `objective_regions` list
+ returned by your `init_problem`.
+
+### Expansion bases
+
+The `basis` field returned by `init_problem` is a subclass of the
+`Basis` base class implemented by `meep.adjoint.` You can write
+your own subclass to define an arbitrary basis, or use one of the
+built-in basis sets provided by `meep.adjoint.` A good default
+choice is `FiniteElementBasis`:
+
+```py
+ basis = FiniteElementBasis(lx=4, ly=4, density=4)
+```
+
+which defines a basis of localized
+functions over a rectangle of dimensions `l_x`×`l_y` with `density`
+elements per unit length.
+
+## 2. **Running your problem**: Visualization, single-point calculations, and full iterative optimization
+
+Having implemented your script, you can execute it as a python script to
+run various calculations specified by command-line options.
+
+### 2a. Visualizing your geometry
+
+### 2b. Evaluating the objective function value and gradient at a single design point
+
+### 2c. Testing gradient components by finite-differencing
+
+### 2d. Running many single-point calculations in parallel: `ParallelDesignTester`
+
+### 2e. Full iterative optimization
+
+### 2f. Running many single-point calculations in parallel: `ParallelDesignTester`
+
+
+## 4. Built-in command-line options
+
+The following command-line options are defined by the `OptimizationProblem`
+base class and are available in all `mp.adjoint` sessions.
+
+### Options affecting Meep timestepping
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--res` | resolution
+| `--dpml` | PML thickness (-1 --> autodetermined)
+| `--fcen` | center frequency
+| `--df` | frequency width
+| `--source_mode` | mode index of eigenmode source
+| `--dft_reltol` | convergence threshold for end of timestepping
+| `--dft_timeout` | max runtime in units of last_source_time
+| `--dft_interval` | meep time DFT convergence checks in units of last_source_time
+
+### Options affecting outputs from Meep computations
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--nfreq` | number of output frequencies
+| `--full_dfts` | compute DFT fields over full volume
+| `--complex_fields` | force complex fields
+| `--filebase` | base name of output files
+
+### Options specifying initial values for basis-function coefficients
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--betafile` | file of expansion coefficients
+| `--beta` | set value of expansion coefficient
+| `--eps_design` | functional expression for initial design permittivity
+
+### Options describing the calculation to be done
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--eval_objective` | evaluate objective function value
+| `--eval_gradient` | evaluate objective function value and gradient
+| `--gradient_qname` | name of objective quantity to differentiate via adjoint method
+| `--fd_order` | finite-difference order (0,1,2)
+| `--fd_index` | index of differentiation variable
+| `--fd_rel_delta` | relative finite-difference delta
+| `--optimize` | perform automated design optimization
+
+### Options affecting optimization
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--alpha` | gradient descent relaxation parameter
+| `--min_alpha` | minimum value of alpha
+| `--max_alpha` | maximum value of alpha
+| `--boldness` | sometimes you just gotta live a little
+| `--timidity` | can\'t be too careful in this dangerous world
+| `--max_iters` | max number of optimization iterations
+
+### Options configurating adjoint-solver options
+
+| Option | Description |
+| --------------------------------- | -------------- |
+| `--verbose` | produce more output
+| `--concise` | produce less output
+| `--visualize` | produce visualization graphics
+| `--label_source_regions` | label source regions in visualization plots
+| `--logfile` | log file name
+| `--pickle_data` | save state to binary data file
+| `--animate_component` | plot time-domain field component
+| `--animate_interval` | meep time between animation frames
+
+--8<-- "doc/docs/AdjointSolver/AdjointLinks.md"
diff --git a/doc/docs/AdjointSolver/Visualization.md b/doc/docs/AdjointSolver/Visualization.md
new file mode 100644
index 000000000..e16302943
--- /dev/null
+++ b/doc/docs/AdjointSolver/Visualization.md
@@ -0,0 +1,317 @@
+--8<-- "doc/docs/AdjointSolver/AdjointDocumentationStyleHeader.md"
+
+---
+# `meep.adjoint.Visualization`:
+# Easy, standardized Meep visualization routines
+# that *just work*
+---
+
+
+
+
+The `meep.adjoint` python module implementing
+Meep adjoint solver includes
+an extensive set of tools for graphical visualization
+of Meep sessions via
+[`matplotlib`][MatPlotLib].
+Although these tools are packaged with and used by
+`meep.adjoint`, they are independent of the adjoint solver
+and can be used for convenient visualization of *any* python-driven
+Meep session---indeed, it is a specific
+goal of the module that the visualization routines can
+be invoked, with no arguments, on any Meep geometry
+and will do something reasonable and useful in every case.
+
+These and other objectives of the initiative are described in the
+first section below, after which we present examples to illustrate
+the three key types of plot produced by the module and how the
+default options can be customized.
+
+> :bookmark: **tl;dr documentation for `meep.adjoint.visualize`**
+>
+ 0. **Motivation:** Plotting is *much* too hard in Meep, especially
+ at the lower end where you just want something quick and easy
+ to let you sanity-check your geometry.
+ This module attempts to
+ provide some simple functions that should *just work* in most cases
+ to provide some graphical intuition for what you're working with.
+ The spirit of the module is to provide functions that can be called
+ *with no customization options whatsoever* and will do useful
+ relevant things by default, but which can also be customized
+ in cases where you *do* want to take the time to spruce up the
+ output.
+>
+ 1. **Visualizing geometries before timestepping.** Call `visualize_sim` to view your simulation prior to timestepping:
+>
+ ```py
+ from meep.adjoint import visualize_sim
+
+ sim=mp.simulation(...)
+
+ # insert code for adding DFT cells here
+
+ visualize_sim(sim)
+ ```
+>
+
+>
+ 2. **Visualizing frequency-domain fields after timestepping.** Call `visualize_sim` (same function, same no arguments!) after timestepping to view
+ frequency-domain (DFT) fields superposed atop the geometry;
+ the function automatically looks inside the `simulation` to determine what
+ frequency-domain fields are available to plot.
+>
+ ```py
+ sim.run(until=200)
+ visualize_sim(sim)
+ ```
+>
+
+>
+>
+ 3. **Real-time visualization of time-domain fields during timestepping.** Use `AFEClient` to define
+ a custom-made [step function](https://meep.readthedocs.io/en/latest/Scheme_User_Interface/#run-and-step-functions)
+ that intermittently updates a plot of time-domain fields in real time as your
+ meep timestepping run progresses.
+>
+ ```py
+ from meep.adjoint import AFEClient
+
+ step_funcs=[ AFEClient(sim, ['Ez'], interval=2) ]
+ sim.run(*step_funcs, until=200)
+ ```
+>
+ Here are movies of what this produces for the
+ [holey waveguide](Overview#the-holey-waveguide) geometry,
+ showing an incident eigenmode field propagating across a waveguide
+ with and without a hole in it.
+>
+ + Without hole: [`Ez2WithoutHole.h264`](images/Ez2WithoutHole.h264)
+>
+ + With hole: [`Ez2WithHole.h264`](images/Ez2WithHole.h264)
+
+__________________________________________________
+
+
+
+> :bookmark:{.summary} **`table of contents`**
+>
+> 1. **Motivation and Philsophy**
+>
+> a,b. *(a couple of longer-winded discussions culminating in) the*
+>
+> c. *The motivating desiderata for `meep.adjoint.Visualization`*
+>
+> 2. **Easy visualization of simulation geometries before timestepping: `visualize_sim`**
+>
+> 3. **Easy real-time animation of time-domain fields during
+> Meep computations:`AFEClient` and `AFEServer`**
+>
+> 4. **Easy visualization of frequency-domain fields---in isolation or
+> superposed atop your material geometry**
+>
+> 5. **Customization:***
+>
+
+
+
+## Motivation and philosophy
+
+
+>:bulb: Skip down to the [**`tl;dr` version**](tldr-version)
+ of the motivation and goals of this initiative,
+ or click to expand the tabs below for a
+ longer-winder account.
+
+
+??? note "**The error-checking power of visual sanity checks for text-defined geometries**"
+
+ The fact that every aspect of a Meep
+ calculation---from the [material geometry][GeometricObject],
+ to the [absorbing boundaries][PML],
+ to the placement and orientation
+ of the [sources][Source],
+ to the positions and sizes of the
+ [flux monitors][FluxSpectra]---is
+ defined by lines of text in script files offers a welcome
+ change of pace from the immersive GUI experience of commercial
+ CAD tools, as well as a natural platform for complex,
+ 0ustomized calculations that would be unwieldy or impossible
+ to launch from a point-and-click pull-down interface.
+ Nonetheless, as convenient as it is to specify geometries
+ non-graphically, it's *also* quite useful---arguably even
+ *essential*---to review a graphical representation
+ of the geometry you input as interpreted by
+ Meep, both to confirm that your input
+ was processed as you expected and to catch the sorts of
+ inadvertent errors---a flux monitor in the wrong place,
+ a source pointing the wrong direction---that are
+ all too easy to overlook in your python script,
+ but instantly obvious upon visual inspection.
+
+ In addition to the usefulness of visual sanity checks
+ on your input geometry, it can be helpful to look at
+ graphical visualizations of the electric and magnetic
+ fields computed by Meep---both
+ the time-domain and the frequency-domain fields---and
+ in particular to review how the field distributions
+ evolve in space in the presence of your material geometry.
+ This is true even though, in many cases, the ultimate
+ desired output of your calculation will be non-graphical,
+ sucha a number---a power flux, an eigenmode coefficient,
+ etc---whose computation, again, will be described by
+ lines of text you write in a script file. If the calculation
+ works correctly as you wrote it and you get
+ a reasonable numerical output on your first try,
+ well, bully for you! You have no need of graphical
+ visualization aids. For the rest of us---who can't
+ figure out why our flux is 10 orders of magnitude
+ too small until we look at the time-domain power flow and
+ realize we set things up backward, so that all the
+ excitation power flows promptly *out of*, not into,
+ the computational cell---visualization tools can be
+ a lifesaver.
+
+
+??? note "Wait, isn't this a solved problem?"
+
+ So we're agreed that visualization---of both inputs (geometries)
+ and outputs (fields)---is a good thing. But isn't this a
+ solved problem? After all, Meep offers
+ plenty of routines for retrieving as much raw data on
+ geometries and fields as any user could want---and, once one
+ has the raw data, it's just a matter of choosing from among the
+ infinitude of available tools for plotting and visualization.
+ Indeed, already within the Meep documentation
+ itself one can find many different types of visualization generated
+ by many different types of tool. What more is left to say?
+
+
+??? note "Wanted: Quick, easy, *standard* visualization paradigms that *just work* on any geometry"
+
+ Despite the arguably non-issue status of the situation, I would argue
+ that the current situation is suboptimal in at least two ways.
+
+ + **(a)** The absence of a single, canonical solution means that, as a Meep user,
+ every time you feel the urge to visualize something you have to spend some time and effort
+ figuring out how you are going to do it---and then going through the hassle of setting that up.
+ In my experience, this was especially true when it came to making movies of the evolution
+ time-domain fields---although the problem had been solved with documentation a couple of times
+ in the manual, still somehow it felt like a major chore to set this up every time
+ I wanted it.
+
+ + **(b)** On a different note, the absence of generally-accepted visualization protocols
+ means that everybody's visualizations look different. This is not necessarily tragic, and
+ we wouldn't want to enforce sterile conformity, but it might be nice if there were *some*
+ notion of "canonical Meep visualization format" to serve as a common
+ language.
+
+
+>:bulb: The motivating desiderata for `meep.adjoint.Visualization`
+
+ With all of that by way of backstory, here were the guidelines that motivated
+ the design of these routines.
+
+ 1. They should **fully address the three everyday situations that cry out most
+ urgently for visualization support**, namely
+
+ - *Simulation geometry before timestepping.* A static image of the simulation geometry at the begininng of the calculation,
+ before timestepping has begun. The important targets for visualization are the following items and their
+ relative positioning vis-a-vis one another: **(a)** the material geometry, **(b)** the absorbers,
+ **(c)** the location and nature of sources, **(d)** the position and size of any DFT cells in which
+ frequency-domain fields are to be computed.
+
+ - *Time-domain fields during timestepping.* An animation or movie showing the time-domain fields evolving
+ throughout the timestepping. For one thing, this tells you instantly when you have e.g. oriented your
+ source backwards or otherwise misconfigured something obvious---you sit and watch all of
+ your "input" power go flying right out of the cell to the left instead of proceeding rightward
+ into the center like you intended---whereas the *frequency*-domain plots produced even by
+ such totally misbegotten simulations can be pretty effective at obscuring even obvious
+ culprits.
+
+ Also, after all, the ability to observe essentially time-domain phenomena is one of the
+ key reasons to use an FDTD tool in the first place; many frequency-domain problems
+ can be solved much more efficiently by frequency-domain solvers than by something like
+ Meep, but they don't produce nice movies showing e.g.
+ reflection and diffraction of wave packets from interfaces. So it seems particularly
+ appropriate for a time-domain-solver like Meep to offer
+ built-in functionality for animating time-domain field evolutions.
+
+ - *Frequency-domain fields after timestepping.* A static image showing the amplitudes of frequency-domain fields
+ in regions of DFT cells, with the spatial distributiono of Poynting flux plotted for DFT flux cells.
+
+ >:bulb: **A central goal of this module is to provide reliable, easy-to-use routines for
+ each of the aforementioned purposes.**
+
+ 2. *They should offer a hierarchy of customizability options* for users at different
+ places on the willingness-to-suffer-for-their-art spectrum.
+ Of course, the goals outlined above---routines that are *easy to use* and *just work*---suggest
+ an implementation strategy in which the tools make all the artistic decisions for you; the easiest
+ calling convention to remember is the one with no parameters. However, one of the nice things
+ about the [`matplotlib`][MatPlotLib] backend is how very *extremely*
+ customizable it is, and it shouldn't be hard to pass at least some of those options
+ on to users of `meep.adjoint.Visualization.`
+ In view of these considerations, the goal of the routines is to offer
+ a three-tiered hierarchy of customizability levels.
+
+ a. **Zero-effort calling model for lazy forgetful slackers**SZ
+
+ For users who can't be bothered to memorize or look up library function
+ prototypes---or who just want quick and dirty visualization at the
+ prototyping stage with the option to clean up later---it should be
+ possible to call the visualization routines on just the `simulation`
+ entity alone *with no configuration arguments whatsoever,*
+ with the module making reasonable choices to produce helpful graphics
+ in a standardized format.
+
+ b. **Minimal-effort calling model for the slightly more motivated**
+
+ For users who fall basically into the previous bucket but are
+ motivated to change one or another option that otherwise
+ obstructs the visualization (like, label font size too large),
+ there should be an easy way to configure global settings once,
+ after which the user can proceed lazily to call all visualization
+ routines with no customization arguments.
+
+ c. **Full-customization calling model for diligent nitpickers:**
+
+ Finally, for users who *do* want to invest time and effort to
+ configure plotting options to their liking, the visualization
+ module should provide copious customization options to allow
+ as much as possible of the visualization to be customized.
+
+
+## 2. **Easy visualization of simulation geometries before timestepping: `visualize_sim`**
+
+
+As described above, the first visualization task we consider
+is that of double-checking the various geometric inputs
+supplied to Meep when creating a
+[`simulation`](TheSimulationClass):
+the computational cell and material geometry, the absorbing boundaries, the exciting sources,
+and any DFT cells for which we requested tabulation of frequency-domain fields. Also described
+above was the imperative that, to ensure a maximally memorable calling convention,
+the visualization routine be callable on *just the `simulation` class itself with no
+other parameters required.*
+
+These objectives are satisfied by the `visualize_sim` routine in `meep.adjoint.Visualization`,
+whose calling convention is precision-engineered to be as painless as possible: after
+creating a Meep `simulation`
+through a call like `sim=mp.simulation(geometry=...)`,
+you can simply say
+
+```py
+ visualize_sim(sim)
+```
+
+and voilà! A new plot window sprouts open to reveal your simulation
+in all its glory, looking something like this:
+
+
+
+## During timestepping: Real-time visualization of time-domain field evolution
+
+## After timestepping: Visualizing frequency-domain Poynting fluxes and energy densities
+
+## Customization
+
+--8<-- "doc/docs/AdjointSolver/AdjointLinks.md"
diff --git a/doc/docs/AdjointSolver/images/AdjointToyEz0.png b/doc/docs/AdjointSolver/images/AdjointToyEz0.png
new file mode 100644
index 000000000..51cf0cf7a
Binary files /dev/null and b/doc/docs/AdjointSolver/images/AdjointToyEz0.png differ
diff --git a/doc/docs/AdjointSolver/images/AdjointToyGeometry1.png b/doc/docs/AdjointSolver/images/AdjointToyGeometry1.png
new file mode 100644
index 000000000..7b740fb0b
Binary files /dev/null and b/doc/docs/AdjointSolver/images/AdjointToyGeometry1.png differ
diff --git a/doc/docs/AdjointSolver/images/AdjointToyGeometry2.png b/doc/docs/AdjointSolver/images/AdjointToyGeometry2.png
new file mode 100644
index 000000000..7fb91bb75
Binary files /dev/null and b/doc/docs/AdjointSolver/images/AdjointToyGeometry2.png differ
diff --git a/doc/docs/AdjointSolver/images/AdjointToydEz_FD.png b/doc/docs/AdjointSolver/images/AdjointToydEz_FD.png
new file mode 100644
index 000000000..beb9c7866
Binary files /dev/null and b/doc/docs/AdjointSolver/images/AdjointToydEz_FD.png differ
diff --git a/doc/docs/AdjointSolver/images/Ez2WithHole.h264 b/doc/docs/AdjointSolver/images/Ez2WithHole.h264
new file mode 100644
index 000000000..f1a5744f3
Binary files /dev/null and b/doc/docs/AdjointSolver/images/Ez2WithHole.h264 differ
diff --git a/doc/docs/AdjointSolver/images/Ez2WithoutHole.h264 b/doc/docs/AdjointSolver/images/Ez2WithoutHole.h264
new file mode 100644
index 000000000..3a20ace20
Binary files /dev/null and b/doc/docs/AdjointSolver/images/Ez2WithoutHole.h264 differ
diff --git a/doc/docs/AdjointSolver/images/HoleCloakBGeometry.png b/doc/docs/AdjointSolver/images/HoleCloakBGeometry.png
new file mode 100644
index 000000000..056cb887f
Binary files /dev/null and b/doc/docs/AdjointSolver/images/HoleCloakBGeometry.png differ
diff --git a/mkdocs.yml b/mkdocs.yml
index 99e14e724..4cba8c96b 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -78,6 +78,12 @@ pages:
- 'Tutorial/Mode Decomposition': 'Scheme_Tutorials/Mode_Decomposition.md'
- 'Tutorial/Casimir Forces': 'Scheme_Tutorials/Casimir_Forces.md'
- 'Guile and Scheme Information': 'Guile_and_Scheme_Information.md'
+- 'Adjoint Module for Automated Design Optimization':
+ - 'Overview': 'AdjointSolver/Overview.md'
+ - 'Reference Manual': 'AdjointSolver/ReferenceManual.md'
+ - 'Example Gallery': 'AdjointSolver/ExampleGallery.md'
+ - 'Visualization': 'AdjointSolver/Visualization.md'
+ - 'Implementation Notes': 'AdjointSolver/AdjointImplementationNotes.md'
- 'C++ Interface':
- 'Developer Information': 'C++_Developer_Information.md'
- 'Chunks and Symmetry': 'Chunks_and_Symmetry.md'