get_epsilon_grid function for evaluating ε on user-specified grid #1522
Conversation
|
Got it working. |
|
If you pass a list of normal geometries without a material grid, you get a segfault. gdb stops at Seems like there's a problem with the material not being assigned the proper type after it traverses the tree at Make sure the geometry list is properly saved and the right tree is stored when you call |
|
Thanks for looking into this. I had also used Lines 698 to 699 in 1d59e6b In For reference, the The strange thing is: why does this PR work with a
import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
n = 3.4 # index of waveguide
w = 1 # width of waveguide
r = 1 # inner radius of ring
pad = 4 # padding between waveguide and edge of PML
dpml = 2 # thickness of PML
sxy = 2*(r+w+pad+dpml) # cell size
def ring_resonator(p):
rr = (p.x**2+p.y**2)**0.5
if (rr > r) and (rr < r+w):
return mp.Medium(index=n)
return mp.air
geometry = [mp.Block(center=mp.Vector3(),
size=mp.Vector3(sxy,sxy),
material=ring_resonator)]
sim = mp.Simulation(cell_size=mp.Vector3(sxy,sxy),
geometry=geometry,
resolution=20)
cell_size = mp.Vector3(sxy,sxy,0)
x = np.linspace(-3,3,50)
y = np.linspace(-3,3,50)
z = np.array([0])
eps_grid = sim.get_epsilon_grid(x,y,z)
plt.figure()
plt.pcolor(x,y,eps_grid,cmap='Blues')
plt.gca().axis('equal')
plt.gca().set_aspect('equal','box')
plt.title("sampling material function using {} x {} grid".format(len(x),len(x)))
plt.savefig('ring_matfunc_eps_grid_nx{}_ny{}.png'.format(len(x),len(y)),dpi=120,bbox_inches='tight')For reference, the test which produces the segmentation fault is: import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
mp.verbosity(3)
n = 3.4 # index of waveguide
w = 1 # width of waveguide
r = 1 # inner radius of ring
pad = 4 # padding between waveguide and edge of PML
dpml = 2 # thickness of PML
sxy = 2*(r+w+pad+dpml) # cell size
c1 = mp.Cylinder(radius=r+w,
center=mp.Vector3(),
material=mp.Medium(index=n))
c2 = mp.Cylinder(radius=r,
center=mp.Vector3(),
material=mp.air)
sim = mp.Simulation(cell_size=mp.Vector3(sxy, sxy),
geometry=[c1, c2],
resolution=40)
x = np.linspace(-0.5*sxy,0.5*sxy,50)
y = np.linspace(-0.5*sxy,0.5*sxy,50)
z = np.array([0])
## triggers segmentation fault
eps_grid = sim.get_epsilon_grid(x,y,z) |
|
Output from |
|
Similar to |
|
Or just factor |
|
Once we debug the issues with the Lines 325 to 351 in 32788d7 This piece of code is called whenever there is a material grid present in a tree. The first if statement traverses through the tree and will interpolate as needed. The second if statement is a special case, which is called if the tree is done being traversed and the We need to account for the final case, where the current point doesn't have any material grids (e.g. it could be another medium or the it will require some debugging and trial and error. |
|
I tried refactoring The main parts of the Separately (but perhaps related), I noticed a (minor?) difference in the function
Line 1723 in d8ea59c However, Line 1951 in d8ea59c |
|
Actually, this seems to be working now. The fix involved only a single change: within the function init_libctl(_default_material, _ensure_periodicity, &gv,
cell_size, cell_center, &gobj_list);
geom_epsilon geps(gobj_list, mlist, gv.pad().surroundings());This is similar to the Line 1734 in d8ea59c The test case involving the
|
oskooi
changed the title
|
If I revert to passing the volume of only the user-specified subregion to the +std::complex<double> find_array_min_max(int n, const double *data) {
+ double min_val = data[0], max_val = data[0];
+ for (int i = 1; i < n; ++i) {
+ if (data[i] < min_val)
+ min_val = data[i];
+ if (data[i] > max_val)
+ max_val = data[i];
+ }
+ return std::complex<double>(min_val,max_val);
+}
void get_epsilon_grid(geometric_object_list gobj_list,
material_type_list mlist,
material_type _default_material,
@@ -2574,13 +2585,29 @@ void get_epsilon_grid(geometric_object_list gobj_list,
meep::grid_volume gv,
vector3 cell_size,
vector3 cell_center,
- int nx, double *x,
- int ny, double *y,
- int nz, double *z,
+ int nx, const double *x,
+ int ny, const double *y,
+ int nz, const double *z,
double *grid_vals) {
+ std::complex<double> min_max;
+ double min_val[3], max_val[3];
+ for (int n = 0; n < 3; ++n) {
+ if (((n == 0) ? nx : ((n == 1) ? ny : nz)) > 1) {
+ min_max = find_array_min_max(((n == 0) ? nx : ((n == 1) ? ny : nz)),
+ ((n == 0) ? x : ((n == 1) ? y : z)));
+ min_val[n] = real(min_max); max_val[n] = imag(min_max);
+ }
+ }
+ const meep::volume vol(meep::vec(min_val[0],min_val[1],min_val[2]),
+ meep::vec(max_val[0],max_val[1],max_val[2]));
+ master_printf("vol: (%0.3f,%0.3f,%0.3f), (%0.3f,%0.3f,%0.3f)\n",vol.in_direction_min(meep::X),vol.in_direction_min(meep::Y),vol.in_direction_min(meep::Z),vol.in_direction_max(meep::X),vol.in_direction_max(meep::Y),vol.in_direction_max(meep::Z));
+ meep::volume vol2 = gv.pad().surroundings();
+ master_printf("vol2: (%0.3f,%0.3f,%0.3f), (%0.3f,%0.3f,%0.3f)\n",vol2.in_direction_min(meep::X),vol2.in_direction_min(meep::Y),vol2.in_direction_min(meep::Z),vol2.in_direction_max(meep::X),vol2.in_direction_max(meep::Y),vol2.in_direction_max(meep::Z));
init_libctl(_default_material, _ensure_periodicity, &gv,
cell_size, cell_center, &gobj_list);
- geom_epsilon geps(gobj_list, mlist, gv.pad().surroundings());
+ dim = gv.dim;
+ // geom_epsilon geps(gobj_list, mlist, gv.pad().surroundings());
+ geom_epsilon geps(gobj_list, mlist, vol); |
|
I have discovered some strange behavior which may provide some clues for debugging. For a test case involving two overlapping This behavior can be demonstrated by running the same script 20 separate times and printing out how many values in the array returned by A clue to what could be causing this behavior is shown in the details of the geometry tree produced by the verbose output from three separate runs for the two test cases of incorrect and correct results: incorrect 1 incorrect 2 incorrect 3 correct 1 correct 2 correct 3 The difference in the geometry trees is highlighted by the lines with These results were generated using this test script: import meep as mp
import numpy as np
mp.verbosity(3)
n = 3.4 # index of waveguide
w = 1 # width of waveguide
r = 1 # inner radius of ring
pad = 4 # padding between waveguide and edge of PML
dpml = 2 # thickness of PML
sxy = 2*(r+w+pad+dpml) # cell size
c1 = mp.Cylinder(radius=r+w,
center=mp.Vector3(),
material=mp.Medium(index=n),
height=mp.inf)
c2 = mp.Cylinder(radius=r,
center=mp.Vector3(),
material=mp.air,
height=mp.inf)
sim = mp.Simulation(cell_size=mp.Vector3(sxy, sxy),
geometry=[c1,c2],
resolution=20)
x = np.linspace(-2.0,2.0,5)
y = np.linspace(-2.0,2.0,5)
z = np.array([0])
eps_grid = sim.get_epsilon_grid(x,y,z)
print("num. grid points with epsilon > 1.0: {}".format(np.where(eps_grid > 1.0)[0].shape)) |
|
Turns out the random results were due to uninitialized values of the user-specified volume in the z coordinate direction (as part of the The latest commit fixes this bug by initializing the coordinates of the empty dimensions of the user-specified volume to |
…noComp#1522) * get_epsilon_grid function for evaluating ε on user-specified grid * refactor fragment_stats::init_libctl out of class fragment_stats * clean up * initialize coordinates of user-specified volume with empty dimensions to zero * remove if statment from get_epsilon_grid * modify find_array_min_max to accept arguments by reference rather than return complex double * fixes
Closes #1509. Follows the outline suggested by @stevengj in #1509 (comment). Also includes documentation.
This feature is mostly working as demonstrated in the results shown below involving a cylinder as a
MaterialGrid(with the projection step turned off viabeta = 0) at a resolution of50which is sampled at three different resolutions (20,50, and100).Ideally,
Simulation.init_sim()should not be necessary (as demonstrated in the example below) sinceget_material_ptjust requires ageom_epsilonobject which itself does not involve initializing the Yee grid viastructure::set_materials. However, whenever thegeometrycontainsPrismobjects (including those loaded from a GDS file), it turns out thatSimulation.init_sim()is necessary otherwise there is a segmentation fault due to uninitialized materials. Also, this PR does not seem to work when thegeometryconsists ofGeometricObjectssuch asCylinders, etc. It doesn't seem to be obvious why this is the case. This PR should therefore not be merged until these two issues are resolved.