plot geometry for dispersive materials without initializing structure object

[oskooi]

Closes #1824.

Removed the frequency parameter from plot2D since this cannot be passed to get_epsilon_grid. As a result, this API change is not backwards compatible. As a workaround, we can leave the frequency parameter in but just mention in the documentation that it is simply ignored.

cc @joamatab

[codecov-commenter]

Codecov Report

Merging #1827 (00f7b48) into master (03ea7ad) will decrease coverage by 0.09%.
The diff coverage is 63.63%.

@@            Coverage Diff             @@
##           master    #1827      +/-   ##
==========================================
- Coverage   73.53%   73.44%   -0.10%     
==========================================
  Files          17       17              
  Lines        4879     4891      +12     
==========================================
+ Hits         3588     3592       +4     
- Misses       1291     1299       +8     

Impacted Files	Coverage Δ
python/visualization.py	45.72% <62.85%> (-1.05%)	⬇️
python/simulation.py	76.87% <80.00%> (+0.16%)	⬆️

[oskooi]

For the record, there was a single test failure unrelated to this PR that went away after rerunning. The symmetry unit test of test_material_grid.py was failing a check in the sixth decimal digit for single precision:

======================================================================
FAIL: test_symmetry (__main__.TestMaterialGrid)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/runner/work/meep/meep/build/meep-1.22.0-beta/_build/sub/python/../../../python/tests/test_material_grid.py", line 161, in test_symmetry
    self.assertAlmostEqual(tran_nosym, tran_sym)
AssertionError: 20.27522531963953 != 20.275228328448193 within 7 places (3.0088086617752197e-06 difference)

It might therefore be useful to increase the tolerance slightly for this test to avoid this failure in the future.

[oskooi]

To demonstrate that this feature is working correctly, here is a 2d test case which contains a variety of different simulation objects: two types of geometric objects (Cylinder and MaterialGrid), a source, PML in a single direction, and a DFT monitor. plot2D is called with three different resolution values of 20, 80, and 320 for sampling the ε grid. A single-threaded run takes about ~0.5 s because the structure object is never initialized.

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

mp.verbosity(3)

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

design_region_resolution = 500
design_shape = mp.Vector3(1.0,1.0,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
design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad

matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                          mp.air,
                          mp.Medium(index=3.5),
                          weights=design_params,
                          do_averaging=True,
                          beta=1000,
                          eta=0.5)

geometry = [mp.Cylinder(center=mp.Vector3(0.35,0.1),
                        radius=0.1,
                        height=mp.inf,
                        material=mp.Medium(index=1.5)),
            mp.Block(center=mp.Vector3(-0.15,-0.2),
                     size=mp.Vector3(0.2,0.24,mp.inf),
                     material=mp.Medium(index=2.5)),
            mp.Block(center=mp.Vector3(-0.2,0.2),
                     size=mp.Vector3(0.4,0.4,mp.inf),
                     material=matgrid)]

sources = [mp.Source(mp.ContinuousSource(1.0),
                     center=mp.Vector3(),
                     size=mp.Vector3(0.2,0.2,0),
                     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))

grid_resolution = 320
sim.plot2D(resolution=grid_resolution)
plt.title('resolution = {}'.format(grid_resolution))
plt.savefig('ring_matgrid_res{}'.format(grid_resolution),
            dpi=150,
            bbox_inches='tight')

output

