artifacts from grid discretization

[oskooi]

There are unexpectedly large discretization artifacts when using a amplitude function for the Source object. This effect can be demonstrated using an example based on "Method 3" of Tutorial/Custom Source/Stochastic Dipole Emission in Light Emitting Diodes involving a cosine Fourier series (m=7) source amplitude function. The radiated flux from the LED structure at the pulse center frequency shows large variations due to tiny changes in the resolution.

This effect may be due to the way the source amplitude function is interpolated onto the Yee grid.

import numpy as np
import meep as mp
from meep.materials import Ag

resolution = 50

dpml = 1.0
dair = 0.9
hrod = 0.6
wrod = 0.8
dsub = 5.4
dAg = 0.4

sx = 1.5
sy = dpml+dair+hrod+dsub+dAg

cell_size = mp.Vector3(sx,sy)

pml_layers = [mp.PML(direction=mp.Y,
                     thickness=dpml,
                     side=mp.High)]

fcen = 1.0
df = 0.2
ndipole = int(sx*resolution)
run_time = 200

geometry = [mp.Block(material=mp.Medium(index=3.45),
                     center=mp.Vector3(0,0.5*sy-dpml-dair-hrod-0.5*dsub),
                     size=mp.Vector3(mp.inf,dsub,mp.inf)),
            mp.Block(material=Ag,
                     center=mp.Vector3(0,-0.5*sy+0.5*dAg),
                     size=mp.Vector3(mp.inf,dAg,mp.inf)),
            mp.Block(material=mp.Medium(index=3.45),
                     center=mp.Vector3(0,0.5*sy-dpml-dair-0.5*hrod),
                     size=mp.Vector3(wrod,hrod,mp.inf))]

def src_amp_func(n):
    def _src_amp_func(p):
        if n == 0:
            return 1/np.sqrt(sx)
        else:
            return np.sqrt(2/sx) * np.cos(2*n*np.pi*(p.x+0.5*sx)/sx)
    return _src_amp_func

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                     component=mp.Ez,
                     center=mp.Vector3(0,-0.5*sy+dAg+0.5*dsub),
                     size=mp.Vector3(sx,0),
                     amp_func=src_amp_func(7))]

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    k_point=mp.Vector3(),
                    boundary_layers=pml_layers,
                    geometry=geometry,
                    sources=sources)

flux_mon = sim.add_flux(fcen, 0, 1,
                        mp.FluxRegion(center=mp.Vector3(0,0.5*sy-dpml),size=mp.Vector3(sx)))

sim.run(until=run_time)

flux = mp.get_fluxes(flux_mon)
print("flux:, {}".format(flux[0]))

resolution = 50

flux:, 2.8759697666698216e-18

resolution = 51

flux:, 8.409680174807422e-05

resolution = 52

flux:, 3.1591784346399953e-18

resolution = 53

flux:, 7.981497205970976e-05

resolution = 54

flux:, 1.657503397760138e-18

resolution = 55

flux:, 0.00010086263373814502

Additional notes: this effect does not depend on the value of m or the presence of the Ag back reflector.

[oskooi]

I think I have tracked down the cause of this bug: the source amplitude function (amp_func) returns incorrect (?) values at the cell edge whenever a k_point is specified (i.e., periodic boundary conditions). This can be demonstrated by plotting the source amplitude function via the get_source routine over the entire line on which it is defined (including the "ghost" pixels at the cell edge). At the left cell edge (x=0), when k_point is 0 the amplitude function is not sinusoidal but rather discontinuous (region marked by the dotted red circle in the first figure below). By comparison, the amplitude function is sinusoidal when k_point is not specified. The source amplitude function is the same for both cases everywhere else away from the cell edge.

1. k_point=mp.Vector3()

2. no k_point

[stevengj]

I think the problem is that the cosine sources are not periodic (for odd n). Arguably, in this case you should be using the sine and cosine sources (with twice the frequency, so that they are periodic). Not the problem since your frequency is 2πnx/sx.

[stevengj]

Probably you are just seeing discretization slightly breaking the mirror symmetry depending on whether you have an odd or even number of pixels, so that the non-radiative source can radiate due to discretization errors in one case (1e-5 power) and due only to roundoff errors in the other case (1e-17 power).

(This is related to the fact that, if you specify a mirror symmetry plane in x, Meep will actually pad the cell by 1 pixel in some cases to maintain the discrete symmetry.)

[stevengj]

(I think there may be some separate problem in get_source_slice, because an actual glitch that big in the current seems like it would cause a much bigger error than 1e-5 in the power.)

[oskooi]

It is indeed discretization effects related to the number of pixels for the line source (equal to sx, the cell width) which are causing large changes (~13 orders of magnitude) in the radiated flux. In this particular example, an even resolution always results in an x cell dimension of 1.5 but something slightly larger whenever the resolution is odd.

Computational cell is 1.5 x 8.3 x 0 with resolution 50
flux:, 3.8039371409980637e-19
...
Computational cell is 1.5098 x 8.29412 x 0 with resolution 51
flux:, 1.9238223363157157e-06
...
Computational cell is 1.5 x 8.30769 x 0 with resolution 52
flux:, 3.2364410466606466e-19
...
Computational cell is 1.50943 x 8.30189 x 0 with resolution 53
flux:, 1.9478907243015177e-06
...
Computational cell is 1.5 x 8.3 x 0 with resolution 60
flux:, 4.695093780565119e-19
...
Computational cell is 1.50769 x 8.30769 x 0 with resolution 65
flux:, 1.2741839335654913e-06

(Note: The same thing happens when periodic boundaries are not used. Also, no mirror symmetry plane is used in this example and thus the cell is never padded.)

Since this effect is a feature of the grid discretization and not a bug (and not something that we can provide a built-in protection against), it would be useful to document it in Features/Exploiting Symmetry so that users are aware of it.

[oskooi]

Replaced by #1291.