[Meep-discuss] Comparing Poynting-Flux from get_fluxes with time-average

[Karl-Johan Olofsson]

Hello!


I am comparing the flux of the whole simulation cell using meep.get_fluxes with the far-field flux out from a sphere with a radius much larger than the length of the sim. cell. However, I get a complete mismatch of the flux values where I believe they should be close to same and I can't figure out the error.




Here is a scaled down version of the code:


sx=sy=sz=4

dpml=1.6

symmetry=[mp.Mirror(mp.X),mp.Mirror(mp.Y,phase=-1)]

pml_layer=[mp.PML(dpml)]

padding=0.1

fcen=2
df=1.2
nfreq=4



source2=[mp.Source(mp.ContinuousSource(frequency=fcen,start_time=0,end_time=1),
component=mp.Ey,
center=mp.Vector3(0,0,0))]

sim=mp.Simulation(cell_size=
mp.Vector3(sx+2*dpml,sy+2*dpml,sz+2*dpml),
symmetries=symmetry,
sources=source2,
dimensions=3,
boundary_layers=pml_layer,
resolution=10)



fluxregion1=mp.FluxRegion( #region -y to calculate flux from
center=mp.Vector3(0,-sy/2+padding,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y,
weight=-1)

fluxregion2=mp.FluxRegion( #region y to calculate flux from
center=mp.Vector3(0,sy/2-padding,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y)

fluxregion3=mp.FluxRegion( # region -x to calculate flux from
center=mp.Vector3(-sx/2+padding,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X,
weight=-1)

fluxregion4=mp.FluxRegion( #region x to calculate flux from
center=mp.Vector3(sx/2-padding,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X)

fluxregion5=mp.FluxRegion( #z-bottom region to calculate flux from
center=mp.Vector3(0,0,sz/2-padding),
size=mp.Vector3(sx-padding*2,sz-padding*2,0),
direction=mp.Z)

fluxregion6=mp.FluxRegion( #z-top region to calculate flux from
center=mp.Vector3(0,0,-sz/2+padding),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z,
weight=-1)

flux_total=sim.add_flux(fcen,df,nfreq,
fluxregion1,fluxregion2,fluxregion3,fluxregion4,fluxregion5,fluxregion6)
--------------
nearfieldregion1=mp.Near2FarRegion( #nearfield -z. above pyramid.
center=mp.Vector3(0,-sz/2+padding,0),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z,
weight=-1)

nearfieldregion6=mp.Near2FarRegion( #nearfield z. below pyramid.
center=mp.Vector3(0,sz/2+padding,0),
size=mp.Vector3(sx-padding*2,sy-padding*2,0),
direction=mp.Z)

nearfieldregion2=mp.Near2FarRegion(
center=mp.Vector3(sx/2,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X)

nearfieldregion3=mp.Near2FarRegion(
center=mp.Vector3(-sx/2,0,0),
size=mp.Vector3(0,sy-padding*2,sz-padding*2),
direction=mp.X,
weight=-1)

nearfieldregion4=mp.Near2FarRegion(
center=mp.Vector3(0,sy/2,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y)

nearfieldregion5=mp.Near2FarRegion(
center=mp.Vector3(0,-sy/2,0),
size=mp.Vector3(sx-padding*2,0,sz-padding*2),
direction=mp.Y,
weight=-1)

nearfield=sim.add_near2far(fcen,df,nfreq,nearfieldregion1,nearfieldregion2,nearfieldregion3,
nearfieldregion4,nearfieldregion5,nearfieldregion6)

-------------

sim.run(until=10)

P_tot_ff = np.zeros(nfreq)
r=100
npts=10
Px=0
Py=0
Pz=0
for n in range(npts):
angleN=(n/npts)*math.pi*2
for m in range(npts):
angleM=(m/npts)*math.pi
ff=sim.get_farfield(nearfield, mp.Vector3(r*math.sin(angleM)*math.cos(angleN),r*math.sin(angleM)*math.sin(angleN),-r*math.cos(angleM)))

i=0
for k in range(nfreq):
Px=(ff[i+1]*np.conjugate(ff[i+5])-ff[i+2]*np.conjugate(ff[i+4]))
Py=(ff[i+2]*np.conjugate(ff[i+3])-ff[i]*np.conjugate(ff[i+5]))
Pz=(ff[i]*np.conjugate(ff[i+4])-ff[i+1]*np.conjugate(ff[i+3]))
Pr=(Px)*math.cos(angleN)*math.sin(angleM)+(Py)*math.sin(angleN)*math.sin(angleM)-(Pz)*math.cos(angleM)
surface_Element=math.pow(r,2)*math.pow(angle,2)*math.pow((1/npts),2)
P_tot_ff[k] += surface_Element*(1/2)*np.real(Pr)
print(n,m,k,P_tot_ff[k])
i=i+6

flux_tot_out = mp.get_fluxes(flux_total)

P_tot_ff= [0.00021628 0.00140266 0.00369171 0.00201077]
flux_tot_out = [0.031127320105939483, 0.19604245927203143, 0.3616725444974181, 0.20233403694380303]

P_tot_ff should be the radial component of the flux summed over the whole sphere. So in my mind, P_tot_ff and flux_tot_out should be very close to the same, but they obviously are not?

I would be very grateful for some help.

Warm Greetings,
Karl-Johan

[Ardavan Oskooi]

Post by Karl-Johan Olofsson
I am comparing the flux of the whole simulation cell using
meep.get_fluxes with the far-field flux out from a sphere with a
radius much larger than the length of the sim. cell. However, I get a
complete mismatch of the flux values where I believe they should be
close to same and I can't figure out the error.

Note that the far fields are computed from the near fields using an
analytic transformation involving the Green's function (i.e., the far
field at a single point in space will consist of contributions from all
the related near fields weighted by the equivalent surface current at
that point; for details, see Section 4.2.1 of our book chapter:
https://arxiv.org/abs/1301.5366). Thus, there is not a simple scaling
relationship between the near and far fields which is why their flux
values are different.

[Steven G. Johnson]

Post by Karl-Johan Olofsson
I am comparing the flux of the whole simulation cell using meep.get_fluxes with the far-field flux out from a sphere with a radius much larger than the length of the sim. cell. However, I get a complete mismatch of the flux values where I believe they should be close to same and I can't figure out the error.

Note that the far fields are computed from the near fields using an analytic transformation involving the Green's function (i.e., the far field at a single point in space will consist of contributions from all the related near fields weighted by the equivalent surface current at that point; for details, see Section 4.2.1 of our book chapter: https://arxiv.org/abs/1301.5366 <https://arxiv.org/abs/1301.5366>). Thus, there is not a simple scaling relationship between the near and far fields which is why their flux values are different.

By Poynting's theorem, the flux spectrum (integrated around a closed surface) should be the same whether you calculate it from the near or far fields (unless there are sources or absorbers in between), with slight differences due to discretization errors. If the two values are not converging to one another as you increase resolution, my first guess would be that you have a bug in your code calculating the flux values.

If you can't figure it out, you might file an issue on github with a minimal example (e.g. a small 2d simulation).