compute material volume averages of a MaterialGrid using analytic formulation #1568
Conversation
|
The volume averages for ε and ε-1 are now computed by evaluating u just once at the voxel center. However, there is a large discrepancy in the results for the 2d photonic-crystal test case computed using the latest commit compared to the master branch:
The results obtained using this branch do not appear to be converging with second-order accuracy as would be expected. (Another problem is that these results could only be generated using the serial Meep because parallel Meep produced a segmentation fault in |
|
I would first test that you are doing the quadrature correctly. I would compare your 1d integral to (1) a 2d integral (over a ball) using the linear approximation of u(x,y) and (2) a 2d integral over a ball using the full material-grid u(x,y). (1) should agree exactly with your 1d integral, up to the quadrature tolerance, while (2) should agree approximately (becoming more and more accurate as the Meep resolution increases). Note that to an integral over a ball of radius R, you do a 2d integral over r=0..R and θ=0..2π, multiplying your integrand by a factor of r (as usual for cylindrical integrals). |
|
I fixed a bug related to how the |
|
I fixed another bug related to how the average ε and ε-1 were being computed and it seems to be working now. The results shown in the left figure below for a test problem involving a 2d photonic crystal is showing quadratic convergence in the limit of infinite resolution and matching the results from the master branch. This test problem involved using β=1000 for the projection of the
Currently, subpixel smoothing is not supported for overlapping grids. |
oskooi
changed the title
|
Support for overlapping |
|
It would be good to add the result for a geometric object to the same plot. (It would still be nice to compare the analytical integral against a brute-force integral over a sphere.) |
|
Adding the result for a
In this convergence plot, the resolution is successively doubled (10, 20, 40, 80, 160, 320) and the "exact" value in each case used to compute the relative error is at a resolution of 640. In the two plots above, the resolution is increased in fixed-size increments of 20 pixels/µm which is why the data is more noisy. |
|
Comparing the analytic integral (1d) to a brute-force integral of the material grid over a circular voxel (the second suggestion from @stevengj's comment above) produces results which seem to become more similar as the resolution is increased. This is demonstrated in the figure below which is a plot of the resonant frequency vs. resolution for the two cases.
For additional comparison of the two methods, here are the actual values for the average ε at 10 different interface voxels at a resolution of 20. The data consists of two columns: the first column is the xy coordinates of the voxel center and the second column is the average ε. analytic integral (1d) brute-force integral (circle, 2d) Note that the analytic integral produces a binary set of values (1.0 or 12.25) whereas the brute-force integral produces continuously varying values in the range [1.0,12.25]. Results for the analytic integral can be generated from the brute-force integral by "snapping" average ε values of 6.625 (corresponding to η=0.5) to ε=1 and those larger to ε=12.25. For additional reference, here is the convergence plot as well as the code for the brute-force integration from
struct matgrid_volavg {
vector3 med1_eps_diag; // diagonal of epsilon tensor from medium 1
vector3 med2_eps_diag; // diagonal of epsilon tensor from medium 2
vector3 cen; // voxel center
geom_box_tree tp;
int oi;
material_data *md;
};
#ifdef CTL_HAS_COMPLEX_INTEGRATION
static cnumber matgrid_ceps_func(int n, number *x, void *mgva_) {
matgrid_volavg *mgva = (matgrid_volavg *)mgva_;
vector3 p;
p.x = mgva->cen.x + x[0]*cos(x[1]);
p.y = mgva->cen.y + x[0]*sin(x[1]);
p.z = 0;
double uval = matgrid_val(p, mgva->tp, mgva->oi, mgva->md);
double u_proj = tanh_projection(uval, mgva->md->beta, mgva->md->eta);
vector3 med1_eps_diag = mgva->med1_eps_diag;
vector3 med2_eps_diag = mgva->med2_eps_diag;
double eps1 = (med1_eps_diag.x + med1_eps_diag.y + med1_eps_diag.z)/3;
double eps2 = (med2_eps_diag.x + med2_eps_diag.y + med2_eps_diag.z)/3;
cnumber ret;
ret.re = (1-u_proj)*eps1 + u_proj*eps2;
ret.im = (1-u_proj)/eps1 + u_proj/eps2;
return x[0]*ret;
}
#else
static number matgrid_eps_func(int n, number *x, void *mgva_) {
matgrid_volavg *mgva = (matgrid_volavg *)mgva_;
vector3 p;
p.x = mgva->cen.x + x[0]*cos(x[1]);
p.y = mgva->cen.y + x[0]*sin(x[1]);
p.z = 0;
double uval = matgrid_val(p, mgva->tp, mgva->oi, mgva->md);
double u_proj = tanh_projection(uval, mgva->md->beta, mgva->md->eta);
vector3 med1_eps_diag = mgva->med1_eps_diag;
vector3 med2_eps_diag = mgva->med2_eps_diag;
double eps1 = (med1_eps_diag.x + med1_eps_diag.y + med1_eps_diag.z)/3;
double eps2 = (med2_eps_diag.x + med2_eps_diag.y + med2_eps_diag.z)/3;
double eps_interp = (1-u_proj)*eps1 + u_proj*eps2;
return x[0]*eps_interp;
}
static number matgrid_inveps_func(int n, number *x, void *mgva_) {
matgrid_volavg *mgva = (matgrid_volavg *)mgva_;
vector3 p;
p.x = mgva->cen.x + x[0]*cos(x[1]);
p.y = mgva->cen.y + x[0]*sin(x[1]);
p.z = 0;
double uval = matgrid_val(p, mgva->tp, mgva->oi, mgva->md);
double u_proj = tanh_projection(uval, mgva->md->beta, mgva->md->eta);
vector3 med1_eps_diag = mgva->med1_eps_diag;
vector3 med2_eps_diag = mgva->med2_eps_diag;
double eps1 = (med1_eps_diag.x + med1_eps_diag.y + med1_eps_diag.z)/3;
double eps2 = (med2_eps_diag.x + med2_eps_diag.y + med2_eps_diag.z)/3;
double epsinv_interp = (1-u_proj)/eps1 + u_proj/eps2;
return x[0]*epsinv_interp;
}
#endif
// fallback meaneps using libctl's adaptive cubature routine
void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3],
const meep::volume &v, double tol, int maxeval) {
symmetric_matrix chi1p1, chi1p1_inv;
material_type material;
vector3 p = vec_to_vector3(v.center());
boolean inobject;
material =
(material_type)material_of_unshifted_point_in_tree_inobject(p, restricted_tree, &inobject);
material_data *md = material;
meep::vec gradient(zero_vec(v.dim));
matgrid_volavg mgva;
if (md->which_subclass == material_data::MATERIAL_GRID) {
geom_box_tree tp;
int oi;
tp = geom_tree_search(p, restricted_tree, &oi);
gradient = matgrid_grad(p, tp, oi, md);
mgva.tp = tp;
mgva.oi = oi;
mgva.md = md;
mgva.cen = p;
}
else {
gradient = normal_vector(meep::type(c), v);
}
if (md->which_subclass == material_data::MATERIAL_GRID) {
mgva.med1_eps_diag = md->medium_1.epsilon_diag;
mgva.med2_eps_diag = md->medium_2.epsilon_diag;
xmin[0] = 0;
xmin[1] = 0;
xmin[2] = 0;
xmax[0] = v.diameter()/2;
xmax[1] = 2*meep::pi;
xmax[2] = 0;
double vol = meep::pi * v.diameter()/2 * v.diameter()/2;
#ifdef CTL_HAS_COMPLEX_INTEGRATION
cnumber ret = cadaptive_integration(matgrid_ceps_func, xmin, xmax, 2, (void *)&mgva, 0, tol, maxeval,
&esterr, &errflag);
meps = ret.re/vol;
minveps = ret.im/vol;
#else
meps = adaptive_integration(matgrid_eps_func, xmin, xmax, 2, (void *)&mgva, 0, tol, maxeval, &esterr,
&errflag) / vol;
minveps = adaptive_integration(matgrid_inveps_func, xmin, xmax, 2, (void *)&mgva, 0, tol, maxeval, &esterr,
&errflag) / vol;
#endif
}
|
|
The mismatch between the brute-force and analytic integrands makes no sense to me — your filter radius is 60 pixels, so the u(x) function is very slowly varying on the scale of the pixel and a a 1st-order Taylor expansion should be nearly exact. I wouldn't bother with comparing convergence plots for now. (Almost anything you do, even no averaging at all, will converge eventually to the same result.) The key thing is to make sure you are computing the ε integral correctly for boundary pixels (ones where u passes through 0.5 within the pixel). The first step might be to compare your first-order taylor expansion to the actual bilinear-interpolated u(x) value for a few points in a boundary pixel. |
src/meepgeom.cpp
Outdated
| xmin[2] = 0; | ||
| xmax[0] = v.diameter()/2; | ||
| xmax[1] = 0; | ||
| xmax[2] = 0; |
There was a problem hiding this comment.
cadaptive_integration will only look at at the first component of xmin and xmax since your integral has dimension 1.
|
I think I know what is causing the mismatch between the brute-force and analytic integrands. The magnitude of the gradient ||∇u|| is actually off by a factor of 1/Δs2 where Δs = {Δx,Δy} is the voxel size of the
The missing 1/Δs2 factor is due to the mapping function used for the bilinear interpolation of the Lines 690 to 697 in cf14e30 In this mapping, Δs=1. The fix involves simply scaling the gradients by the correct 1/Δs2 values. |
…,dz using chain rule
|
After fixing a bug involving the scaling of the gradients, the average ε values at the interface voxels at a resolution of 1000 are nearly identical for the analytic 1d integral vs. brute force method: analytic integral (1d) brute-force integral (circle, 2d) relative error (as a percentage) |
|
For reference, here are the convergence plots comparing the analytic integral in 1d to the three other methods: (1) brute-force 2d integral (in polar coordinates) of the In the limit of infinite resolution, the results for the analytic integral (in blue) are converging with second-order accuracy (right figure) and are practically the same as methods (1) and (2) (left figure).
|
|
You also need to include the derivative of the As a test of this, use a square material grid (e.g. 100x100), but apply it to a rotated rectangular object. |
|
To include the derivative of the The function itself consists of the product of three things: (1) a 3x3 projection matrix Perhaps we will want to modify the function header for |
|
What you need is a function that computes the vector-Jacobian product of the gradient of the matgrid_val function with the Jacobian of the to_geom_box_coords function. This should really be a separate function, I think. Something like: vector3 to_geom_object_coords_VJP(vector3 v, geometric_object o) {
switch (o.which_subclass) {
default: {
vector3 po = {0, 0, 0};
return po;
}
case GEOM SPHERE: {
number radius = o.subclass.sphere_data->radius;
return vector3_scale(0.5 / radius, v);
}
/* case GEOM CYLINDER:
NOT YET IMPLEMENTED */
case GEOM BLOCK: {
vector3 size = o.subclass.block_data->size;
if (size.x != 0.0) v.x /= size.x;
if (size.y != 0.0) v.y /= size.y;
if (size.z != 0.0) v.z /= size.z;
return matrix3x3_transpose_vector3_mult(o.subclass.block_data->projection_matrix, v);
}
/* case GEOM PRISM:
NOT YET IMPLEMENTED */
}
}That is, let |
… of the geometric object using vector-Jacobian product
|
After adding the function test problem import numpy as np
from scipy.ndimage import gaussian_filter
import argparse
import meep as mp
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')
parser.add_argument('-design_region_res',
type=int,
default=50,
help='design region resolution (default: 50)')
parser.add_argument('-sigma',
type=float,
default=3.0,
help='standard deviation for Gaussian filter kernel')
parser.add_argument('-beta',
type=float,
default=1000.0,
help='turn-on rate of hyperbolic tangent (default: 1000.0)')
args = parser.parse_args()
cell_size = mp.Vector3(1.0,1.0,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 = args.design_region_res
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=args.sigma, output=np.double)
eta = 0.5
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
mp.Medium(index=3.5),
weights=filtered_design_params,
do_averaging=True,
beta=args.beta,
eta=eta)
geometry = [mp.Block(center=mp.Vector3(),
size=mp.Vector3(2*rad,3*rad,0),
e1=mp.Vector3(1,0,0).rotate(mp.Vector3(z=1),np.radians(30)),
e2=mp.Vector3(0,1,0).rotate(mp.Vector3(z=1),np.radians(30)),
material=matgrid)]
sim = mp.Simulation(resolution=args.res,
cell_size=cell_size,
geometry=geometry)
sim.init_sim()analytic integral (1d) brute-force integral (circle, 2d) relative error (as a percentage) edit
|
|
Some of those relative errors seem pretty large (e.g. 15%, 7%, etc.). Are the "bad" gradients localized to some region (or regions with similar features)? It might be good to visualize the relative error around the contours to identify "hot spots". It might uncover a bug. |
|
I'm not sure why there are those "outlier" voxels with relatively large errors. However, I verified that those maximum relative errors are reduced with increasing resolution of the Meep Yee grid and the |
From the small dataset above though, I'm not so sure those are outliers. There's 6/20 entries with errors >=5%. Rather, it seems indicative of a bug. In addition to my earlier suggestions (visualizing where these values are clustered) maybe try checking more samples. |
|
I should have mentioned that there are in fact over 5000 interface voxels at this resolution, of which I have shown only the first 20 entries. Note that in the results I posted earlier in this PR for a circle (i.e., not rotated and stretched), the range of relative errors is also large: [~0.001%, ~1.3%]. As I mentioned previously, the key point is that the maximum errors are reduced with increasing resolution and also the resonant frequencies computed using the 1d analytical integral and the brute force numerical integration over a circle (with diameter equal to the voxel width) are converging to the same result with increasing resolution. If there were indeed a bug somewhere, neither of these two results would be possible. edit: Could there be a problem not with this PR but with how the |
…mulation (NanoComp#1568) * compute material volume averages of a MaterialGrid using analytic formulation * fixes * evaluate weight of matgrid once at the voxel center * compute weights array at voxel center using geom_box coordinates * flip weights for two materials when computing averaged quantities * add support for overlapping grids * compute gradients with respect to x,y,z coordinates rather than dx,dy,dz using chain rule * fix bug in scaling of gradient in z direction * a few simplifications * compute gradient of the weights array with respect to the coordinates of the geometric object using vector-Jacobian product * check if geometric_object is not NULL in to_geom_object_coords_VJP * add chain rule for the map_lattice_coordinates function for default material * refactor to minimize copy-pasted lines * fix bug in matgrid_ceps_func * fixes and tweaks
|
The above tests describe the benefits of subpixel smoothing but don't really highlight the benefits of the hybrid levelset also introduced by this PR. As discussed in #1780, a simple example might involve perturbing a single pixel of a quasi FP cavity, and analyzing how the objective function (e.g. transmission) changes. I coded up the example here. The simulation looks like this: and the script perturbs a parameter on the left of the block. The results are below:
As expected, the new hybrid levelset method preserves a smoothly varying function of |
|
Why is the sharp transition in the results for the "traditional density method" around |
Good question. In short, it's because of the smoothing filter. In this particular example, the raw, unfiltered design parameters look something like Where But, as readily seen above, a value of |
Initial attempt to compute the volume averages for ε and ε-1 as part of subpixel smoothing using an analytic formulation based on a 1d line integral through a spherical voxel. This is intended to replace the existing approach involving an adaptive quadrature scheme over a cubic voxel.
Currently only supports 2d. It compiles but for a test case involving a 2d photonic crystal (same as the one used in #1539) the fields blow up.