Artifacts around slanted prism objects

[danielwboyce]

See comments on NanoComp/meep#1129. Artifacts and fuzziness appear along the prism boundary with nonzero values of sidewall_angle. See these comments in particular (1)(2)(3).

Some strange artifacts are found in the epsilon profile of a hexagonal prism in the xy plane of a 2d cell whenever the sidewall_angle is not zero as shown below.

import meep as mp
import numpy as np
import math
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

cell_size = mp.Vector3(2.5,2.5)

resolution = 50

vertices = [mp.Vector3(-1,0),
            mp.Vector3(-0.5,math.sqrt(3)/2),
            mp.Vector3(0.5,math.sqrt(3)/2),
            mp.Vector3(1,0),
            mp.Vector3(0.5,-math.sqrt(3)/2),
            mp.Vector3(-0.5,-math.sqrt(3)/2)]

geometry = [mp.Prism(vertices,                     
                     height=1.5,
                     center=mp.Vector3(),
                     sidewall_angle=np.radians(15),
                     material=mp.Medium(index=3.5))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    geometry=geometry)

sim.init_sim()
plt.figure()
sim.plot2D()
plt.savefig('prism_sidewall.png')

sidewall_angle=0

sidewall_angle=np.radians(15)

sidewall_angle=np.radians(27)

Note that the size of the cross section is shrinking as the sidewall_angle increases. This is correct given the intersection of the prism with the z=0 plane.

Artifacts along the prism boundaries also appear in 3d.

import meep as mp
import numpy as np
import math
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

cell_size = mp.Vector3(3,3,3)

resolution = 50

vertices = [mp.Vector3(-1,0),
            mp.Vector3(-0.5,math.sqrt(3)/2),
            mp.Vector3(0.5,math.sqrt(3)/2),
            mp.Vector3(1,0),
            mp.Vector3(0.5,-math.sqrt(3)/2),
            mp.Vector3(-0.5,-math.sqrt(3)/2)]

geometry = [mp.Prism(vertices,
                     axis=mp.Vector3(0,0,1),
                     height=1.5,
                     center=mp.Vector3(),
                     sidewall_angle=np.radians(10),
                     material=mp.Medium(index=3.5))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    geometry=geometry)

sim.run(mp.at_beginning(mp.output_epsilon),until=0)

$ h5topng -o eps-xy.png -z 75 prism_sidewall-eps-000000.00.h5
$ h5topng -o eps-xz.png -y 75 prism_sidewall-eps-000000.00.h5
$ h5topng -o eps-yz.png -x 75 prism_sidewall-eps-000000.00.h5

eps-xy.png

eps-xz.png

eps-yz.png

When the sidewall angle is 0, no artifacts appear along the boundaries.

eps-xy.png

eps-xz.png

eps-yz.png

[danielwboyce]

As a quick check I added this test to utils/test-prism.c (constructing the same prism defined in the above tests) to see if libctl is having issues with this geometry for some reason.

/************************************************************************/
/* 8th unit test: quick test of hexagon test prism to see if the        */
/* vertices are calculated correctly                                    */
/************************************************************************/
int test_hex_prism() {
  void *m = NULL;

  vector3_list nodes;
  nodes.num_items = 6;
  nodes.items = (vector3 *) malloc(nodes.num_items * sizeof(vector3));
  nodes.items[0] = make_vector3(-1.0, 0.0, 0.0);
  nodes.items[1] = make_vector3(-0.5, sqrt(3.0) / 2.0, 0.0);
  nodes.items[2] = make_vector3(0.5, sqrt(3.0) / 2.0, 0.0);
  nodes.items[3] = make_vector3(1.0, 0.0, 0.0);
  nodes.items[4] = make_vector3(0.5, -sqrt(3.0) / 2.0, 0.0);
  nodes.items[5] = make_vector3(-0.5, -sqrt(3.0) / 2.0, 0.0);

  double height = 1.5;
  vector3 zhat = make_vector3(0, 0, 1);

  double ten_degree_sidewall = 10.0 * K_PI / 180.0;
  geometric_object hex_ten_degree_sidewall_geom_object = make_slanted_prism(m, nodes, nodes.num_items, height, zhat,
                                                                            ten_degree_sidewall);
  prism *hex_ten_degree_sidewall_prism = hex_ten_degree_sidewall_geom_object.subclass.prism_data;

  prism2gnuplot(hex_ten_degree_sidewall_prism, "hex_ten_degree_sidewall_gnu_plot.dat");

  return 0;
}

Plotting the calculated vertices with gnuplot, we see that the vertices look as expected. So it seems that libctl doesn't have an inherent problem with this geometry, but that the problem is the way meep calls libctl or just something in meep internals.