"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/python/examples/solve-cw.py" (5 Jun 2019, 3195 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 latest Fossies "Diffs" side-by-side code changes report for "solve-cw.py": 1.9.0_vs_1.10.0.

    1 import meep as mp
    2 import numpy as np
    3 from numpy import linalg as LA
    4 import matplotlib.pyplot as plt
    5 
    6 n = 3.4
    7 w = 1
    8 r = 1
    9 pad = 4
   10 dpml = 2
   11 
   12 sxy = 2*(r+w+pad+dpml)
   13 cell_size = mp.Vector3(sxy,sxy)
   14 
   15 pml_layers = [mp.PML(dpml)]
   16 
   17 nonpml_vol = mp.Volume(mp.Vector3(), size=mp.Vector3(sxy-2*dpml,sxy-2*dpml))
   18 
   19 geometry = [mp.Cylinder(radius=r+w, material=mp.Medium(index=n)),
   20             mp.Cylinder(radius=r)]
   21 
   22 fcen = 0.118
   23 df = 0.08
   24 
   25 src = [mp.Source(mp.ContinuousSource(fcen,fwidth=df),
   26                  component=mp.Ez,
   27                  center=mp.Vector3(r+0.1)),
   28        mp.Source(mp.ContinuousSource(fcen,fwidth=df),
   29                  component=mp.Ez,
   30                  center=mp.Vector3(-(r+0.1)),
   31                  amplitude=-1)]
   32 
   33 symmetries = [mp.Mirror(mp.X,phase=-1),
   34               mp.Mirror(mp.Y,phase=+1)]
   35 
   36 sim = mp.Simulation(cell_size=cell_size,
   37                     geometry=geometry,
   38                     sources=src,
   39                     resolution=10,
   40                     force_complex_fields=True,
   41                     symmetries=symmetries,
   42                     boundary_layers=pml_layers)
   43 
   44 num_tols = 5
   45 tols = np.power(10, np.arange(-8.0,-8.0-num_tols,-1.0))
   46 ez_dat = np.zeros((122,122,num_tols), dtype=np.complex_)
   47 
   48 for i in range(num_tols):
   49     sim.init_sim()
   50     sim.solve_cw(tols[i], 10000, 10)
   51     ez_dat[:,:,i] = sim.get_array(vol=nonpml_vol, component=mp.Ez)
   52 
   53 err_dat = np.zeros(num_tols-1)
   54 for i in range(num_tols-1):
   55     err_dat[i] = LA.norm(ez_dat[:,:,i]-ez_dat[:,:,num_tols-1])
   56 
   57 plt.figure(dpi=150)
   58 plt.loglog(tols[:num_tols-1], err_dat, 'bo-');
   59 plt.xlabel("frequency-domain solver tolerance");
   60 plt.ylabel("L2 norm of error in fields");
   61 plt.show()
   62 
   63 eps_data = sim.get_array(vol=nonpml_vol, component=mp.Dielectric)
   64 ez_data = np.real(ez_dat[:,:,num_tols-1])
   65 
   66 plt.figure()
   67 plt.imshow(eps_data.transpose(), interpolation='spline36', cmap='binary')
   68 plt.imshow(ez_data.transpose(), interpolation='spline36', cmap='RdBu', alpha=0.9)
   69 plt.axis('off')
   70 plt.show()
   71 
   72 if np.all(np.diff(err_dat) < 0):
   73     print("PASSED solve_cw test: error in the fields is decreasing with increasing resolution")
   74 else:
   75     print("FAILED solve_cw test: error in the fields is NOT decreasing with increasing resolution")
   76 
   77 sim.reset_meep()
   78 
   79 src = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
   80                  component=mp.Ez,
   81                  center=mp.Vector3(r+0.1)),
   82        mp.Source(mp.GaussianSource(fcen,fwidth=df),
   83                  component=mp.Ez,
   84                  center=mp.Vector3(-(r+0.1)),
   85                  amplitude=-1)]
   86 
   87 sim = mp.Simulation(cell_size=mp.Vector3(sxy,sxy),
   88                     geometry=geometry,
   89                     sources=src,
   90                     resolution=10,
   91                     symmetries=symmetries,
   92                     boundary_layers=pml_layers)
   93 
   94 dft_obj = sim.add_dft_fields([mp.Ez], fcen, fcen, 1, where=nonpml_vol)
   95 
   96 sim.run(until_after_sources=100)
   97 
   98 eps_data = sim.get_array(vol=nonpml_vol, component=mp.Dielectric)
   99 ez_data = np.real(sim.get_dft_array(dft_obj, mp.Ez, 0))
  100 
  101 plt.figure()
  102 plt.imshow(eps_data.transpose(), interpolation='spline36', cmap='binary')
  103 plt.imshow(ez_data.transpose(), interpolation='spline36', cmap='RdBu', alpha=0.9)
  104 plt.axis('off')
  105 plt.show()