$$ \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} $$

Design optimization with the meep adjoint solver: A reference manual


As described in the Overview, 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.

🔖 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. Defining your objective function

    d. Objective regions, extra regions, design region

    e. Selecting an expansion basis

    f. Defining your own expansion basis

  2. *Running your problem, Part 1: *: Single-point calculations

    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

  3. *Running your problem, Part 2: *: Iterative design optimization

    a. Configuring and launching an iterative optimization

    b. Optimization algorithms

  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, 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 with two pure virtual methods that your derived class must override to describe the specifics of your design problem.

1b. Mandatory class-method overrides: init_problem and create_sim

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.)

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 communicate to the base class the CommonElementsOfOptimizationGeometries 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

    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.


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.


Optional class-method override: add_args

add_args: Configure problem-specific command-line arguments

The OptimizationProblem base class defines a large number of general-purpose 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

    A structure that has been initialized with all built-in mp.adjoint 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.


1d. Defining your objective function

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.

🔖 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.

1e. Selecting an expansion basis

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:

   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.

1f. Defining your own expansion basis

If none of the pre-implemented basis sets suffices for your needs, you may define your own arbitrary basis by implementing a new subclass of the Basis abstract base class in mp.adjoint. This class abstracts the representation of a spatially-varying permittivity function as a finite linear combination of scalar basis functions, $$ \epsilon(\mathbf{x})=\sum \beta_n b_n(\mathbf{x}).$$ The base class defines several methods related to this task and provides default non-performance-optimized implementations for most of these, which should be fine for most use cases. However, for design optimization of large structures there are two junctures, bookending each iteration of the iterative optimization process, in which a sluggish Basis implementation could potentially slow things down significantly: (1) At the beginning of each iteration, Basis.call() is invoked once for each Yee-grid pixel in the design region. (2) At the end of each iteration, Basis.project() is called to compute the basis-function projection of the objective-function gradient g(\mathbf{x}) \equiv \partial f/\partial \epsilon(\mathbf{x}).

The easiest way to implement your own basis is probably to start with one of the existing subclasses in mp.adjoint, such as SimpleFiniteElementBasis or PlaneWaveBasis, and modify appropriately, subject to the following guidelines.

  • Your subclass must furnish an override for the pure virtual __call__ method in the base class. This method accepts a single input argument p of type mp.Vector3 (coordinates of evaluation point) and returns the value of the full basis-function expansion at p, i.e. \sum \beta_n b_n(\mathbf{p}), with the $\beta_n$ coefficient values taken from thenp.array-valued class fieldbeta_vector` defined in the parent class.

  • The constructor (__init__ method) of your subclass must call the base class constructor, which accepts a single integer-valued parameter dim specifying the number of functions in the basis.

  • The base class provides default implementations of all other class methods, so in principle a bare-minimum subclass need not do anything beyond the above for for a bare-minimum implementation. However, as noted above, for performance reasons you may wish to override project() with an implementation that exploits special features of your basis to accelerate the calculation.

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