From 35ec6d8ad9372eff7f5831625c19fff7cbe4ab6e Mon Sep 17 00:00:00 2001 From: smartalecH Date: Sat, 28 May 2022 20:24:26 -0400 Subject: [PATCH] update tests and add comments to code --- python/adjoint/filters.py | 18 ++++ python/requirements.txt | 1 + python/tests/test_material_grid.py | 80 +++++++++++----- src/meepgeom.cpp | 144 ++++++++++++----------------- 4 files changed, 139 insertions(+), 104 deletions(-) diff --git a/python/adjoint/filters.py b/python/adjoint/filters.py index c7eddd688..4e45dfb19 100644 --- a/python/adjoint/filters.py +++ b/python/adjoint/filters.py @@ -7,6 +7,8 @@ import meep as mp from scipy import special from scipy import signal +import skfmm +from scipy import interpolate def _proper_pad(x,n): ''' @@ -854,3 +856,19 @@ def gray_indicator(x): density-based topology optimization. Archive of Applied Mechanics, 86(1-2), 189-218. ''' return npa.mean(4 * x.flatten() * (1 - x.flatten())) * 100 + +def make_sdf(data): + ''' + Assume the input, data, is the desired output shape + (i.e. 1d, 2d, or 3d) and that it's values are between + 0 and 1. + ''' + # create signed distance function + sd = skfmm.distance(data- 0.5 , dx = 1) + + # interpolate zero-levelset onto 0.5-levelset + x = [np.min(sd.flatten()), 0, np.max(sd.flatten())] + y = [0, 0.5, 1] + f = interpolate.interp1d(x, y, kind='linear') + + return f(sd) diff --git a/python/requirements.txt b/python/requirements.txt index 772aac158..c43d34869 100644 --- a/python/requirements.txt +++ b/python/requirements.txt @@ -7,3 +7,4 @@ numpy parameterized pytest scipy +scikit-fmm diff --git a/python/tests/test_material_grid.py b/python/tests/test_material_grid.py index e8bee3562..1c4c73bbe 100644 --- a/python/tests/test_material_grid.py +++ b/python/tests/test_material_grid.py @@ -13,7 +13,6 @@ except: import adjoint as mpa import numpy as np -from scipy.ndimage import gaussian_filter import unittest rad = 0.301943 @@ -92,8 +91,7 @@ def compute_transmittance(matgrid_symmetry=False): return tran - -def compute_resonant_mode_2d(res,default_mat=False): +def compute_resonant_mode_2d(res,radius=rad,default_mat=False,cylinder=False,cylinder_matgrid=True,do_averaging=True): fcen = 0.3 df = 0.2*fcen sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), @@ -110,24 +108,37 @@ def compute_resonant_mode_2d(res,default_mat=False): x = np.linspace(-0.5*matgrid_size.x,0.5*matgrid_size.x,Nx) y = np.linspace(-0.5*matgrid_size.y,0.5*matgrid_size.y,Ny) xv, yv = np.meshgrid(x,y) - weights = np.sqrt(np.square(xv) + np.square(yv)) < rad - filtered_weights = gaussian_filter(weights, - sigma=3.0, - output=np.double) - + weights = mpa.make_sdf(np.sqrt(np.square(xv) + np.square(yv)) < radius) + if cylinder: + weights = 0.0*weights+1.0 + matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny), mp.air, Si, - weights=filtered_weights, - do_averaging=True, - beta=1000, + weights=weights, + beta=np.inf, eta=0.5) - - geometry = [mp.Block(center=mp.Vector3(), + + # Use a cylinder object, not a block object + if cylinder: + # within the cylinder, use a (uniform) material grid + if cylinder_matgrid: + geometry = [mp.Cylinder(center=mp.Vector3(), + radius=radius, + material=matgrid)] + # within the cylinder, just use a normal medium + else: + geometry = [mp.Cylinder(center=mp.Vector3(), + radius=radius, + material=Si)] + # use a block object + else: + geometry = [mp.Block(center=mp.Vector3(), size=mp.Vector3(matgrid_size.x,matgrid_size.y,0), material=matgrid)] sim = mp.Simulation(resolution=res, + eps_averaging=do_averaging, cell_size=cell_size, default_material=matgrid if default_mat else mp.Medium(), geometry=geometry if not default_mat else [], @@ -136,7 +147,7 @@ def compute_resonant_mode_2d(res,default_mat=False): h = mp.Harminv(mp.Hz, mp.Vector3(0.3718,-0.2076), fcen, df) sim.run(mp.after_sources(h), - until_after_sources=200) + until_after_sources=300) try: for m in h.modes: @@ -147,13 +158,12 @@ def compute_resonant_mode_2d(res,default_mat=False): return freq - -def compute_resonant_mode_mpb(resolution=1024): +def compute_resonant_mode_mpb(resolution=512): geometry = [mp.Cylinder(rad, material=Si)] geometry_lattice = mp.Lattice(cell_size) ms = mpb.ModeSolver(num_bands=1, - k_points=[k_point], + k_points=[mp.cartesian_to_reciprocal(k_point, geometry_lattice)], geometry=geometry, geometry_lattice=geometry_lattice, resolution=resolution) @@ -165,23 +175,51 @@ def compute_resonant_mode_mpb(resolution=1024): class TestMaterialGrid(unittest.TestCase): def test_subpixel_smoothing(self): + res = 25 + + def subpixel_test_matrix(radius,do_averaging): + freq_cylinder = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=True, cylinder_matgrid=False, do_averaging=do_averaging) + freq_matgrid = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=False, do_averaging=do_averaging) + freq_matgrid_cylinder = compute_resonant_mode_2d(res, radius, default_mat=False, cylinder=True, cylinder_matgrid=True, do_averaging=do_averaging) + return [freq_cylinder,freq_matgrid,freq_matgrid_cylinder] + + # when smoothing is off, all three tests should be identical to machine precision + no_smoothing = subpixel_test_matrix(rad,False) + self.assertEqual(no_smoothing[0], no_smoothing[1]) + self.assertEqual(no_smoothing[1], no_smoothing[2]) + + # when we slightly perturb the radius, the results should be the same as before. + no_smoothing_perturbed = subpixel_test_matrix(rad+0.01/res,False) + self.assertEqual(no_smoothing, no_smoothing_perturbed) + + # when smoothing is on, the simple material results should be different from the matgrid results + smoothing = subpixel_test_matrix(rad,True) + self.assertNotEqual(smoothing[0], smoothing[1]) + self.assertAlmostEqual(smoothing[1], smoothing[2], 3) + + # when we slighty perturb the radius, the results should all be different from before. + smoothing_perturbed = subpixel_test_matrix(rad+0.01/res,True) + self.assertNotEqual(smoothing, smoothing_perturbed) + # "exact" frequency computed using MPB freq_ref = compute_resonant_mode_mpb() - - res = [25, 50] + print(freq_ref) + res = [50, 100] freq_matgrid = [] for r in res: - freq_matgrid.append(compute_resonant_mode_2d(r)) + freq_matgrid.append(compute_resonant_mode_2d(r,cylinder=False)) # verify that the frequency of the resonant mode is # approximately equal to the reference value self.assertAlmostEqual(freq_ref, freq_matgrid[-1], 2) + print("results: ",freq_ref,freq_matgrid) # verify that the relative error is decreasing with increasing resolution # and is better than linear convergence because of subpixel smoothing self.assertLess(abs(freq_matgrid[1]-freq_ref)*(res[1]/res[0]), abs(freq_matgrid[0]-freq_ref)) - freq_matgrid_default_mat = compute_resonant_mode_2d(res[0], True) + # ensure that a material grid as a default material works + freq_matgrid_default_mat = compute_resonant_mode_2d(res[0], rad, True) self.assertAlmostEqual(freq_matgrid[0], freq_matgrid_default_mat) def test_symmetry(self): diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index eed88d58c..9f08db4d1 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -362,8 +362,14 @@ cvector3 to_geom_object_coords_VJP(cvector3 v, const geometric_object *o) { /* case geometric_object::CYLINDER: NOT YET IMPLEMENTED */ case geometric_object::BLOCK: { - vector3 v_real = cvector3_re(v); - vector3 v_imag = cvector3_im(v); + /* In order to leverage the underlying libctl infrastructure + *and* the dual library that makes computing derivatives so easy, + me perform a bit of a trick: we use the *complex* libctl vector, + storing the real and dual parts of the dual library and *manually* + propagate the dual through existing libraries as needed. + */ + vector3 v_real = cvector3_re(v); // real part + vector3 v_imag = cvector3_im(v); // "dual" part vector3 size = o->subclass.block_data->size; if (size.x != 0.0) {v_real.x /= size.x; v_imag.x /= size.x;} @@ -379,6 +385,9 @@ cvector3 to_geom_object_coords_VJP(cvector3 v, const geometric_object *o) { } cvector3 material_grid_grad(vector3 p, material_data *md, const geometric_object *o) { + /* computes the actual spatial gradient at point `p` + for the specified material grid `md`. */ + if (!is_material_grid(md)) {meep::abort("Invalid material grid detected.\n"); } cvector3 gradient = cvector_zero(); @@ -468,6 +477,14 @@ void map_lattice_coordinates(double &px, double &py, double &pz) { } cvector3 matgrid_grad(vector3 p, geom_box_tree tp, int oi, material_data *md) { + /* loops through all the material grids at a current point + and computes the final *spatial* gradient (w.r.t. x,y,z) + after all appropriate transformations (e.g. due to + overlapping grids). Calls the helper function, `material_grid_grad`, + which is what actually computes the spatial gradient. + */ + + // check for proper overlapping grids if (md->material_grid_kinds == material_data::U_MIN || md->material_grid_kinds == material_data::U_PROD) meep::abort("%s:%i:matgrid_grad does not support overlapping grids with U_MIN or U_PROD\n",__FILE__,__LINE__); @@ -493,6 +510,7 @@ cvector3 matgrid_grad(vector3 p, geom_box_tree tp, int oi, material_data *md) { ++matgrid_val_count; } + // compensate for overlapping grids if (md->material_grid_kinds == material_data::U_MEAN) gradient = cvector3_scale(1.0/matgrid_val_count,gradient); @@ -1250,12 +1268,11 @@ duals::duald get_material_grid_fill(meep::ndim dim, duals::duald d, double r, du return -1.0; // garbage fill } else { if (dim == meep::D1) - rel_fill = (r-d)/(2*r); - else if (dim == meep::D2 || dim == meep::Dcyl){ - rel_fill = (1/(r*r*meep::pi)) * (r*r*acos(d/r)-d*sqrt(r*r-d*d)); - } + return (r-d)/(2*r); + else if (dim == meep::D2 || dim == meep::Dcyl) + return (1/(r*r*meep::pi)) * (r*r*acos(d/r)-d*sqrt(r*r-d*d)); else if (dim == meep::D3) - rel_fill = (((r-d)*(r-d))/(4*meep::pi*r*r*r))*(2*r+d); + return (((r-d)*(r-d))/(4*meep::pi*r*r*r))*(2*r+d); } return rel_fill; @@ -1593,20 +1610,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3] 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)); - double uval = 0; - // TODO cleanup and remove - 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); - //uval = matgrid_val(p, tp, oi, md)+this->u_p; - } - else { - gradient = normal_vector(meep::type(c), v); - } - + meep::vec gradient = normal_vector(meep::type(c), v); get_material_pt(material, v.center()); material_epsmu(meep::type(c), material, &chi1p1, &chi1p1_inv); material_gc(material); @@ -1634,73 +1638,47 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3] integer errflag; double meps, minveps; - if (md->which_subclass == material_data::MATERIAL_GRID) { - number xmin[1], xmax[1]; - matgrid_volavg mgva; - mgva.dim = v.dim; - mgva.ugrad_abs = meep::abs(gradient); - mgva.uval = uval; - mgva.rad = v.diameter()/2; - mgva.beta = md->beta; - mgva.eta = md->eta; - mgva.eps1 = (md->medium_1.epsilon_diag.x+md->medium_1.epsilon_diag.y+md->medium_1.epsilon_diag.z)/3; - mgva.eps2 = (md->medium_2.epsilon_diag.x+md->medium_2.epsilon_diag.y+md->medium_2.epsilon_diag.z)/3; - xmin[0] = -v.diameter()/2; - xmax[0] = v.diameter()/2; -#ifdef CTL_HAS_COMPLEX_INTEGRATION - cnumber ret = cadaptive_integration(matgrid_ceps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval, - &esterr, &errflag); - meps = ret.re; - minveps = ret.im; -#else - meps = adaptive_integration(matgrid_eps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval, &esterr, - &errflag); - minveps = adaptive_integration(matgrid_inveps_func, xmin, xmax, 1, (void *)&mgva, 0, tol, maxeval, &esterr, - &errflag); -#endif + integer n; + number xmin[3], xmax[3]; + vector3 gvmin, gvmax; + gvmin = vec_to_vector3(v.get_min_corner()); + gvmax = vec_to_vector3(v.get_max_corner()); + xmin[0] = gvmin.x; + xmax[0] = gvmax.x; + if (dim == meep::Dcyl) { + xmin[1] = gvmin.z; + xmin[2] = gvmin.y; + xmax[1] = gvmax.z; + xmax[2] = gvmax.y; } else { - integer n; - number xmin[3], xmax[3]; - vector3 gvmin, gvmax; - gvmin = vec_to_vector3(v.get_min_corner()); - gvmax = vec_to_vector3(v.get_max_corner()); - xmin[0] = gvmin.x; - xmax[0] = gvmax.x; - if (dim == meep::Dcyl) { - xmin[1] = gvmin.z; - xmin[2] = gvmin.y; - xmax[1] = gvmax.z; - xmax[2] = gvmax.y; - } - else { - xmin[1] = gvmin.y; - xmin[2] = gvmin.z; - xmax[1] = gvmax.y; - xmax[2] = gvmax.z; - } - if (xmin[2] == xmax[2]) - n = xmin[1] == xmax[1] ? 1 : 2; - else - n = 3; - double vol = 1; - for (int i = 0; i < n; ++i) - vol *= xmax[i] - xmin[i]; - if (dim == meep::Dcyl) vol *= (xmin[0] + xmax[0]) * 0.5; - eps_ever_negative = 0; - func_ft = meep::type(c); + xmin[1] = gvmin.y; + xmin[2] = gvmin.z; + xmax[1] = gvmax.y; + xmax[2] = gvmax.z; + } + if (xmin[2] == xmax[2]) + n = xmin[1] == xmax[1] ? 1 : 2; + else + n = 3; + double vol = 1; + for (int i = 0; i < n; ++i) + vol *= xmax[i] - xmin[i]; + if (dim == meep::Dcyl) vol *= (xmin[0] + xmax[0]) * 0.5; + eps_ever_negative = 0; + func_ft = meep::type(c); #ifdef CTL_HAS_COMPLEX_INTEGRATION - cnumber ret = cadaptive_integration(ceps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, - &esterr, &errflag); - meps = ret.re / vol; - minveps = ret.im / vol; + cnumber ret = cadaptive_integration(ceps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, + &esterr, &errflag); + meps = ret.re / vol; + minveps = ret.im / vol; #else - meps = adaptive_integration(eps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr, - &errflag) / vol; - minveps = adaptive_integration(inveps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr, - &errflag) / vol; + meps = adaptive_integration(eps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr, + &errflag) / vol; + minveps = adaptive_integration(inveps_func, xmin, xmax, n, (void *)this, 0, tol, maxeval, &esterr, + &errflag) / vol; #endif - } + if (eps_ever_negative) // averaging negative eps causes instability minveps = 1.0 / (meps = eps(v.center())); {