enable support for subpixel smoothing of MaterialGrid #1503
Conversation
|
Here is a plot of the relative error vs. resolution for the
In terms of runtime performance, the larger the resolution of the |
Unfortunately, this isn't true. Remember, the material grid stores continuous values from 0 to 1, up to whatever precision you specify. In theory, you could use a binary array... but this not only kills your optimization, but the SWIG code still casts this as a 64-bit float array. More importantly, the whole point of TO is to leverage continuous optimization methods on a large scale. That being said, the storage requirements aren't terrible from meep's perspective. Most of the other data arrays are far bigger. Rather, the optimizer needs to be able to handle an array of that size. Normally MMA algorithms (e.g. CCSA in NLOPT) should work just fine. However, if you are using MPI, then this algorithm runs on each process, which means this array is copied to each process. You could of course get around this by only optimizing on one processor and broadcasting the results, but you just need to make sure your bookkeeping is correct. |
|
It occurs to me that we may be able to do better if we implement the projection step as part of the MaterialGrid rather than in Python. The basic reason for this is that, before projection, we have a nice smooth (low-pass-filtered) function, for which bilinear interpolation is a good description. That is, suppose g∈[0,1]ⁿ is our MaterialGrid data (before projection) and L(g,𝐱) ∈ [0,1] denotes the bilinear interpolation of g at 𝐱. Then, we define ε at each point in space (at "infinite" resolution) by ε(𝐱) = P(ε₁,ε₂,L(g,𝐱)) where P is our projection function onto [ε₁,ε₂]. Our mean ⟨ε⟩ etcetera is then defined by integrals of ε(𝐱), which we can sample as finely in space as we want (e.g. at 30× the resolution, or using specialized quadrature schemes). (Another way of putting it is that we should interpolate and then project, rather than project on the grid and then interpolate.) |
|
In particular, the MaterialGrid object would have two additional parameters β and η (for the tanh projection, taken from Python). Then |
|
For reference, the hyperbolic tangent function that is to be replicated in C is: meep/python/adjoint/filters.py Line 583 in 7067536 We will need to use the |
|
The point is that this then decouples the MaterialGrid resolution from the Meep resolution. The MaterialGrid effectively defines an discontinuous ε(x) at infinite resolution, which you can then either do subpixel averaging on (to get the Meep ε) or sample at any desired resolution to ship to the foundry (by upsampling & then projecting), so that the Meep calculation is consistent (to the degree allowed by the Meep resolution) with the foundry design. (Conversely, if you project and then interpolate, then the MaterialGrid defines a continuous ε(x) at infinite resolution, which is not what you want.) Note also that you would test convergence by keeping the MaterialGrid resolution fixed and increasing the Meep resolution. |
|
The MaterialGrid should have a high enough resolution that the "corners" induced in the interpolated+projected process (at infinite resolution you will have "corners" at the boundaries between MaterialGrid pixels) are shallow enough that they don't affect convergence much. For a test case, I would probably use a smooth MaterialGrid input, similar to what we plan to use in the adjoint solver with filters. You can either get this by taking your cylinder and filtering it. or just by putting in u = Gaussian "bump" or similar. Use a large β (say 1000) to get a binary structure, and η=0.5. |
|
This seems to be working and producing quadratic convergence using a fixed resolution for the The filtered design which is then interpolated is shown in the figure inset. Note that the boundaries have been blurred using SciPy's Two new parameters
import numpy as np
from scipy.ndimage import gaussian_filter
import argparse
import meep as mp
import matplotlib
import time
parser = argparse.ArgumentParser()
parser.add_argument('-res',
type=int,
default=20,
help='resolution (default: 20 pixels/um)')
parser.add_argument('-geom_type',
type=int,
choices=[1,2],
default=1,
help='type of geometry: 1: Cylinder (default), 2: material grid')
args = parser.parse_args()
cell_size = mp.Vector3(1,1,0)
rad = 0.301943
if args.geom_type == 1:
geometry = [mp.Cylinder(radius=rad,
center=mp.Vector3(),
height=mp.inf,
material=mp.Medium(index=3.5))]
else:
design_shape = mp.Vector3(1,1,0)
design_region_resolution = 50
Nx = int(design_region_resolution*design_shape.x)
Ny = int(design_region_resolution*design_shape.y)
x = np.linspace(-0.5*cell_size.x,0.5*cell_size.x,Nx)
y = np.linspace(-0.5*cell_size.y,0.5*cell_size.y,Ny)
xv, yv = np.meshgrid(x,y)
design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad
filtered_design_params = gaussian_filter(design_params, sigma=3.0, mode='nearest', output=np.double)
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
mp.Medium(index=3.5),
design_parameters=filtered_design_params,
grid_type='U_SUM',
do_averaging=True,
beta=1000,
eta=0.5)
geometry = [mp.Block(center=mp.Vector3(),
size=mp.Vector3(design_shape.x,design_shape.y,0),
material=matgrid)]
fcen = 0.3
df = 0.2*fcen
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=mp.Hz,
center=mp.Vector3(-0.1057,0.2094,0))]
k_point = mp.Vector3(0.3892,0.1597,0)
sim = mp.Simulation(resolution=args.res,
cell_size=cell_size,
geometry=geometry,
sources=sources,
k_point=k_point,
subpixel_maxeval=100000,
subpixel_tol=1e-5)
start = time.time()
sim.init_sim()
print("set_epsilon:, {} s".format(time.time()-start))
h = mp.Harminv(mp.Hz, mp.Vector3(0.3718,-0.2076), fcen, df)
sim.run(mp.after_sources(h),
until_after_sources=200)
for m in h.modes:
print("harminv:, {}, {}, {}".format(args.res,m.freq,m.Q)) |
|
Looks good! It might be nice to have...
And of course, we still need to add gradient support. |
|
This is an interesting approach to getting subpixel smoothing working with I'm just wondering if it would be useful to keep a more general interface for |
I'm not quite sure which three fields you're referring to? The proposed change would tightly couple the material grid to a particular projection operation, but even then you could turn this projection off (set β=0) and do the projection separately in any way that you want (albeit without subpixel smoothing). |
If Note that you can still parameterize however you want before you pass into the material grid. You can do as many filters, thresholds, or function fits as you have time for. |
|
It also occurs to me that we can now compute the normal vector directly via by averaging the gradient ∇u of the interpolated field before projection, similar to the discussion #1229 (since it's essentially a level-set function), rather than by spherical quadrature. (We could also do the ε averaging in a more clever way too.) |
I was thinking of different kinds and combinations of projection / filtering operations, along the lines of what @smartalecH mentioned. It initially struck me as less flexible to hardcode one of these operations into the Meep
Is there an approach that would allow for subpixel smoothing in the It seems like we just need a way to provide Meep with a continuous (or at least effectively continuous) representation of a design here? Above, my very high resolution binarized design might be one example of an effectively continuous design. The smooth level set variable (after interpolation) might be the continuous representation in the three field case. Is there an approach that could support both of these cases? |
Yes. Just use a MaterialGrid with β=0 at a high resolution, which still works with this PR (at a cost in memory). |
Why do you need a helper function here, instead of just directly accessing the fields, i.e. |
Say, for example you already have a design that works pretty well, but maybe you just want to eek out a bit of performance with another round of optimization. As @stevengj mentioned, you could still use sub-pixel averaging on this, provided the original design has a really high resolution and you turn off the threshold. Or you could do a mini topology optimization problem, where you want to find the new set of design parameters that, after passing through a filter and projection, correspond to your initial starting guess. Because of the nonlinear projection, it's not going to be a convex problem, but should be pretty straightfoward to solve (I've done this quite a few times with meep's adjoint framework). In fact, you can use the existing autograd libraries to efficiently compute the gradients. Then, with this new set of design parameters, you can continue to optimize with the projection step and take advantage of the subpixel averaging without having to hypersample your image (lots of memory) because you are passing it through a filter first. This methodology can be applied to other use cases too. |
Yeah, performing re-optimization is certainly an interesting use case there. As I mentioned above, I think what you're describing is one way of creating an effectively continuous representation in the three field parameterization: the bilinear interpolation of a relatively coarse filtered density array will give Meep a smooth level set-like distribution that can be thresholded. That's all good, but what strikes me as odd is the choice of hardcoding this into the Meep interface, when you could instead have a Meep adjoint module that is much more general. Not all users may be optimizing with a three field parameterization. A user could have some other function for generating their geometry (either that they want to optimize or that they just want to simulate with the benefit of subpixel smoothing). Being required to pass this geometry to Meep as a very high resolution numpy array seems less flexible and may not be the most convenient. For example, perhaps the geometry isn't coming from an array in the first place. You could still support the three field use case you describe above but with an interface that is more extensible to enable other approaches. I think what I am imagining here actually looks very similar to the material function interface that already exists? You give Meep a way of querying density values at arbitrary (finely) sampled points. Maybe you need to assume some well-behavedness or smoothness of the function, but that seems fair. You would also need an interface to pass the sensitivity / gradient at those points back out, but you basically already have this. You could then build the three field-specific scheme with the bilinear interpolation on top of this API, but still allow flexibility for other use cases. |
You already have a way of querying density values at arbitrarily finely sampled points with the MaterialGrid (which currently offers linear interpolation from the grid). The basic issue here is precisely that this is a poor representation of discontinuous functions unless you have a very fine grid resolution, or unless the MaterialGrid is treated as level-set function (i.e. an interface at u(x)=η). Incorporating the (optional) projection into the MaterialGrid precisely corresponds to allowing you to treat it as a level set (with optional smoothing β and offset η), which is a very general way to describe discontinuous materials that still allows efficient subpixel smoothing and other operations. It is not specific to a "three-field" parameterization that you keep referring to.
We already have this feature, and you can already do subpixel smoothing with it. The basic problem is that it is extremely slow because you have to call back to Python many times per pixel. If there are other flexible interfaces we could implement to specify the geometry, we can certainly consider them. But I think level sets are a useful general feature to have. |
Thanks for the clarification. I incorrectly thought that subpixel smoothing was not supported in that interface. I suppose there's also the issue of getting gradients if one wanted to do optimization using that but, more importantly, I can imagine that it is slow. Just a thought... but could an alternative interface be one where Meep provides the coordinates of the points where it wants to perform the evaluations? This might allow one to take advantage of vectorization and other speed ups on the Python side, depending on how the function is computed. |
It's pretty hard to do subpixel smoothing that way, because the subpixel quadrature points required would be hard to compute in advance (e.g. with vectorized code) without massive memory usage. |
* enable support for subpixel smoothing of MaterialGrid * apply hyperbolic tangent projection/thresholding after bilinear interpolation * fixes * minor fix to if statement
Closes #1500.
After some recent consultation with @smartalecH, it turns out enabling subpixel smoothing for the
MaterialGridis actually much easier than expected. All that is required is to make sure that the parameterfallbackis set totrueinsidegeom_epsilon::eff_chi1inv_matrixin order to use the existing adaptive quadrature routine from libctl:meep/src/meepgeom.cpp
Line 874 in 7067536
Because the routine
geom_epsilon::fallback_chi1inv_rowinvokesgeom_epsilon::get_material_pt(material_type &material, const meep::vec &r)to obtain ε at a given pointrfor which the use case ofmaterialas a material grid is already supported, no additional functions need to be added. In fact, setting it up this way means that subpixel smoothing also automatically supports symmetries for the material grid.This PR simply adds a new
do_averagingboolean property to theMaterialGridobject (default isFalse) which is passed from the Python interface intogeom_epsilon::eff_chi1inv_matrixin C/C++. This means that the existing parameterssubpixel_maxevalandsubpixel_tolfor theSimulationconstructor which are used to specify the properties of the quadrature scheme can be used as-is.The convergence test described in #1500 (shown below) verifies that the quadrature scheme for the
MaterialGridis second-order accurate (green line). This requires making sure that the resolution of theMaterialGridis suffficiently higher than Meep's Yee grid (in this example, the ratio is100). If the resolution of theMaterialGridis small relative to the Meep resolution (e.g., 2X), the convergence is closer to first-order accurate.The next feature to add is the gradient of the subpixel smoothing for back propagation.