Using MPI version 3.1, 1 processes
     cylinder, center = (0.35,0.1,0)
          radius 0.1, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (2.25,2.25,2.25)
     block, center = (-0.15,-0.2,0)
          size (0.2,0.24,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (6.25,6.25,6.25)
     block, center = (-0.2,0.2,0)
          size (0.4,0.4,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
Geometric-object bounding-box tree:
     box (-0.5..0.5, -0.5..0.5, 0..0)
          bounding box (-0.4..0, 0..0.4, -5e+19..5e+19)
          shift object by (0, 0, 0)
          block, center = (-0.2,0.2,0)
               size (0.4,0.4,1e+20)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (-0.25..-0.05, -0.32..-0.08, -5e+19..5e+19)
          shift object by (0, 0, 0)
          block, center = (-0.15,-0.2,0)
               size (0.2,0.24,1e+20)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (0.25..0.45, 0..0.2, -5e+19..5e+19)
          shift object by (0, 0, 0)
          cylinder, center = (0.35,0.1,0)
               radius 0.1, height 1e+20, axis (0, 0, 1)
Geometric object tree has depth 1 and 3 object nodes (vs. 3 actual objects)

Elapsed run time = 0.5704 s

[smartalecH]

Is it just as fast when plotting a 2D cross section of a 3D simulation?

[oskooi]

Is it just as fast when plotting a 2D cross section of a 3D simulation?

Taking the example above and extending it to 3d, the difference in runtimes at a resolution of 320 for this PR compared to master is 0.5773 s vs. 87.5444 s. This translates to a speedup of a factor of more than 150. Using the revamped plot2D in this PR to generate a 2D cross section of a 3D cell seems to be just as fast as the 2D cell shown previously (both roughly ~0.5 s).

PR 1827

Using MPI version 3.1, 1 processes
     cylinder, center = (0.35,0.1,0)
          radius 0.1, height 0.5, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (2.25,2.25,2.25)
     block, center = (-0.15,-0.2,0)
          size (0.2,0.24,0.8)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (6.25,6.25,6.25)
     block, center = (-0.2,0.2,0)
          size (0.4,0.4,0.3)
          axes (1,0,0), (0,1,0), (0,0,1)
Geometric-object bounding-box tree:
     box (-0.5..0.5, -0.5..0.5, 0..0)
          bounding box (-0.4..0, 0..0.4, -0.15..0.15)
          shift object by (0, 0, 0)
          block, center = (-0.2,0.2,0)
               size (0.4,0.4,0.3)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (-0.25..-0.05, -0.32..-0.08, -0.4..0.4)
          shift object by (0, 0, 0)
          block, center = (-0.15,-0.2,0)
               size (0.2,0.24,0.8)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (0.25..0.45, 0..0.2, -0.25..0.25)
          shift object by (0, 0, 0)
          cylinder, center = (0.35,0.1,0)
               radius 0.1, height 0.5, axis (0, 0, 1)
Geometric object tree has depth 1 and 3 object nodes (vs. 3 actual objects)

Elapsed run time = 0.5773 s

master

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
time for choose_chunkdivision = 3.6012e-05 s
Working in 3D dimensions.
Computational cell is 1 x 1 x 1 with resolution 320
     cylinder, center = (0.35,0.1,0)
          radius 0.1, height 0.5, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (2.25,2.25,2.25)
     block, center = (-0.15,-0.2,0)
          size (0.2,0.24,0.8)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (6.25,6.25,6.25)
     block, center = (-0.2,0.2,0)
          size (0.4,0.4,0.3)
          axes (1,0,0), (0,1,0), (0,0,1)
Geometric-object bounding-box tree:
     box (-0.503125..0.503125, -0.503125..0.503125, -0.503125..0.503125)
          bounding box (-0.4..0, 0..0.4, -0.15..0.15)
          shift object by (0, 0, 0)
          block, center = (-0.2,0.2,0)
               size (0.4,0.4,0.3)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (-0.25..-0.05, -0.32..-0.08, -0.4..0.4)
          shift object by (0, 0, 0)
          block, center = (-0.15,-0.2,0)
               size (0.2,0.24,0.8)
               axes (1,0,0), (0,1,0), (0,0,1)
          bounding box (0.25..0.45, 0..0.2, -0.25..0.25)
          shift object by (0, 0, 0)
          cylinder, center = (0.35,0.1,0)
               radius 0.1, height 0.5, axis (0, 0, 1)
Geometric object tree has depth 1 and 3 object nodes (vs. 3 actual objects)
subpixel-averaging is 17.152% done, 19.3209 s remaining
subpixel-averaging is 30.5805% done, 9.08038 s remaining
subpixel-averaging is 43.1867% done, 5.26365 s remaining
subpixel-averaging is 59.4996% done, 2.72321 s remaining
subpixel-averaging is 77.5109% done, 1.16075 s remaining
subpixel-averaging is 91.4617% done, 0.374134 s remaining
subpixel-averaging is 17.1149% done, 19.3729 s remaining
subpixel-averaging is 30.3951% done, 9.1621 s remaining
subpixel-averaging is 42.9373% done, 5.31743 s remaining
subpixel-averaging is 59.0481% done, 2.77456 s remaining
subpixel-averaging is 77.1773% done, 1.18323 s remaining
subpixel-averaging is 91.0708% done, 0.392198 s remaining
subpixel-averaging is 17.1217% done, 19.3657 s remaining
subpixel-averaging is 30.5299% done, 9.1024 s remaining
subpixel-averaging is 43.2002% done, 5.26053 s remaining
subpixel-averaging is 59.4356% done, 2.7301 s remaining
subpixel-averaging is 77.4469% done, 1.16553 s remaining
subpixel-averaging is 91.2022% done, 0.386007 s remaining
time for set_epsilon = 86.1741 s
-----------

Elapsed run time = 87.5444 s

[stevengj]

Removed the frequency parameter from plot2D since this cannot be passed to get_epsilon_grid.

Why can't we update get_epsilon_grid to accept this?

[oskooi]

Why can't we update get_epsilon_grid to accept this?

This seems like it is possible without much additional effort since get_epsilon_grid can just call a modified version of the recently added function eff_chi1inv_row_disp (from #1801) to compute the trace of the dispersive ε tensor rather than the current function chi1p1 which computes the instantaneous ε.

[stevengj]

See how array slice gets a scalar quantity from the trace of inv(ε):

meep/src/array_slice.cpp

Lines 379 to 387 in 03ea7ad

	else if (cS[i] == Dielectric) {
	complex<double> tr(0.0, 0.0);
	for (int k = 0; k < data->ninveps; ++k) {
	tr += (fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx, frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[1 + 2 * k], frequency) +
	fc->s->get_chi1inv_at_pt(iecs[k], ieds[k], idx + ieos[2 * k] + ieos[1 + 2 * k],
	frequency));
	if (abs(tr) == 0.0) tr += 4.0; // default inveps == 1

[oskooi]

Pull out and refactor:

meep/src/meepgeom.cpp

Lines 2568 to 2588 in 03ea7ad

	// loop over all the tensor components
	for (int i = 0; i < 9; i++) {
	std::complex<double> a,b;
	// compute first part containing conductivity
	vector3 dummy;
	dummy.x = dummy.y = dummy.z = 0.0;
	double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i);
	a = std::complex<double>(1.0, conductivityCur / (freq));
	// compute lorentzian component
	b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i);
	for (const auto &mm_susc: mm->E_susceptibilities) {
	meep::lorentzian_susceptibility sus = meep::lorentzian_susceptibility(
	mm_susc.frequency, mm_susc.gamma, mm_susc.drude);
	double sigma = vec_to_value(mm_susc.sigma_diag, mm_susc.sigma_offdiag, i);
	b += sus.chi1(freq, sigma);
	}
	// elementwise multiply
	tensor[i] = a * b;
	}

[stevengj]

Refactor that code into new function get_chi1_tensor_disp(meep::component c, std::complex<double> tensor[9], const meep::vec &r, double freq, geom_epsilon *geps), which is then called by eff_chi1inv_row_disp. Your code will call get_chi1_tensor_disp directly, and then compute the trace/3 as the average "scalar" ε (this is the average eigenvalue).

[oskooi]

Added support for dispersive materials following the approach we discussed above. This involved reinstating the frequency parameter of plot2D and adding a new frequency parameter to get_epsilon_grid. The test case below for an arbitrary collection of dispersive and non-dispersive material geometries demonstrates a speedup of more than a factor of 600 times compared to master.

import meep as mp
from meep.materials import SiN, Co
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

mp.verbosity(3)

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

design_region_resolution = 500
design_shape = mp.Vector3(1.0,1.0,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.logical_and(np.sqrt(np.square(xv) + np.square(yv)) > rad,
                               np.sqrt(np.square(xv) + np.square(yv)) < rad+w,
                               dtype=np.double)

matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
                          mp.air,
                          mp.Medium(index=3.5),
                          weights=design_params,
                          do_averaging=True,
                          beta=1000,
                          eta=0.5)

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=SiN),
            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=Co)]

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=400,
                    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))

grid_resolution = 400
frequency = 0
sim.plot2D(resolution=grid_resolution,
           frequency=frequency,
           output_plane=mp.Volume(center=mp.Vector3(),
                                  size=mp.Vector3(1,1,0)))
plt.title('frequency = {}, resolution = {}'.format(frequency,grid_resolution))
plt.savefig('ring_matgrid_3d_dispersive_res{}_freq{}_master'.format(grid_resolution,frequency),
            dpi=150,
            bbox_inches='tight')