Error when doing adjoint optimization with fabrication constraint

[satadrudas]

Hi, I have been playing around with the codes here on different approaches to adjoint optimization. Currently exploring the fabrication constraint part as given in the tutorial.

There is a part of the code as follows where a condition is given that , if beta>=beta_thresh, beta=0 :

        # execute a single forward run before the start of each
        # epoch and manually set the initial epigraph variable to
        # slightly larger than the largest value of the objective
        # function over the six wavelengths and the lengthscale
        # constraint (final epoch only).
        t0 = opt(
            [
                mapping(
                    x[1:],
                    eta_i,
                    beta if beta < beta_thresh else 0,
                ),
            ],
            need_gradient=False,
        )

Could someone please explain why is the beta considered zero there? Shouldn't it be some large value to binarize it?. Besides beta=0 gives the following warning:

The weights parameter of MaterialGrid must be in the range [0,1]

also another error that I get which I think I am getting because of beta=0:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[9], line 74
     59     solver.add_inequality_mconstraint(
     60         lambda rr, xx, gg: glc(
     61             rr,
   (...)
     66         tol_lw,
     67     )
     69 # execute a single forward run before the start of each
     70 # epoch and manually set the initial epigraph variable to
     71 # slightly larger than the largest value of the objective
     72 # function over the six wavelengths and the lengthscale
     73 # constraint (final epoch only).
---> 74 t0,_ = opt(
     75     [
     76         mapping(
     77             x[1:],
     78             eta_i,
     79             beta if beta < beta_thresh else 0,
     80         ),
     81     ],
     82     need_gradient=False,
     83 )
     85 x[0] = np.amax(t0)
...
File [~/anaconda3/envs/meep/lib/python3.11/site-packages/meep/__init__.py:3040](https://file+.vscode-resource.vscode-cdn.net/Users/satadrudas/invdes/Nanophotonic_Design_Optimization_Using_Ising_Machine/MEEP%20codes/~/anaconda3/envs/meep/lib/python3.11/site-packages/meep/__init__.py:3040), in fields.step(self)
   3039 def step(self):
-> 3040     return _meep.fields_step(self)

RuntimeError: meep: simulation fields are NaN or Inf

The main body of the code:

if __name__ == "__main__":
    input_flux, input_flux_data = straight_waveguide()

    algorithm = nlopt.LD_MMA

    # number of design parameters
    n = Nx * Ny

    # initial design parameters
    x = np.ones((n,)) * 0.5
    x[Si_mask.flatten()] = 1.0  # set the edges of waveguides to silicon
    x[SiO2_mask.flatten()] = 0.0  # set the other edges to SiO2

    # lower and upper bounds for design weights
    lb = np.zeros((n,))
    lb[Si_mask.flatten()] = 1.0
    ub = np.ones((n,))
    ub[SiO2_mask.flatten()] = 0.0

    # insert epigraph variable initial value (arbitrary) and bounds into the
    # design array. the actual value is determined by the objective and
    # constraint functions below.
    x = np.insert(x, 0, 1.2)
    lb = np.insert(lb, 0, -np.inf)
    ub = np.insert(ub, 0, +np.inf)

    objfunc_history = []
    epivar_history = []
    cur_iter = [0]

    beta_thresh = 64 # threshold beta above which to use subpixel smoothing
    betas = [8, 16, 32, 64, 128, 256]
    max_evals = [80, 80, 100, 120, 120, 100]
    tol_epi = np.array([1e-4] * 2 * len(frqs))  # R, 1-T
    tol_lw = np.array([1e-8] * 2)  # line width, line spacing

    for beta, max_eval in zip(betas, max_evals):
        solver = nlopt.opt(algorithm, n + 1)
        solver.set_lower_bounds(lb)
        solver.set_upper_bounds(ub)
        solver.set_min_objective(f)
        solver.set_maxeval(max_eval)
        solver.set_param("dual_ftol_rel", 1e-7)
        solver.add_inequality_mconstraint(
            lambda rr, xx, gg: c(
                rr,
                xx,
                gg,
                eta_i,
                beta,
                False if beta < beta_thresh else True,
            ),
            tol_epi,
        )
        solver.set_param("verbosity", 1)

        opt = mode_converter_optimization(
            input_flux,
            input_flux_data,
            True,  # use_damping
            False if beta < beta_thresh else True,  # use_epsavg
            beta,
        )

        # apply the minimum linewidth constraint
        # only in the final epoch to an initial
        # binary design from the previous epoch.
        if beta == betas[-1]:
            res = np.zeros(2)
            grd = np.zeros((2, n + 1))
            t = glc(res, x, grd, beta)
            solver.add_inequality_mconstraint(
                lambda rr, xx, gg: glc(
                    rr,
                    xx,
                    gg,
                    beta,
                ),
                tol_lw,
            )

        # execute a single forward run before the start of each
        # epoch and manually set the initial epigraph variable to
        # slightly larger than the largest value of the objective
        # function over the six wavelengths and the lengthscale
        # constraint (final epoch only).
        t0 = opt(
            [
                mapping(
                    x[1:],
                    eta_i,
                    beta if beta < beta_thresh else 0,
                ),
            ],
            need_gradient=False,
        )
        t0 = np.concatenate((t0[0][0], t0[0][1]))
        t0_str = '[' + ','.join(str(tt) for tt in t0) + ']'
        x[0] = np.amax(t0)
        x[0] = 1.05 * (max(x[0], t) if beta == betas[-1] else x[0])
        print(f"data:, {beta}, {t0_str}, {x[0]}")

        x[:] = solver.optimize(x)

[oskooi]

There is a part of the code as follows where a condition is given that , if beta>=beta_thresh, beta=0. Could someone please explain why is the beta considered zero there? Shouldn't it be some large value to binarize it?

The reason for setting beta = 0 in the MaterialGrid definition for large beta is described in "important item" six of the tutorial:

The subpixel-smoothing feature of the `MaterialGrid` is necessary whenever the β parameter of the projection 
function is large (> ~50) and thus the design is binary (or nearly so). Without subpixel smoothing, then the 
gradients are nearly zero except for u≈η±1/β where the derivatives and second derivatives blow up, causing 
optimization algorithms to break down. When subpixel smoothing is enabled (do_averaging=True), the 
weights are projected internally using the beta parameter. For this reason, any preprocessing (i.e., mapping) 
of the weights outside of the MaterialGrid should apply only a filter to the weights but must not perform any 
projection.

An upcoming paper will describe the details of the subpixel-smoothing feature for the MaterialGrid.