S.GetEpsilon varies from run to run

[sbyrnes321]

When I run S.GetEpsilon, often it is a nice fourier-series approximation of my actual epsilon as expected: (Blue is input epsilon, green is from S.GetEpsilon.)

...but more often it is weirdly scaled and shifted:

It randomly varies from run to run!!! Isn't that weird?

I compiled from source, but this happens whether or not I use FFTW, or BLAS, or LAPACK, or leave PTHREAD_LIB blank. So I don't think it's a threading-related error, although I'm no expert.

I haven't tried this in the lua frontend, only python.

Code:

from __future__ import division, print_function

import matplotlib.pyplot as plt
import numpy as np

import S4

period = 1.
s = 0.6*period #grating bar width
a = 0.4*period
eBar=10.234
thickness = 0.63*period #HCG thickness
lam_vac = 1.3*period #incident wavelength

S = S4.New(Lattice=period, NumBasis=100)

S.SetMaterial(Name = 'mymaterial', Epsilon = eBar)
S.SetMaterial(Name = 'Vacuum', Epsilon = 1)

S.AddLayer(Name='AirAbove', Thickness=0, Material='Vacuum')
S.AddLayer(Name='grating', Thickness=thickness, Material='Vacuum')
S.AddLayerCopy(Name='AirBelow', Thickness=0, Layer='AirAbove')

S.SetRegionRectangle(Layer='grating', Material='mymaterial', Center=(0,0), Angle=0, Halfwidths=(s/2,1))

S.SetExcitationPlanewave(IncidenceAngles=(0,0), sAmplitude = 1,
                         pAmplitude = 0, Order = 0)

S.SetFrequency(lam_vac)

xlist = np.linspace(-period/2, period/2, num=200)
zlist = np.linspace(-0.5*thickness, 1.5*thickness, num=3)
plt.figure()
plt.plot(xlist, [eBar if abs(x)<s/2 else 1 for x in xlist])
plt.plot(xlist, [S.GetEpsilon(x,0,thickness/2).real for x in xlist])

and the makefile is

BLAS_LIB = 
LAPACK_LIB = 
LUA_INC = -I/usr/include/lua5.2
LUA_LIB = -llua5.2 -ldl -lm
FFTW3_INC =
FFTW3_LIB = 
PTHREAD_INC = 
PTHREAD_LIB = 
CHOLMOD_INC = 
CHOLMOD_LIB = 
MPI_INC =
MPI_LIB =
CXX = g++
CC  = gcc
CFLAGS += -O3 -fPIC

but again, I get the same error with other makefiles I've tried.

[victorliu]

Ok, that is totally strange behavior... You might try defining S4_TRACE and DUMP_MATRICES etc in the compilation flags to get more detailed tracing info.

On Wed, Sep 24, 2014 at 9:29 AM, Steven Byrnes [email protected] wrote:

When I run S.GetEpsilon, often it is a nice fourier-series approximation of my actual epsilon as expected: (Blue is input epsilon, green is from S.GetEpsilon.)

[image: example that works] https://cloud.githubusercontent.com/assets/1857674/4391583/8aba32f8-4406-11e4-9cf8-8029492ce3b8.png

...but more often it is weirdly scaled and shifted:

[image: example that did not work] https://cloud.githubusercontent.com/assets/1857674/4391586/8eb6199e-4406-11e4-9ac1-f3ade316de22.png

It randomly varies from run to run!!! Isn't that weird?

I compiled from source, but this happens whether or not I use FFTW, or BLAS, or LAPACK, or leave PTHREAD_LIB blank. So I don't think it's a threading-related error, although I'm no expert.

I haven't tried this in the lua frontend, only python.

Code:

from future import division, print_function

import matplotlib.pyplot as plt import numpy as np

import S4

period = 1. s = 0.6_period #grating bar width a = 0.4_period eBar=10.234 thickness = 0.63_period #HCG thickness lam_vac = 1.3_period #incident wavelength

S = S4.New(Lattice=period, NumBasis=100)

S.SetMaterial(Name = 'mymaterial', Epsilon = eBar) S.SetMaterial(Name = 'Vacuum', Epsilon = 1)

S.AddLayer(Name='AirAbove', Thickness=0, Material='Vacuum') S.AddLayer(Name='grating', Thickness=thickness, Material='Vacuum') S.AddLayerCopy(Name='AirBelow', Thickness=0, Layer='AirAbove')

S.SetRegionRectangle(Layer='grating', Material='mymaterial', Center=(0,0), Angle=0, Halfwidths=(s/2,1))

S.SetExcitationPlanewave(IncidenceAngles=(0,0), sAmplitude = 1, pAmplitude = 0, Order = 0)

