"discrete" option in python/visualization.plot_eps

[simbilod]

Small PR for support for "discrete" plotting options in python/visualization.plot_eps

This is useful if a geometry has both large and small index contrasts. Below, for instance, we have ncore=3.45, nclad=1.44, and nfibercore=1.4467. sim.plot2D() with default or contour options cannot resolve the fiber core:

eps_parameters = {}
sim.plot2D(eps_parameters=eps_parameters)

eps_parameters = {}
eps_parameters['contour'] = True
sim.plot2D(eps_parameters=eps_parameters)

In this PR, the eps_parameters option discrete (a boolean like contour) allows plotting of user-provided permittivity ranges in solid colors:

# Make an iterable with some permittivities of the simulation geometry
epsilons = []
for item in geometry:
    epsilons.append(item.material.epsilon(1 / 1.55))
epsilons = np.array(epsilons)[:, 0, 0]

eps_parameters = {}
eps_parameters['discrete'] = True
eps_parameters['discrete_eps_levels'] = epsilons
sim.plot2D(eps_parameters=eps_parameters)

This is controlled by new optional entries of the eps_parameters dict:

• discrete: if True (and contour False), enables this option. Default is False
• discrete_eps_levels: an iterable with the permittivities to highlight. The code actually sets the coloring thresholds at the midpoints of the unique, sorted values of this iterable. If left unset, will use the minimum and maximum of permittivity data across the simulation region, and threshold at the average.
• discrete_eps_colors: an iterable with the colors to use for each permittivity level (sorted in increasing order). Needs to be larger than the iterable corresponding to unique(discrete_eps_levels), otherwise will throw an AssertionError. If left, unset, will generate len(discrete_eps_levels) colors equally distributed from the cmap item, which by default is the binary matplotlib colormap.

Examples with all options are below. The effect of subpixel averaging can be seen.

Using another colormap:

epsilons = []
for item in geometry:
    epsilons.append(item.material.epsilon(1 / 1.55))
epsilons = np.array(epsilons)[:, 0, 0]

eps_parameters = {}
eps_parameters['discrete'] = True
eps_parameters['discrete_eps_levels'] = epsilons
eps_parameters['cmap'] = 'viridis'
sim.plot2D(eps_parameters=eps_parameters)

Using an iterable of colors:

epsilons = []
for item in geometry:
    epsilons.append(item.material.epsilon(1 / 1.55))
epsilons = np.array(epsilons)[:, 0, 0]

eps_parameters = {}
eps_parameters['discrete'] = True
eps_parameters['discrete_eps_levels'] = epsilons
eps_parameters['discrete_eps_colors'] = ['red', 'blue', 'green', 'yellow']
sim.plot2D(eps_parameters=eps_parameters)

[codecov-commenter]

Codecov Report

Merging #1867 (610e823) into master (1e43687) will decrease coverage by 0.32%.
The diff coverage is 14.81%.

@@            Coverage Diff             @@
##           master    #1867      +/-   ##
==========================================
- Coverage   73.40%   73.08%   -0.33%     
==========================================
  Files          17       17              
  Lines        4896     4923      +27     
==========================================
+ Hits         3594     3598       +4     
- Misses       1302     1325      +23     

Impacted Files	Coverage Δ
python/visualization.py	44.19% <14.81%> (-1.47%)	⬇️

[simbilod]

Looking back I realize one could probably do all this processing outside of the plot function, and just pass the right norm and cmap arguments. But maybe this is useful for others

[oskooi]

It should be possible to generate a contour plot (via setting eps_parameters['contour'] to True) for any number of ε values present in the simulation. In the example you showed above, it seems strange why the fiber core does not appear in the image. The user shouldn't need to pass in a separate list of ε values. If necessary, this list should be determined automatically from the default_material property of the Simulation object as well as the geometry objects.

