filters used by adjoint solver sometimes produce slightly invalid results

[oskooi]

The filtering functions used by the adjoint solver sometimes produce invalid results. This can be demonstrated in the following example in which a 2d array of binary values (0 and 1) representing a ring resonator is passed to the conical filter (meep.adjoint.conic_filter). The filtered array which is returned has extrema which lie slightly outside the valid range of [0, 1]:

ring_weights:, min=-5.050993619243203e-16, max=1.0000000000000007

The violations do not go away with changes in the filter radius or array size. While this specific example involves minor violations (on the scale of machine precision) which are corrected by #1881 (which clips the weights input of the MaterialGrid to lie in the range [0, 1]), larger violations can potentially affect the performance of the optimizer if there happens to be a bug in the filters. It would be useful to ensure that the filters always return valid output given valid input. This could be a unit test.

import meep as mp
import meep.adjoint as mpa
import numpy as np

sxy = 5.0
cell_size = mp.Vector3(sxy,sxy,0)

design_region_size = mp.Vector3(1.5,1.5)
design_region_resolution = 128
Nx = int(design_region_resolution*design_region_size.x) + 1
Ny = int(design_region_resolution*design_region_size.y) + 1

filter_radius = 10/design_region_resolution

x = np.linspace(-0.5*design_region_size.x,0.5*design_region_size.x,Nx)
y = np.linspace(-0.5*design_region_size.y,0.5*design_region_size.y,Ny)
xv, yv = np.meshgrid(x,y)

rad = 0.538948295
wdt = 0.194838432
weights = np.where(np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                                  np.sqrt(np.square(xv) + np.square(yv)) < rad+wdt),
                   1.,
                   0.)

def mapping(x):
    filtered_weights = mpa.conic_filter(x,
                                        filter_radius,
                                        design_region_size.x,
                                        design_region_size.y,
                                        design_region_resolution)

    return filtered_weights.flatten()

ring_weights = mapping(weights)
print(f"ring_weights:, min={np.amin(ring_weights)}, max={np.amax(ring_weights)}")

[smartalecH]

This looks like roundoff error. You can always use the npa.clip() function to avoid this.