Add Python support for Prisms

[ChristopherHogan]

• Add prism Python class and conversion function for typemaps
• Move compare_arrays into a utils module so it's available for all tests.
• Replace blocks in python/tests/bend_flux.py with prisms. I had to increase the tolerance from 1e-7 to 1e-2 for the bend_flux.py test to pass. Is that alarming or expected?
  @stevengj @oskooi @HomerReid

[stevengj]

Why not just vertices[i] = v3?

[ChristopherHogan]

Done.

[stevengj]

Is that alarming or expected?

It means that it is discretized slightly differently than before. The real question is whether the discretization error has increased or is just different. A good way to test this might be:

1. re-run the bend test with blocks at as high a resolution as you can manage, and call this the "exact" result.

2. run the bend test at the resolutions 10 and 20 with both blocks and prisms, and compare the answer to the "exact" result — are the errors of similar magnitudes? Or are the prisms much worse?

If the prisms are much worse, then it would suggest a bug in the subpixel averaging of stevengj/libctl#13.

[ChristopherHogan]

It looks to me like the refl flux error is comparable for the Block and Prism, but the trans flux Prism error is a couple orders of magnitude larger than the Block. Here's the full data if someone wants to verify.
bend_flux.xlsx

[stevengj]

Did you look at output-epsilon to make sure that the structures look the same? e.g. maybe the prisms don't extend all the way to the edge of the cell like the blocks do?

[ChristopherHogan]

The epsilon files are slightly different for blocks and prisms at the edges of the geometry. Here is a vertical slice of the epsilon values for the horizontal waveguide made with a block

1.2778157
5.8317757
12       
12       
12       
12       
12       
12       
12       
12       
5.8317757
1.2778157

and here is the same thing for the prism

1.44      
12        
12        
12        
12        
12        
12        
12        
12        
12        
2.57142857

The block has a more uniform fade-in/fade-out. You can see the effect visually too (prism on the left, block on the right):

Are these supposed to be exactly equal? Is this an issue with subpixel averaging?

Additionally, the effect is different where the waveguide meets the edge of the cell. Here is the vertical waveguide where it meets the top of the cell boundary (block on left, prism on right):

When the block is the same size as the cell, it touches the edge, but when the prism is the cell size, it leaves a gap. This is a 16x32 cell, and the prism points at the top are (3, 16) and (4, 16).

Strangely, this doesn't happen with the horizontal waveguide. The points (-8, -11) and (-8, -12) actually touch the cell wall on the left (prism left, block right):

Notice, however, that the prism has the full epsilon value (12, all black) at the boundary, while the block has a blended region.

[ChristopherHogan]

Here's a minimal C++ reproducer representing an inconsistency I found while stepping through the block and prism versions line by line.

int main(int argc, char *argv[])
{
  meep_geom::material_type mat = meep_geom::make_dielectric(12);
  vector3 center = {0, 0, 0};
  vector3 e1 = {1, 0, 0};
  vector3 e2 = {0, 1, 0};
  vector3 e3 = {0, 0, 1};
  vector3 size = {1, 1, 0};
  geometric_object block = make_block(mat, center, e1, e2, e3, size);

  vector3 vertices[4] = {
    {-0.5, 0.5, 0},
    {0.5, 0.5, 0},
    {0.5, -0.5, 0},
    {-0.5, -0.5, 0}
  };
  vector3 axis = {0, 0, 1};
  geometric_object prism = make_prism(mat, vertices, 4, 0, axis);

  vector3 p = {-0.5, 0, 0};

  vector3 nhat_block = normal_to_object(p, block);
  vector3 nhat_prism = normal_to_object(p, prism);

  double tolerance = 1.0e-6;
  if (!vector3_nearly_equal(nhat_block, nhat_prism, tolerance)) {
    printf("FAIL: Normals not equal\n");
    printf("  Block: {%f, %f, %f}\n", nhat_block.x, nhat_block.y, nhat_block.z);
    printf("  Prism: {%f, %f, %f}\n", nhat_prism.x, nhat_prism.y, nhat_prism.z);
  }

  return 0;
}

[ChristopherHogan]

The issue is that normal_to_object(p, prism) returns a zero vector if p is on the prism boundary (e.g., (-0.5, 0) on a 1x1 square prism) which causes it to go to the noavg label in geom_epsilon::eff_chi1inv_matrix (meepgeom.cpp:863):

  normal = unit_vector3(normal_to_fixed_object(vector3_minus(p, shiftby), *o));
  if (normal.x == 0 && normal.y == 0 && normal.z == 0)
    goto noavg; // couldn't get normal vector for this point, punt

