{"action":"closed","issue":{"url":"https://api.github.com/repos/stevengj/meep/issues/597","repository_url":"https://api.github.com/repos/stevengj/meep","labels_url":"https://api.github.com/repos/stevengj/meep/issues/597/labels{/name}","comments_url":"https://api.github.com/repos/stevengj/meep/issues/597/comments","events_url":"https://api.github.com/repos/stevengj/meep/issues/597/events","html_url":"https://github.com/stevengj/meep/issues/597","id":378117740,"node_id":"MDU6SXNzdWUzNzgxMTc3NDA=","number":597,"title":"non-zero kz not passed to get_eigenmode_coefficients for 2d cell","user":{"login":"oskooi","id":7152530,"node_id":"MDQ6VXNlcjcxNTI1MzA=","avatar_url":"https://avatars1.githubusercontent.com/u/7152530?v=4","gravatar_id":"","url":"https://api.github.com/users/oskooi","html_url":"https://github.com/oskooi","followers_url":"https://api.github.com/users/oskooi/followers","following_url":"https://api.github.com/users/oskooi/following{/other_user}","gists_url":"https://api.github.com/users/oskooi/gists{/gist_id}","starred_url":"https://api.github.com/users/oskooi/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/oskooi/subscriptions","organizations_url":"https://api.github.com/users/oskooi/orgs","repos_url":"https://api.github.com/users/oskooi/repos","events_url":"https://api.github.com/users/oskooi/events{/privacy}","received_events_url":"https://api.github.com/users/oskooi/received_events","type":"User","site_admin":false},"labels":[{"id":88929133,"node_id":"MDU6TGFiZWw4ODkyOTEzMw==","url":"https://api.github.com/repos/stevengj/meep/labels/bug","name":"bug","color":"fc2929","default":true}],"state":"closed","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":1,"created_at":"2018-11-07T02:35:55Z","updated_at":"2018-11-10T16:54:40Z","closed_at":"2018-11-10T16:54:40Z","author_association":"CONTRIBUTOR","body":"The k<sub>z</sub> component of the `k_point` is ignored by `get_eigenmode_coefficients` (taken to be as just `0`) for a 2d cell (in *x* and *y*) even though k<sub>z</sub> is non-zero and the simulation itself is 3d.\r\n\r\nAs a demonstration of this bug, the [binary grating example](https://github.com/stevengj/meep/blob/master/python/examples/binary_grating_oblique.py) is modified in the script below to include a non-zero k<sub>z</sub>. The `k_point` had originally been \"in the plane\" with k<sub>z</sub>=0 resulting in a 2d simulation.\r\n\r\n```py\r\nimport meep as mp\r\nimport math\r\nimport cmath\r\nimport numpy as np\r\n\r\nresolution = 40        # pixels/μm                                                                                                            \r\n\r\ndpml = 1.0             # PML thickness                                                                                                        \r\ndsub = 1.0             # substrate thickness                                                                                                  \r\ndpad = 1.0             # length of padding between grating and pml                                                                            \r\ngp = 10.0              # grating period                                                                                                       \r\ngh = 0.5               # grating height                                                                                                       \r\ngdc = 0.5              # grating duty cycle                                                                                                   \r\n\r\nsx = dpml+dsub+gh+dpad+dpml\r\nsy = gp\r\n\r\ncell_size = mp.Vector3(sx,sy,0)\r\n\r\n# replace anisotropic PML with isotropic Absorber to attenuate parallel-directed fields of oblique source                                     \r\nabs_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]\r\n\r\nwvl = 0.5              # center wavelength                                                                                                    \r\nfcen = 1/wvl           # center frequency                                                                                                     \r\ndf = 0.05*fcen         # frequency width                                                                                                      \r\n\r\nng = 1.5\r\nglass = mp.Medium(index=ng)\r\n\r\n# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                \r\ntheta_out = math.radians(20.3)\r\n\r\n# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                               \r\ntheta_in = math.radians(10.7)\r\n\r\n# k (in source medium) with correct length (plane of incidence: XY)                                                                           \r\nk = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)\r\n\r\nsymmetries = []\r\neig_parity = mp.NO_PARITY\r\nif theta_out == 0:\r\n  eig_parity = mp.ODD_Z\r\n  if theta_in == 0:\r\n    k = mp.Vector3(0,0,0)\r\n    symmetries = mp.Mirror(mp.Y)\r\n    eig_parity += mp.EVEN_Y\r\n\r\ndef pw_amp(k,x0):\r\n  def _pw_amp(x):\r\n    return cmath.exp(1j*2*math.pi*k.dot(x+x0))\r\n  return _pw_amp\r\n\r\nsrc_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)\r\nsources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),\r\n                     component=mp.Ez,\r\n                     center=src_pt,\r\n                     size=mp.Vector3(0,sy,0),\r\n                     amp_func=pw_amp(k,src_pt))]\r\n\r\nsim = mp.Simulation(resolution=resolution,\r\n                    cell_size=cell_size,\r\n                    boundary_layers=abs_layers,\r\n                    k_point=k,\r\n                    default_material=glass,\r\n                    sources=sources,\r\n                    symmetries=symmetries)\r\n\r\nrefl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)\r\nrefl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))\r\n\r\nsim.run(until_after_sources=200)\r\n\r\ninput_flux = mp.get_fluxes(refl_flux)\r\ninput_flux_data = sim.get_flux_data(refl_flux)\r\n\r\nsim.reset_meep()\r\n\r\ngeometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),\r\n            mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]\r\n\r\nsim = mp.Simulation(resolution=resolution,\r\n                    cell_size=cell_size,\r\n                    boundary_layers=abs_layers,\r\n                    geometry=geometry,\r\n                    k_point=k,\r\n                    sources=sources,\r\n                    symmetries=symmetries)\r\n\r\nrefl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))\r\nsim.load_minus_flux_data(refl_flux,input_flux_data)\r\n\r\ntran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)\r\ntran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))\r\n\r\nsim.run(until_after_sources=300)\r\n\r\nnm_r = np.floor((math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp) # \\\r\nnumber of reflected orders                                                                                                                    \r\nif theta_in == 0 and theta_out == 0:\r\n  nm_r = nm_r*0.5 # since eig_parity removes degeneracy in y-direction                                                                        \r\nRsum = 0\r\nfor nm in range(int(nm_r)):\r\n  res = sim.get_eigenmode_coefficients(refl_flux, [nm+1], eig_parity=eig_parity)\r\n  r_coeffs = res.alpha\r\n  r_kdom = res.kdom[0]\r\n  Rmode = (abs(r_coeffs[0,0,1])**2)/input_flux[0]\r\n  r_angle = np.sign(r_kdom.y)*math.acos(r_kdom.x/(ng*fcen))\r\n  print(\"refl:, {}, {:.2f}, {:.8f}\".format(nm,math.degrees(r_angle),Rmode))\r\n  Rsum += Rmode\r\n\r\nnm_t = np.floor((math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp) # number\\\r\n of transmitted orders                                                                                                                        \r\nif theta_in == 0 and theta_out == 0:\r\n  nm_t = nm_t*0.5 # since eig_parity removes degeneracy in y-direction                                                                        \r\nTsum = 0\r\nfor nm in range(int(nm_t)):\r\n  res = sim.get_eigenmode_coefficients(tran_flux, [nm+1], eig_parity=eig_parity)\r\n  t_coeffs = res.alpha\r\n  t_kdom = res.kdom[0]\r\n  Tmode = (abs(t_coeffs[0,0,0])**2)/input_flux[0]\r\n  t_angle = np.sign(t_kdom.y)*math.acos(t_kdom.x/fcen)\r\n  print(\"tran:, {}, {:.2f}, {:.8f}\".format(nm,math.degrees(t_angle),Tmode))\r\n  Tsum += Tmode\r\n\r\nprint(\"mode-coeff:, {:.6f}, {:.6f}, {:.6f}\".format(Rsum,Tsum,Rsum+Tsum))\r\n\r\nr_flux = mp.get_fluxes(refl_flux)\r\nt_flux = mp.get_fluxes(tran_flux)\r\nRflux = -r_flux[0]/input_flux[0]\r\nTflux =  t_flux[0]/input_flux[0]\r\nprint(\"poynting-flux:, {:.6f}, {:.6f}, {:.6f}\".format(Rflux,Tflux,Rflux+Tflux))\r\n```\r\nThe output indicates that the cell is 3d which is correct.\r\n```\r\n-----------\r\nInitializing structure...\r\nWorking in 3D dimensions.\r\nComputational cell is 4.5 x 10 x 0.025 with resolution 40\r\n     block, center = (-1.25,0,0)\r\n          size (2,1e+20,1e+20)\r\n          axes (1,0,0), (0,1,0), (0,0,1)\r\n          dielectric constant epsilon diagonal = (2.25,2.25,2.25)\r\n     block, center = (0,0,0)\r\n          size (0.5,5,1e+20)\r\n          axes (1,0,0), (0,1,0), (0,0,1)\r\n          dielectric constant epsilon diagonal = (2.25,2.25,2.25)\r\ntime for set_epsilon = 0.336185 s\r\ntime for set_conductivity = 0.006742 s\r\ntime for set_conductivity = 0.00620389 s\r\ntime for set_conductivity = 0.00629306 s\r\ntime for set_conductivity = 0.00618911 s\r\ntime for set_conductivity = 0.00618219 s\r\ntime for set_conductivity = 0.0062089 s\r\n-----------\r\nMeep: using complex fields.\r\n...\r\n```\r\nNote that the `eig_parity` is set to `NO_PARITY` due to the non-zero k<sub>z</sub>. In this example, k<sub>z</sub> is `1.0408`. All modes computed by MPB via `get_eigenmode_coefficients` should have the same value. However, the k<sub>z</sub> values of these modes are always 0.\r\n```\r\n..\r\nDominant planewave for band 4: (1.998494,-0.077596,0.000000)\r\n...\r\nDominant planewave for band 26: (1.900688,0.622404,0.000000)\r\n..\r\n```\r\nThis indicates that `get_eigenmode_coefficients` is not correctly inferring the dimensions of the cell.\r\n\r\nAs an additional note, calling `get_eigenmode` separately with a non-zero k<sub>z</sub> does yield modes with the same k<sub>z</sub>. This is demonstrated in the following script.\r\n\r\n```py\r\nimport meep as mp\r\nimport math\r\n\r\nresolution = 50        # pixels/μm                                                                                                            \r\n\r\ndpml = 2               # PML thickness                                                                                                        \r\ndmat = 10              # length of material                                                                                                   \r\n\r\nsx = dpml+dmat+dpml\r\nsy = 10\r\n\r\ncell_size = mp.Vector3(sx,sy,0)\r\npml_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]\r\n\r\nwvl = 0.5              # center wavelength                                                                                                    \r\nfcen = 1/wvl           # center frequency                                                                                                     \r\n\r\nng = 1.5\r\nglass = mp.Medium(index=ng)\r\n\r\nnum_band = 6\r\n\r\n# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                \r\ntheta_out = math.radians(20.3)\r\n\r\n# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                               \r\ntheta_in = math.radians(10.7)\r\n\r\n# k (in source medium) with correct length (plane of incidence: XY)                                                                           \r\nk = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)\r\n\r\nsim = mp.Simulation(resolution=resolution,\r\n                    cell_size=cell_size,\r\n                    boundary_layers=pml_layers,\r\n                    k_point=k)\r\n\r\nsim.init_sim()\r\n\r\nmon_pt = mp.Vector3(0.5*sx-dpml-0.2*dmat,0,0)\r\nmon_EigenmodeData = sim.get_eigenmode(fcen, mp.X, mp.Volume(center=mon_pt,size=mp.Vector3(0,sy,0)), num_band, k, verbose=True)\r\n```\r\nThe output shows a `KPOINT` with the correct k<sub>z</sub> value.\r\n```\r\n-----------\r\nInitializing structure...\r\nWorking in 3D dimensions.\r\nComputational cell is 14 x 10 x 0.02 with resolution 50\r\ntime for set_epsilon = 0.914184 s\r\ntime for set_conductivity = 0.0185812 s\r\ntime for set_conductivity = 0.0202799 s\r\ntime for set_conductivity = 0.018465 s\r\ntime for set_conductivity = 0.0202739 s\r\ntime for set_conductivity = 0.0190289 s\r\ntime for set_conductivity = 0.019552 s\r\n-----------\r\nMeep: using complex fields.\r\nKPOINT: 2.76475, 0.522404, 1.04081\r\n    near maximum in trace\r\n...\r\n```"}}
{"actor":{"display_login":"stevengj"}}