NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.25k stars 626 forks source link

Prism does weird things when placed with "center=" kwarg #2144

Open sequoiap opened 2 years ago

sequoiap commented 2 years ago

When defining geometries, Prisms don't seem to build properly when offset with the "center" argument. The z-value changes to some weird and whacky value, even when the z-offset is zero. This has the effect of rendering weird shapes for the dielectric that don't match the defined Prism when using plot2D.

Script to reproduce the problem:

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

mp.quiet(quietval=True)
RESOLUTION = 50
OFFSET_Y = 0.5

Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)

WG_W0 = 0.5
WG_W1 = 1.5
WG_RUN_L = 2
WG_TAPER_L = 10

Sx = 2*WG_RUN_L + WG_TAPER_L
Sy = 5
cell_size = (Sx, Sy)

PML_THICKNESS = 1.0
pml_layers = [mp.PML(PML_THICKNESS)]

FCEN = 1/1.55
WIDTH = 0.2
FWIDTH = WIDTH * FCEN
source_center  = mp.Vector3(-Sx / 2 + PML_THICKNESS + 0.5, OFFSET_Y, 0)
source_size    = mp.Vector3(0,1,0)
kpoint = mp.Vector3(1,0,0)
src = mp.GaussianSource(frequency=FCEN, fwidth=FWIDTH)
source = [mp.EigenModeSource(src,
                    eig_band=1,
                    direction=mp.NO_DIRECTION,
                    eig_kpoint=kpoint,
                    size=source_size,
                    center=source_center)]

taper_vtx = [
    mp.Vector3(-WG_TAPER_L/2, WG_W0/2),
    mp.Vector3(WG_TAPER_L/2, WG_W1/2),
    mp.Vector3(WG_TAPER_L/2, -WG_W1/2),
    mp.Vector3(-WG_TAPER_L/2, -WG_W0/2),
]

# Taper (prism) with center specified
geometry_bug = [
    mp.Block(center=mp.Vector3(-(WG_RUN_L/2 + WG_TAPER_L/2), OFFSET_Y), material=Si, size=mp.Vector3(WG_RUN_L, WG_W0)),
    mp.Prism(taper_vtx, mp.inf, center=mp.Vector3(0, OFFSET_Y), material=Si),
    mp.Block(center=mp.Vector3((WG_RUN_L/2 + WG_TAPER_L/2), OFFSET_Y), material=Si, size=mp.Vector3(WG_RUN_L, WG_W1)),
]

# Taper (prism) with predefined vertices shift
geometry_ok = [
    mp.Block(center=mp.Vector3(-(WG_RUN_L/2 + WG_TAPER_L/2), OFFSET_Y), material=Si, size=mp.Vector3(WG_RUN_L, WG_W0)),
    mp.Prism([taper + mp.Vector3(0, OFFSET_Y) for taper in taper_vtx], mp.inf, material=Si),
    mp.Block(center=mp.Vector3((WG_RUN_L/2 + WG_TAPER_L/2), OFFSET_Y), material=Si, size=mp.Vector3(WG_RUN_L, WG_W1)),
]

sim1 = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=geometry_bug,
                    sources=source,
                    default_material=SiO2,
                    resolution=RESOLUTION)

sim2 = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=geometry_ok,
                    sources=source,
                    default_material=SiO2,
                    resolution=RESOLUTION)

if mp.am_really_master():
    fig, axs = plt.subplots(1,2,figsize=(16,6))
    sim1.plot2D(ax=axs[0])
    sim2.plot2D(ax=axs[1])
    plt.savefig("prismbug.png")

    print(geometry_bug[1].vertices)
    # Output:
    # [Vector3<-5.0, 0.75, -5e+19>, Vector3<5.0, 1.25, -5e+19>, Vector3<5.0, -0.25, -5e+19>, Vector3<-5.0, 0.25, -5e+19>]

Result of plot2D (buggy on left, workaround on right):

image

OS Info: Red Hat Enterprise Linux Server Version 7.9 PyMeep Version: 1.23.0 Conda Version: 4.9.2

Yaraslaut commented 2 years ago

Issue arise because center(centroid) of the prism has infinite position along z direction, if you will change mp.inf to some finite number everything will work fine.

This should fix it https://github.com/Yaraslaut/meep/commit/2eabb4fb34ff76816012a5751723118bab8cd738

stevengj commented 2 years ago

Thanks @Yaraslaut, can you submit a PR?