By contrast, calling normal_to_object(p, block) with the same point returns {0, 0, 1}, and subpixel averaging is performed in the pixel surrounding p. Hopefully @HomerReid can take a look at this.

[HomerReid]

@ChristopherHogan, this is officially the most insanely useful bug report in history. It would have taken me hours to track down the issue you identified. It is fixed in stevengj/libctl#14. I verified that your reproducer code fails in the current libctl master, but succeeds with the fix. I also updated the test-prisms unit test in libctl to test normals for points on prism surfaces; the unit test fails in master but succeeds with the fix.

[ChristopherHogan]

With stevengj/libctl#14, the prism and block runs are accurate to a tolerance of 1e-7. However, to get this accuracy I had to change the scheme bend-flux.ctl example by removing the infinity occurrences in the block data, re-running bend-flux.ctl to get the reference results, and removing mp.inf from the python block data. Without these steps the prism run only matches the block run to a tolerance of 1e-2. Here is a reproducer for that case (if it matters):

int main(int argc, char *argv[])
{
  meep_geom::material_type mat = meep_geom::make_dielectric(12);
  vector3 center = {0, -11.5, 0};
  vector3 e1 = {1, 0, 0};
  vector3 e2 = {0, 1, 0};
  vector3 e3 = {0, 0, 1};
  vector3 size = {1e20, 1, 1e20};
  geometric_object block = make_block(mat, center, e1, e2, e3, size);

  vector3 vertices[4] = {
    {-8, -11, 0},
    {8, -11, 0},
    {8, -12, 0},
    {-8, -12, 0}
  };
  vector3 axis = {0, 0, 1};
  geometric_object prism = make_prism(mat, vertices, 4, 0, axis);

  vector3 p = {6.9, -12, 0};

  vector3 nhat_block = normal_to_object(p, block);
  vector3 nhat_prism = normal_to_object(p, prism);

  double tolerance = 1.0e-6;
  if (!vector3_nearly_equal(nhat_block, nhat_prism, tolerance)) {
    printf("FAIL: Normals not equal\n");
    printf("  Block: {%f, %f, %f}\n", nhat_block.x, nhat_block.y, nhat_block.z);
    printf("  Prism: {%f, %f, %f}\n", nhat_prism.x, nhat_prism.y, nhat_prism.z);
  }

  return 0;
}

[oskooi]

Changing the z size of the blocks from infinity to 0 is not correct since the structure is supposed to be 2d and thus infinitely extended in the z direction. A size of 0 corresponds to a single pixel; this is not the same as a size of infinity. Also, a z size of infinity makes computing the normal unambiguous: it must lie in the xy plane.

[ChristopherHogan]

Ok, I reverted the scheme example back to using infinite sized blocks. The question is, how do we reproduce that geometry with prisms? I tried this:

no_bend_vertices = [
  mp.Vector3(-mp.inf, wvg_ycen - 0.5 * w),
  mp.Vector3(+mp.inf, wvg_ycen - 0.5 * w),
  mp.Vector3(+mp.inf, wvg_ycen + 0.5 * w),
  mp.Vector3(-mp.inf, wvg_ycen + 0.5 * w)
]

geometry = [mp.Prism(no_bend_vertices, mp.inf, material=mp.Medium(epsilon=12))]

but it produced a completely black epsilon file (all ones).

[ChristopherHogan]

I should also add that the reason I had replaced the infinite sizes with zeros was because that's what #341 did. I guess it needs to be fixed too.

[ChristopherHogan]

Extending the prism a little past the cell and giving it a height of 100 made the no-bend initialization run accurate to the block geometry up to 1e-7. Changing the height to 100 on the bend run helped a little, but it's still only accurate to 1e-3. Hopefully this is good enough. The only thing I noticed is a couple strange artifacts in the corner of the block bend (left), whereas the prism (right) is completely smooth.

[stevengj]

How do the errors in the bend case (compared to running at a very high resolution) compare for the block vs prism case now? Are they both on the same order?

[ChristopherHogan]

Yes, it looks better now.
bend_flux2.xlsx

[ChristopherHogan]

Ardavan requested that I update the mode-decomposition.py example to use prisms. I have a branch where that is working, but I can't open a PR until this is merged so it would be nice to get this merged before Thursday so there's time to rebase. That is, if you're happy with the numbers I posted above in the bend_flux2 spreadsheet.