"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/python/tests/mode_coeffs.py" (16 Mar 2019, 5600 Bytes) of package /linux/privat/meep-1.10.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "mode_coeffs.py": 1.8.0_vs_1.9.0.

    1 from __future__ import division
    2 
    3 import unittest
    4 import numpy as np
    5 import meep as mp
    6 
    7 
    8 class TestModeCoeffs(unittest.TestCase):
    9 
   10     def run_mode_coeffs(self, mode_num, kpoint_func, nf=1, resolution=15):
   11 
   12         w = 1   # width of waveguide
   13         L = 10  # length of waveguide
   14 
   15         Si = mp.Medium(epsilon=12.0)
   16 
   17         dair = 3.0
   18         dpml = 3.0
   19 
   20         sx = dpml + L + dpml
   21         sy = dpml + dair + w + dair + dpml
   22         cell_size = mp.Vector3(sx, sy, 0)
   23 
   24         prism_x = sx + 1
   25         prism_y = w / 2
   26         vertices = [mp.Vector3(-prism_x, prism_y),
   27                     mp.Vector3(prism_x, prism_y),
   28                     mp.Vector3(prism_x, -prism_y),
   29                     mp.Vector3(-prism_x, -prism_y)]
   30 
   31         geometry = [mp.Prism(vertices, height=mp.inf, material=Si)]
   32 
   33         boundary_layers = [mp.PML(dpml)]
   34 
   35         # mode frequency
   36         fcen = 0.20  # > 0.5/sqrt(11) to have at least 2 modes
   37         df   = 0.5*fcen
   38 
   39         source=mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=df),
   40                                   eig_band=mode_num,
   41                                   size=mp.Vector3(0,sy-2*dpml,0),
   42                                   center=mp.Vector3(-0.5*sx+dpml,0,0),
   43                                   eig_match_freq=True,
   44                                   eig_resolution=2*resolution)
   45 
   46         sim = mp.Simulation(resolution=resolution,
   47                             cell_size=cell_size,
   48                             boundary_layers=boundary_layers,
   49                             geometry=geometry,
   50                             sources=[source],
   51                             symmetries=[mp.Mirror(mp.Y, phase=1 if mode_num % 2 == 1 else -1)])
   52 
   53         xm = 0.5*sx - dpml  # x-coordinate of monitor
   54         mflux = sim.add_mode_monitor(fcen, df, nf, mp.ModeRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml)))
   55         mode_flux = sim.add_flux(fcen, df, nf, mp.FluxRegion(center=mp.Vector3(xm,0), size=mp.Vector3(0,sy-2*dpml)))
   56 
   57         # sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(-0.5*sx+dpml,0), 1e-10))
   58         sim.run(until_after_sources=200)
   59 
   60         ##################################################
   61         # If the number of analysis frequencies is >1, we
   62         # are testing the unit-power normalization
   63         # of the eigenmode source: we observe the total
   64         # power flux through the mode_flux monitor (which
   65         # equals the total power emitted by the source as
   66         # there is no scattering in this ideal waveguide)
   67         # and check that it agrees with the prediction
   68         # of the eig_power() class method in EigenmodeSource.
   69         ##################################################
   70         if nf>1:
   71             power_observed=mp.get_fluxes(mode_flux)
   72             freqs=mp.get_flux_freqs(mode_flux)
   73             power_expected=[source.eig_power(f) for f in freqs]
   74             return freqs, power_expected, power_observed
   75 
   76         modes_to_check = [1, 2]  # indices of modes for which to compute expansion coefficients
   77         res = sim.get_eigenmode_coefficients(mflux, modes_to_check, kpoint_func=kpoint_func)
   78 
   79         self.assertTrue(res.kpoints[0].close(mp.Vector3(0.604301, 0, 0)))
   80         self.assertTrue(res.kpoints[1].close(mp.Vector3(0.494353, 0, 0), tol=1e-2))
   81         self.assertTrue(res.kdom[0].close(mp.Vector3(0.604301, 0, 0)))
   82         self.assertTrue(res.kdom[1].close(mp.Vector3(0.494353, 0, 0), tol=1e-2))
   83 
   84         mode_power = mp.get_fluxes(mode_flux)[0]
   85 
   86         TestPassed = True
   87         TOLERANCE = 5.0e-3
   88         c0 = res.alpha[mode_num - 1, 0, 0] # coefficient of forward-traveling wave for mode #mode_num
   89         for nm in range(1, len(modes_to_check)+1):
   90             if nm != mode_num:
   91                 cfrel = np.abs(res.alpha[nm - 1, 0, 0]) / np.abs(c0)
   92                 cbrel = np.abs(res.alpha[nm - 1, 0, 1]) / np.abs(c0)
   93                 if cfrel > TOLERANCE or cbrel > TOLERANCE:
   94                     TestPassed = False
   95 
   96         self.sim = sim
   97 
   98         # test 1: coefficient of excited mode >> coeffs of all other modes
   99         self.assertTrue(TestPassed, msg="cfrel: {}, cbrel: {}".format(cfrel, cbrel))
  100         # test 2: |mode coeff|^2 = power
  101         self.assertAlmostEqual(mode_power / abs(c0**2), 1.0, places=1)
  102 
  103         return res
  104 
  105     def test_modes(self):
  106         self.run_mode_coeffs(1, None)
  107         res = self.run_mode_coeffs(2, None)
  108 
  109         # Test mp.get_eigenmode and EigenmodeData
  110         vol = mp.Volume(center=mp.Vector3(5), size=mp.Vector3(y=7))
  111         emdata = self.sim.get_eigenmode(0.2, mp.X, vol, 2, mp.Vector3())
  112         self.assertEqual(emdata.freq, 0.2)
  113         self.assertEqual(emdata.band_num, 2)
  114         self.assertTrue(emdata.kdom.close(res.kdom[1]))
  115 
  116         eval_point = mp.Vector3(0.7, -0.2, 0.3)
  117         ex_at_eval_point = emdata.amplitude(eval_point, mp.Ex)
  118         hz_at_eval_point = emdata.amplitude(eval_point, mp.Hz)
  119         self.assertAlmostEqual(ex_at_eval_point, 0.45358518109307083+0.5335421986481814j)
  120         self.assertAlmostEqual(hz_at_eval_point, 3.717865162096829-3.1592989829386298j)
  121 
  122     def test_kpoint_func(self):
  123 
  124         def kpoint_func(freq, mode):
  125             return mp.Vector3()
  126 
  127         self.run_mode_coeffs(1, kpoint_func)
  128 
  129     def test_eigensource_normalization(self):
  130         f, p_exp, p_obs=self.run_mode_coeffs(1, None, nf=51, resolution=15)
  131         #self.assertAlmostEqual(max(p_exp),max(p_obs),places=1)
  132         max_exp, max_obs=max(p_exp), max(p_obs)
  133         self.assertLess(abs(max_exp-max_obs), 0.5*max(abs(max_exp),abs(max_obs)))
  134 
  135 
  136 if __name__ == '__main__':
  137     unittest.main()