S.SetFrequency(lam_vac)

xlist = np.linspace(-period/2, period/2, num=200) zlist = np.linspace(-0.5_thickness, 1.5_thickness, num=3) plt.figure() plt.plot(xlist, [eBar if abs(x)<s/2 else 1 for x in xlist]) plt.plot(xlist, [S.GetEpsilon(x,0,thickness/2).real for x in xlist])

and the makefile is

BLAS_LIB = LAPACK_LIB = LUA_INC = -I/usr/include/lua5.2 LUA_LIB = -llua5.2 -ldl -lm FFTW3_INC = FFTW3_LIB = PTHREAD_INC = PTHREAD_LIB = CHOLMOD_INC = CHOLMOD_LIB = MPI_INC = MPI_LIB = CXX = g++ CC = gcc CFLAGS += -O3 -fPIC

but again, I get the same error with other makefiles I've tried.

— Reply to this email directly or view it on GitHub https://github.com/victorliu/S4/issues/12.

[sbyrnes321]

Probably it's from reading unitialized memory? But I can't figure it out so far.

BTW, I can reproduce this in the lua frontend.

[victorliu]

Can you please post the Lua test code?

[sbyrnes321]

AFAICT the Windows binary download works perfectly and has no issues. Obviously I'm building it wrong somehow. I give up, I'm using the Windows binary from now on.

For what it's worth, here is the lua test code. It works in the windows binary, but in the version I built, it prints six numbers that vary from run to run, and then the last line GetPoyntingFlux gives a segfault for good measure. :-P

S = S4.NewSimulation()
S:SetLattice({1,0}, {0,0}) -- 1D lattice
S:SetNumG(28)

S:AddMaterial("Silicon", {10.234,0}) -- real and imag parts
S:AddMaterial("Vacuum", {1,0})

S:AddLayer('AirAbove', 0, 'Vacuum')
S:AddLayer('Slab', 0.63, 'Vacuum')
S:SetLayerPatternRectangle('Slab', 'Silicon', {0,0}, 0, {0.3, 0})
S:AddLayerCopy('AirBelow', 0, 'AirAbove')

S:SetExcitationPlanewave({0,0}, {0,0}, {1,0})
S:UsePolarizationDecomposition()

freq = 1.3
S:SetFrequency(freq)
print(S:GetEpsilon({0,0,0.3}))
print(S:GetEpsilon({0.1,0,0.3}))
print(S:GetEpsilon({0.2,0,0.3}))
print(S:GetEpsilon({0.3,0,0.3}))
print(S:GetEpsilon({0.4,0,0.3}))
print(S:GetEpsilon({0.5,0,0.3}))

forw_r, back_r, forw_i, back_i = S:GetPoyntingFlux('AirAbove',0)

[olfduh]

I also get some interesting behaviour for epsilon depending on the lattice size. In my understanding of the FMM, there are only some problems when the period is a multiple of the wavelengh. For a simple 1D lattice I get good results if I set the period to 1, but something completely odd when I change it to 0.9 or to 2.3:

I already posted that on NanoHub, not sure where I can contact you guys best :)

Is it possible, that the period for the Fourier expansion of the permittivity is always 1, regardless of what I define in the Python script? If I increase the NumBasis by a factor of 2it has the same effect as decreasing the period by a factor of 2. Multiplying the period and NumBasis by the same number results in the exact same permittivity profile.

Here is the python script I used:

from __future__ import division, print_function

import matplotlib.pyplot as plt
import numpy as np

import S4

period = 2.3
s = 0.5*period #grating bar width
thickness = .5
freq_vac = 1/.633 #incident wavelength

S = S4.New(Lattice=((period,0),(0,0)), NumBasis=200)

S.SetMaterial(Name = 'mymaterial', Epsilon = 2**2)
S.SetMaterial(Name = 'Vacuum', Epsilon = 1)

S.AddLayer(Name='AirAbove', Thickness=0, Material='Vacuum')
S.AddLayer(Name='grating', Thickness=thickness, Material='mymaterial')
S.AddLayerCopy(Name='AirBelow', Thickness=0, Layer='AirAbove')

S.SetRegionRectangle(Layer='grating', Material='Vacuum', Center=(0,0), Angle=0, Halfwidths=(s/2,1))

S.SetExcitationPlanewave(IncidenceAngles=(0,0), sAmplitude = 1,
                         pAmplitude = 0, Order = 0)

S.SetFrequency(freq_vac)

xlist = np.linspace(0, period, num=200)
plt.figure()
plt.plot(xlist, [S.GetEpsilon(x,0,thickness/2).real for x in xlist])