From 5f2d84f11569fc1a0cf6ef3c563205d305d450b5 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 27 Feb 2020 18:18:08 -0800 Subject: [PATCH] Python test for subpixel smoothing of prisms (#1135) * Python test for subpixel smoothing of prisms * add test involving marching squares algorithm * whoops, wrong link to conda pack --- doc/docs/FAQ.md | 2 +- doc/docs/Installation.md | 4 + doc/docs/Python_User_Interface.md | 2 +- python/Makefile.am | 1 + python/tests/data/prism_vertices.npz | Bin 0 -> 11496 bytes python/tests/prism.py | 125 +++++++++++++++++++++++++++ 6 files changed, 132 insertions(+), 2 deletions(-) create mode 100644 python/tests/data/prism_vertices.npz create mode 100644 python/tests/prism.py diff --git a/doc/docs/FAQ.md b/doc/docs/FAQ.md index 40eb76e05..c48dce5cd 100644 --- a/doc/docs/FAQ.md +++ b/doc/docs/FAQ.md @@ -109,7 +109,7 @@ Note also that, as a consequence of the above analysis, ε must go to a positive ### Why are there strange peaks in my reflectance/transmittance spectrum when modeling planar or periodic structures? -There are two possible explanations: (1) the simulation run time may be too short and your results have not sufficiently [converged](#checking-convergence), or (2) you may be using a higher-dimensional cell with multiple periods (a supercell) which introduces unwanted additional modes due to band folding. Modeling flat/planar structures typically requires a 1d cell and periodic structures a single unit cell in 2d/3d. For more details, see Section 4.6 ("Sources in Supercells") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). Note that a 1d cell must be along the $z$ direction with only the Ex and Hy field components permitted. +There are two possible explanations: (1) the simulation run time may be too short or the resolution may be too low and thus your results have not sufficiently [converged](#checking-convergence), or (2) you may be using a higher-dimensional cell with multiple periods (a supercell) which introduces unwanted additional modes due to band folding. One indication that band folding is present is that small changes in the resolution produce large and unexpected changes in the frequency spectra. Modeling flat/planar structures typically requires a 1d cell and periodic structures a single unit cell in 2d/3d. For more details, see Section 4.6 ("Sources in Supercells") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). Note that a 1d cell must be along the $z$ direction with only the Ex and Hy field components permitted. ### How do I model the solar radiation spectrum? diff --git a/doc/docs/Installation.md b/doc/docs/Installation.md index 71287fcc9..19d744e26 100644 --- a/doc/docs/Installation.md +++ b/doc/docs/Installation.md @@ -123,6 +123,10 @@ print(mp.__version__) This will show something like `1.11.0-1-g415bc8eb` where the first three digits (`1.11.0`) refer to a stable tarball release, the following digit is the number of commits after this stable release, and the eight characters following the `g` in the final string refer to the commit hash. +### Non-Networked Systems + +To install the PyMeep Conda package on a [non-networked system](https://docs.anaconda.com/anaconda/user-guide/tasks/install-packages/#installing-packages-on-a-non-networked-air-gapped-computer), using the bz2 tarball of the [official release](https://anaconda.org/conda-forge/pymeep/files) or [nightly build](https://anaconda.org/simpetus/pymeep/files) will *not* work without the dependencies. A possible workaround is [Conda-Pack](https://github.com/conda/conda-pack). + Installation on Linux ------------------------- diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index d3bcf0b87..a4e1e3259 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -1852,7 +1852,7 @@ with the following input parameters: + `omega`: optional frequency point over which the average eigenvalue of the dielectric and permeability tensors are evaluated (defaults to 0). -For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional omega parameter (defaults to 0). +For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional `omega` parameter (defaults to 0). **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for simulations involving Bloch-periodic boundaries (via `k_point`) will result in arrays that have slightly *different* dimensions than e.g. `get_array(center=meep.Vector3(), size=cell_size, component=meep.Dielectric`, etc. (i.e., the slice spans the entire cell volume `cell_size`). Neither of these approaches is "wrong", they are just slightly different methods of fetching the boundaries. The key point is that if you pass the same value for the `size` parameter, or use the default, the slicing routines always give you the same-size array for all components. You should *not* try to predict the exact size of these arrays; rather, you should simply rely on Meep's output. diff --git a/python/Makefile.am b/python/Makefile.am index b5edbc501..97a042bdb 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -65,6 +65,7 @@ TESTS = \ $(TEST_DIR)/n2f_periodic.py \ $(TEST_DIR)/oblique_source.py \ $(TEST_DIR)/physical.py \ + $(TEST_DIR)/prism.py \ $(TEST_DIR)/pw_source.py \ $(TEST_DIR)/refl_angular.py \ $(TEST_DIR)/ring.py \ diff --git a/python/tests/data/prism_vertices.npz b/python/tests/data/prism_vertices.npz new file mode 100644 index 0000000000000000000000000000000000000000..79e96e70c1ec772a7a1d7e845b2cac58a45304e2 GIT binary patch literal 11496 zcmb`NPi&l56~(8d&?zU6d{U$#|n@lDi?_bC=!aOBNnK2lL|o!QjmZZ zQjy34SW(3rHW3Ri+KnGUky=nNHA!kaIF2VC{~wQ6ux94F=XXB&4&KnZil%4Y``-87 zz31L@@3W`Qop^krGy3_}AHOqo{WJe~^(!CibVmQ47=HJ{^JhQ(^30Ea^TM~jG4<`o zIzQ~Ze(I&^%ip_n>I;*nzIgHTr_M~Cy7Gr=LH2X7cPaXC{9-`Tzc$JbnI+zx?}jzw}>x>*-(YFBb4`e>eL1tMQ+Qqc6TS zZ+vy)#@D)6W;3olKU!@#-&pGz>wU&2m#)uz>Eu#t|zL(Q`HS)h$}L>OI=e zZjE!Tm2T`)`(G{3gdOWsaSr33H1EAzUL9?B=0NK&mRM&{u`?JugRyhM+Q&ej^Tybj z`Pdm8=Uno9?9BeKGv~q1VA-c~e(WqiRqTxYurv1Mlkb_wbF9O_@h)X@eOJFSv& z+V%S7X5l?pWtIGTq|f@F<^`=vg2}S5J;ZCHjbi4~F+(cz+bndc8P*r<2Aj z7~X?}NAn(s_l(1PaLixue0a}zWV<^x<_+ZUyTnO__t*#CgW)~uE&o*f6YH)S!+X>d zmBQcfK35O$1>S?sIV^Y;l1pGeX-BHZ{iu=r+637)(UpQZrCx~6}#gH z{DPmd-?CrD^8ydx#nCS1-8UaU-TCBjufu;$m-5_)(4`F5_@G@%c=xb8*`%L$QOvvto4n_7R(PXzwL3d z3*+I%n;s_{F&-><$K%QlRDWf?36H<2`WdK;xoxcMVJwyyi$%_k5e5HL_0)D*p(+RN zvcU0q&s|wW>HDg-4oaM3$vC2I+1Wj!ZCO=poBbo|`--Of9#<6kjw7OFRb9@$qTzvY zoPW*g6V=2+bwT@NFL!ydzjbv*ZI>5z*(<0YY8Tm{tFCF6JBvz|hz4@20A~ZTYuXa`mEp%YQw0rDW&e=Z^fffr_|ml^x};wreGO2cLJ{#e(H0R7>{n zmeGy;`8W2|6|B3fCEJI8_GRyN-Qn%az7NN$!SL6AF7&VAH!M@JFD*yHC&ztV~ zU+#({uJ4ymHH0{D*9=B2z<+9P5thcfXvJy0)!fbBlUg3s$!6VqR4MMV!`}nh`;PoS zkUbZT!w>DzsCba#qB zb^U$ubIllCpxcW)wO)y5bb;BMkZ(oEH|TJyqA5@ zt*R>bzQf)-J{RZY+~@*&iJh<;I)L7S(PR98U+`1R>-yFQzvAcYcX$C$;0-*2SMUtp zp$pV2V0a6U;Wa$Zyw7nF`(3HT5pjhs#Qxg2i+EphoD#RhagJ-camMM!$rI!a@<`4r zIq$@C3(iaAE%F$7jXan0Ue1d-Z^rv?I#0*@=8bb+&UriM^_=$uU-Nq5?=jsd>_6+y z|JJ2B&(FFR^XEODbusHEx{5u~UHl_^7O?DZEdLnGf2s=y%13RTAb-1V5I>CRYIVEe zx&(gax`nP7K2!G)54kRi^Q^dzqOKyIsJnV+CD_jsgK)K9%%cRu1Zx8>reLT7l`#31*g5tKI z3t)5sj4mkdWzR*`4+HV@J!9n^W921dbRp&qe9NmmR?r17x&THOz_HG}_;;xF+It*y zf%)hH7+pAu+j}H*A&rw)3%bBKx&V&xInPHI7)KYt=mHpB0HX_FbRpTZNau^}Z!G^9 zqYKQBe6&cviND1g-!mDb3-E=06~3Vhk$(>4-zD95929f`KBEipA6)>W3t)61&e_*_ zR}1lmE`ZSmZhONXYbC|AZ->!^hy&jA@b{+Yg`GE5Uufi6%_z$17C z&)^+A%)E^AugZV@0VC*{Q+@)DE|pq`_T08FOT_yu>4g1l}!88 e$3HQAO243)I``4x8lBFG;eT%qKc77_LjMJT)D8>) literal 0 HcmV?d00001 diff --git a/python/tests/prism.py b/python/tests/prism.py new file mode 100644 index 000000000..82b0edbf5 --- /dev/null +++ b/python/tests/prism.py @@ -0,0 +1,125 @@ +from __future__ import division + +import os +import unittest +import numpy as np +import meep as mp + +class TestPrism(unittest.TestCase): + + def prism_marching_squares(self,npts): + resolution = 50 + + cell = mp.Vector3(3,3) + + data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data')) + vertices_file = os.path.join(data_dir, 'prism_vertices.npz') + vertices_obj = np.load(vertices_file) + + ## prism vertices precomputed for a circle of radius 1.0 using + ## marching squares algorithm of skimage.measure.find_contours + ## ref: https://github.com/NanoComp/meep/issues/1060 + vertices_data = vertices_obj["N{}".format(npts)] + vertices = [mp.Vector3(v[0],v[1],0) for v in vertices_data] + + geometry = [mp.Prism(vertices, + height=mp.inf, + material=mp.Medium(epsilon=12))] + + sim = mp.Simulation(cell_size=cell, + geometry=geometry, + resolution=resolution) + + sim.init_sim() + + prism_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps) + + sim.reset_meep() + + geometry = [mp.Cylinder(radius=1.0, + center=mp.Vector3(), + height=mp.inf, + material=mp.Medium(epsilon=12))] + + sim = mp.Simulation(cell_size=cell, + geometry=geometry, + resolution=resolution) + + sim.init_sim() + + cyl_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps) + + print("epsilon-sum:, {} (prism-msq), {} (cylinder), {} (relative error)".format(abs(prism_eps),abs(cyl_eps),abs((prism_eps-cyl_eps)/cyl_eps))) + + return abs((prism_eps-cyl_eps)/cyl_eps) + + + def prism_circle(self,npts,r): + resolution = 50 + + cell = mp.Vector3(3,3) + + ### prism vertices computed as npts equally-spaced points + ### along the circumference of a circle with radius r + angles = 2*np.pi/npts * np.arange(npts) + vertices = [mp.Vector3(r*np.cos(ang),r*np.sin(ang)) for ang in angles] + geometry = [mp.Prism(vertices, + height=mp.inf, + material=mp.Medium(epsilon=12))] + + sim = mp.Simulation(cell_size=cell, + geometry=geometry, + resolution=resolution) + + sim.init_sim() + + prism_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps) + + sim.reset_meep() + + geometry = [mp.Cylinder(radius=r, + center=mp.Vector3(), + height=mp.inf, + material=mp.Medium(epsilon=12))] + + sim = mp.Simulation(cell_size=cell, + geometry=geometry, + resolution=resolution) + + sim.init_sim() + + cyl_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps) + + print("epsilon-sum:, {} (prism-cyl), {} (cylinder), {} (relative error)".format(abs(prism_eps),abs(cyl_eps),abs((prism_eps-cyl_eps)/cyl_eps))) + + return abs((prism_eps-cyl_eps)/cyl_eps) + + + def test_prism(self): + print("Testing Prism object using marching squares algorithm...") + d = [self.prism_marching_squares(92), + self.prism_marching_squares(192), + self.prism_marching_squares(392)] + + self.assertLess(d[1],d[0]) + self.assertLess(d[2],d[1]) + + print("Testing Prism object using circle formula...") + r = 1.0458710786934182 + d = [self.prism_circle(51,r), + self.prism_circle(101,r), + self.prism_circle(201,r)] + + self.assertLess(d[1],d[0]) + self.assertLess(d[2],d[1]) + + r = 1.2896871096581341 + d = [self.prism_circle(31,r), + self.prism_circle(61,r), + self.prism_circle(121,r)] + + self.assertLess(d[1],d[0]) + self.assertLess(d[2],d[1]) + +if __name__ == '__main__': + unittest.main()