(Note that #1827 revamped plot2D so that plot_eps obtains its data directly from the geometric objects rather than the Yee grid which requires the structure object to be initialized. As such, no subpixel smoothing or interpolation to the centered grid is involved which could introduce intermediate ε values.)

[simbilod]

Ah, I'd missed the switch to geometry-defined permittivity visualization. No intermediate epsilons and just pulling the unique levels from the geometry makes this a lot easier.

Now that you mention it I am indeed able to get the contour to work if I explicitely give the levels. But in my edge case it does require the user to provide the levels.

epsilons = []
for item in geometry:
    epsilons.append(item.material.epsilon(1 / 1.55))
epsilons = np.array(epsilons)[:, 0, 0]

eps_parameters = {}
eps_parameters['contour'] = True
eps_parameters['levels'] = np.unique(epsilons)
sim.plot2D(eps_parameters=eps_parameters)

An easy fix would be to force it in the plotting function, e.g.

ax.contour(eps_data, 0, levels=np.unique(eps_data), colors='black', origin='upper', extent=extent, linewidths=eps_parameters['contour_linewidth'])

And then if discrete color levels is still useful it can also pull the epsilons from there.

[oskooi]

Here is a 2d test adapted from the one in #1827 (comment) which consists of various geometric objects with different values of ε. The figure below shows the output from running this script when the contour property of eps_parameters is set to either False or True. The image plot (left figure) is correct but the contour plot is not. The fix, as you suggest, is to specify the levels parameter of the contour function as a list of all the ε values. However, this needs to be done internally within plot_eps of visualization.py rather than as a user-provided input as a new dictionary key. Regardless, the image plot feature is correct for intermediate values of ε and perhaps you probably just need to be using the latest version of plot2D from the master branch.

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

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

design_region_resolution = 500.
design_shape = mp.Vector3(1.,1.,0)
Nx = int(design_region_resolution*design_shape.x) + 1
Ny = int(design_region_resolution*design_shape.y) + 1
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)
rad = 0.201943
w = 0.104283
design_params = np.where(np.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                                        np.sqrt(np.square(xv) + np.square(yv)) < rad+w),
                         1.,
                         0.)

matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                          mp.air,
                          mp.Medium(index=3.5),
                          weights=design_params,
                          do_averaging=False)

geometry = [mp.Cylinder(center=mp.Vector3(0.35,0.1),
                        radius=0.1,
                        height=0.5,
                        material=mp.Medium(index=1.5)),
            mp.Block(center=mp.Vector3(-0.15,-0.2),
                     size=mp.Vector3(0.2,0.24,0.8),
                     material=mp.Medium(index=1.5431)),
            mp.Block(center=mp.Vector3(-0.2,0.2),
                     size=mp.Vector3(0.4,0.4,0.3),
                     material=matgrid),
            mp.Prism(vertices=[mp.Vector3(0.05,0.45),
                               mp.Vector3(0.32,0.22),
	                       mp.Vector3(0.15,0.10)],
                     height=0.5,
                     material=mp.Medium(index=1.3135))]

sources = [mp.Source(mp.ContinuousSource(1.0),
	             center=mp.Vector3(),
                     size=mp.Vector3(0.2,0.2,0.2),
                     component=mp.Ez)]

sim = mp.Simulation(resolution=20,
                    cell_size=cell_size,
	            geometry=geometry,
                    sources=sources,
                    boundary_layers=[mp.PML(thickness=0.1,side=mp.High,direction=mp.Y)])

mon = sim.add_dft_fields([mp.Ez],
                         1.0,
                         0,
                         1,
                         center=mp.Vector3(0.3,-0.2),
                         size=mp.Vector3(0.3,0.1))

sim.plot2D(eps_parameters={'contour':True, 'resolution':400},
           output_plane=mp.Volume(center=mp.Vector3(),
                                  size=mp.Vector3(1.,1.,0)))
plt.title('contour: True')
plt.gca().set_aspect('equal','box')
plt.savefig('plot2D_test_contour.png',
            dpi=150,
            bbox_inches='tight')

[simbilod]

Thank you for taking a look. I do get this output when I run this example, and my software is up-to-date.

I think my problem is simply that, by default, imshow's norm parameter uses a linear scaling of value to color. So for tiny epsilon differences such as for my single-mode fiber and cladding, in conjunction with a larger one, the small difference is unnoticeable.

To avoid overcomplicating the function, discrete colors like I implemented can be done outside by the user, but to improve contour plotting I can do another PR with enforced levels for the contour option.