unit test for conductivity (NanoComp#1857)

* unit test for conductivity

* describe in the docs how to model the attenutation coefficient using conductivity

* Update python/tests/test_conductivity.py

Co-authored-by: Steven G. Johnson <stevenj@mit.edu>

doc/docs/Materials.md

@@ -88,13 +88,15 @@ Meep can keep track of this energy for the Lorentzian polarizability terms but n
 Conductivity and Complex ε
 --------------------------
 
-Often, you only care about the absorption loss in a narrow bandwidth, where you just want to set the imaginary part of ε (or μ) to some known experimental value, in the same way that you often just care about setting a dispersionless real ε that is the correct value in your bandwidth of interest.
+Often, you only care about the absorption loss in a narrow bandwidth, where you just want to set the imaginary part of $\varepsilon$ (or $\mu$) to some known experimental value, in the same way that you often just care about setting a dispersionless real $\varepsilon$ that is the correct value in your bandwidth of interest.
 
 One approach to this problem would be allowing you to specify a constant, frequency-independent, imaginary part of $\varepsilon$, but this has the disadvantage of requiring the simulation to employ complex fields which double the memory and time requirements, and also tends to be numerically unstable. Instead, the approach in Meep is for you to set the conductivity $\sigma_D$ (or $\sigma_B$ for an imaginary part of $\mu$), chosen so that $\mathrm{Im}\, \varepsilon = \varepsilon_\infty \sigma_D / \omega$ is the correct value at your frequency $\omega$ of interest. Note that, in Meep, you specify $f = \omega/2\pi$ instead of $\omega$ for the frequency, however, so you need to include the factor of $2\pi$ when computing the corresponding imaginary part of $\varepsilon$. Conductivities can be implemented with purely real fields, so they are not nearly as expensive as implementing a frequency-independent complex $\varepsilon$ or $\mu$.
 
-For example, suppose you want to simulate a medium with $\varepsilon = 3.4 + 0.101i$ at a frequency 0.42 (in your Meep units), and you only care about the material in a narrow bandwidth around this frequency (i.e. you don't need to simulate the full experimental frequency-dependent permittivity). Then, in Meep, you could use `meep.Medium(epsilon=3.4, D_conductivity=2*math.pi*0.42*0.101/3.4)` in Python or `(make medium (epsilon 3.4) (D-conductivity (* 2 pi 0.42 0.101 (/ 3.4))))` in Scheme; i.e. $\varepsilon_\infty = \mathrm{Re}[\varepsilon] = 3.4$ and $\sigma_D = \omega \, \mathrm{Im}[\varepsilon / \varepsilon_\infty] = (2\pi \, 0.42) \, 0.101 / 3.4$.
+For example, suppose you want to simulate a medium with $\varepsilon = 3.4 + 0.101i$ at a frequency 0.42 (in your Meep units), and you only care about the material in a narrow bandwidth around this frequency (i.e. you don't need to simulate the full experimental frequency-dependent permittivity). Then, in Meep, you could use `meep.Medium(epsilon=3.4, D_conductivity=2*math.pi*0.42*0.101/3.4)` in Python or `(make medium (epsilon 3.4) (D-conductivity (* 2 pi 0.42 0.101 (/ 3.4))))` in Scheme; i.e. $\varepsilon_\infty = \mathrm{Re}[\varepsilon] = 3.4$ and $\sigma_D = \omega \, \mathrm{Im}[\varepsilon] / \varepsilon_\infty = (2\pi \, 0.42) \, 0.101 / 3.4$.
 
-**Note**: the "conductivity" in Meep is slightly different from the conductivity you might find in a textbook, because for computational convenience it appears as $\sigma_D \mathbf{D}$ in our Maxwell equations rather than the more-conventional $\sigma \mathbf{E}$; this just means that our definition is different from the usual electric conductivity by a factor of ε. Also, just as Meep uses the dimensionless relative permittivity for ε, it uses nondimensionalized units of 1/*a* (where *a* is your unit of distance) for the conductivities $\sigma_{D,B}$. If you have the electric conductivity $\sigma$ in SI units and want to convert to $\sigma_D$ in Meep units, you can simply use the formula: $\sigma_D = (a/c) \sigma / \varepsilon_r \varepsilon_0$ where *a* is your unit of distance in meters, *c* is the vacuum speed of light in m/s, $\varepsilon_0$ is the SI vacuum permittivity, and $\varepsilon_r$ is the real relative permittivity.
+You can also use the $\sigma_D$ feature to model the [attenuation coefficient](https://en.wikipedia.org/wiki/Attenuation_coefficient) $\alpha$ (units of e.g. dB/cm) obtained from experimental measurements (i.e., ellipsometry). This involves first [converting $\alpha$ into a complex refractive index](https://en.wikipedia.org/wiki/Mathematical_descriptions_of_opacity#Complex_refractive_index) (which is then converted into a complex permittivity) with imaginary part given by $\lambda_0\alpha/(4\pi)$ where $\lambda_0$ is the vacuum wavelength.
+
+**Note**: the "conductivity" in Meep is slightly different from the conductivity you might find in a textbook, because for computational convenience it appears as $\sigma_D \mathbf{D}$ in our Maxwell equations rather than the more-conventional $\sigma \mathbf{E}$; this just means that our definition is different from the usual electric conductivity by a factor of $\varepsilon$. Also, just as Meep uses the dimensionless relative permittivity for $\varepsilon$, it uses nondimensionalized units of 1/*a* (where *a* is your unit of distance) for the conductivities $\sigma_{D,B}$. If you have the electric conductivity $\sigma$ in SI units and want to convert to $\sigma_D$ in Meep units, you can simply use the formula: $\sigma_D = (a/c) \sigma / \varepsilon_r \varepsilon_0$ where *a* is your unit of distance in meters, *c* is the vacuum speed of light in m/s, $\varepsilon_0$ is the SI vacuum permittivity, and $\varepsilon_r$ is the real relative permittivity.
 
 Nonlinearity
 ------------

python/Makefile.am

@@ -43,13 +43,14 @@ TESTS =                                        \
     $(TEST_DIR)/test_antenna_radiation.py      \
     $(TEST_DIR)/test_array_metadata.py         \
     $(TEST_DIR)/test_bend_flux.py              \
-    $(TEST_DIR)/test_binary_partition_utils.py       \
+    $(TEST_DIR)/test_binary_partition_utils.py \
     $(BINARY_GRATING_TEST)                     \
     $(TEST_DIR)/test_cavity_arrayslice.py      \
     $(TEST_DIR)/test_cavity_farfield.py        \
     $(TEST_DIR)/test_chunk_balancer.py         \
     $(TEST_DIR)/test_chunk_layout.py           \
     $(TEST_DIR)/test_chunks.py                 \
+    $(TEST_DIR)/test_conductivity.py           \
     $(TEST_DIR)/test_cyl_ellipsoid.py          \
     $(TEST_DIR)/test_dft_energy.py             \
     $(TEST_DIR)/test_dft_fields.py             \

python/tests/test_conductivity.py

@@ -0,0 +1,100 @@
+import unittest
+import numpy as np
+import meep as mp
+
+dB_cm_to_dB_um = 1e-4
+
+class TestConductivity(unittest.TestCase):
+
+    def wvg_flux(self, res, att_coeff):
+        """
+        Computes the Poynting flux in a single-mode waveguide at two
+        locations (5 and 10 μm) downstream from the source given the
+        grid resolution res (pixels/μm) and material attenuation
+        coefficient att_coeff (dB/cm).
+        """
+
+        cell_size = mp.Vector3(14.,14.)
+
+        pml_layers = [mp.PML(thickness=2.)]
+
+        w = 1. # width of waveguide
+
+        fsrc = 0.15 # frequency (in vacuum)
+
+        # note: MPB can only compute modes of lossless material systems.
+        #       The imaginary part of ε is ignored and the launched
+        #       waveguide mode is therefore inaccurate. For small values
+        #       of imag(ε) (which is proportional to att_coeff), this
+        #       inaccuracy tends to be insignificant.
+        sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc),
+                                      center=mp.Vector3(-5.),
+                                      size=mp.Vector3(y=10.),
+                                      eig_parity=mp.EVEN_Y+mp.ODD_Z)]
+
+        # ref: https://en.wikipedia.org/wiki/Mathematical_descriptions_of_opacity
+        # Note that this is the loss of a planewave, which is only approximately
+        # the loss of a waveguide mode.  In principle, we could compute the latter
+        # semi-analytically if we wanted to run this unit test to greater accuracy
+        # (e.g. to test convergence with resolution).
+        n_eff = np.sqrt(12.) + 1j * (1/fsrc) * (dB_cm_to_dB_um * att_coeff) / (4 * np.pi)
+        eps_eff = n_eff * n_eff
+        # ref: https://meep.readthedocs.io/en/latest/Materials/#conductivity-and-complex
+        sigma_D = 2 * np.pi * fsrc * np.imag(eps_eff) / np.real(eps_eff)
+
+        geometry = [mp.Block(center=mp.Vector3(),
+                             size=mp.Vector3(mp.inf,w,mp.inf),
+                             material=mp.Medium(epsilon=np.real(eps_eff),
+                                                D_conductivity=sigma_D))]
+
+        sim = mp.Simulation(cell_size=cell_size,
+                            resolution=res,
+                            boundary_layers=pml_layers,
+                            sources=sources,
+                            geometry=geometry,
+                            symmetries=[mp.Mirror(mp.Y)])
+
+        tran1 = sim.add_flux(fsrc,
+                             0,
+                             1,
+                             mp.FluxRegion(center=mp.Vector3(x=0.), size=mp.Vector3(y=10.)))
+
+        tran2 = sim.add_flux(fsrc,
+                             0,
+                             1,
+                             mp.FluxRegion(center=mp.Vector3(x=5.), size=mp.Vector3(y=10.)))
+
+        sim.run(until_after_sources=20)
+
+        return mp.get_fluxes(tran1)[0], mp.get_fluxes(tran2)[0]
+
+
+    def test_conductivity(self):
+        res = 25. # pixels/μm
+
+        # compute the incident flux for a lossless waveguide
+        incident_flux1, incident_flux2 = self.wvg_flux(res, 0.)
+        self.assertAlmostEqual(incident_flux1/incident_flux2, 1., places=2)
+        print("incident_flux:, {} (measured), 1.0 (expected)".format(incident_flux2/incident_flux1))
+
+        # compute the flux for a lossy waveguide
+        att_coeff = 37.46 # dB/cm
+        attenuated_flux1, attenuated_flux2 = self.wvg_flux(res, att_coeff)
+
+        L1 = 5.
+        expected_att1 = np.exp(-att_coeff * dB_cm_to_dB_um * L1)
+        self.assertAlmostEqual(attenuated_flux1/incident_flux2, expected_att1, places=2)
+        print("flux:, {}, {:.6f} (measured), {:.6f} (expected)".format(L1,
+                                                                       attenuated_flux1/incident_flux2,
+                                                                       expected_att1))
+
+        L2 = 10.
+        expected_att2 = np.exp(-att_coeff * dB_cm_to_dB_um * L2)
+        self.assertAlmostEqual(attenuated_flux2/incident_flux2, expected_att2, places=2)
+        print("flux:, {}, {:.6f} (measured), {:.6f} (expected)".format(L2,
+                                                                       attenuated_flux2/incident_flux2,
+                                                                       expected_att2))
+
+
+if __name__ == '__main__':
+    unittest.main()