(I sent this message two days ago but it was too big and was not distributed.
Here I cut it off a little bit and send again)
Hi all,
I'm learning MEEP and here is a very simple problem which gave me a surprise.
Considering a structure composed of 3 planar layers of dielectrics, say,
epsilon=3, 2, 3, respectively, and we want to study the plane wave
transmission/reflection at normal incidence. The layers are inifinite in y-z
plane thus the problem is basically a 1-D problem. Of course we can simulate
it in MEEP using the 1D feature. We can also construct an equivalent 2D
problem by considering a planar waveguide of certain height (say, sy) composed
of two perfect electric conductors placed at +/-sy, and place the 3 planar
dielectric layers inside the waveguide. A plane wave can be launched and
propogate in this waveguide. We can also construct a 3D equivalent problem, by
considering a waveguide with PEC walls at y dirrection and PMC walls at z
direction. Again a planar wave with electric field along y direction should be
able to be launched. In each simulation we use a plane, gaussian source along
y direction, and requires the Ey field at the end of the waveguide decays to
1.0e-6.
The 1D and 2D simulation using MEEP were sucessful, finished very fast (the 2D
case finished at time=2276.5) and gave identical results. However, the 3D
simulation had weird behavior. First, the field at the output end of the
waveguide increased extremely slow. At time=3794 it was still 1.1e-21. Then
the magnitude swings back and forth, and finally it begings to increase
drastically, without a sign of converge. The simulation was finally terminated
manually at time=702429, when the maximum of the field was something like 63036.
The 3 situations are theoretically equivallent to each other, but why the 3D
simulation does not converge? Am I missing anything? I have been reading the
.ctl file very carefully but I did not figure out any error yet. Hope someone
can give me some suggestions.
The simulation was done on a cluster and multiprocess were used for all the
cases. The .ctl file of the 2D and 3D simulation is below. The outupt file
for the 2D case and part of the output file for 3D case are also attached, in
case some one would like to dig into it.
Thanks for any suggestions.
Jingjing Li
2D case .ctl file:
------------
(define-param thick1 5); Thickness of the first layerr (define-param thick2 5);
Thickness of the second last layer (define-param thick_diel 50); Thickness of
the sandwitched dielectric layer (define-param thick_pml 200); Thickness of the
PML (define-param thick_stack 200); Thickness of the stack between the PML and
; the layers
(define-param sx (+ thick_pml
(+ thick_stack
(+ thick1
(+ thick_diel
(+ thick2
(+ thick_stack thick_pml)
)))))
); The size along x-direction. Include the PML, the distance between the
; the PML and the layers, etc.
(define-param sy (* thick1 5)); Size along y direction
(set! geometry-lattice (make lattice (size sx sy no-size)));
(set! geometry (list (make block
(center (- (- 0 (/ thick_diel 2)) (/ thick1 2)) 0 0)
(size thick1 sy infinity)
(material (make dielectric (epsilon 3))))
(make block
(center 0 0 0)
(size thick_diel sy infinity)
(material (make dielectric (epsilon 2.0))))
(make block
(center (+ (/ thick2 2) (/ thick_diel 2)) 0 0)
(size thick2 sy infinity)
(material (make dielectric (epsilon 3))))
));
(set! sources (list (make source
(src (make gaussian-src (frequency 0.01) (fwidth
0.005)))
(component Ey)
(center (+ 100 (- 0 (+ (+ (/ thick_diel 2) (/ thick1
2)) thick_stack))) 0 0)
(size 0 sy infinity)
)));
(set! pml-layers (list (make pml (thickness thick_pml) (direction X)))); (set!
resolution 1); (define trans (add-flux 0.01 0.005 50
(make flux-region (center (- (* 0.5 sx) thick_pml 10) 0 0) (size 0
sy))));
(meep-fields-set-boundary fields High Y Metallic); (meep-fields-set-boundary
fields Low Y Metallic);
(run-sources+
(stop-when-fields-decayed 200 Ey
(vector3 (- (* 0.5 sx) thick_pml 10) 0 0)
1e-6));
(display-fluxes trans);
---------
2D out file:
....
-----------
Initializing structure...
Working in 2D dimensions.
block, center = (-27.5,0,0)
size (5,25,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 3
block, center = (0,0,0)
size (50,25,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 2
block, center = (27.5,0,0)
size (5,25,1e+20)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 3
time for set_epsilon = 0.057966 s
-----------
field decay(t = 200.5): 2.20104730041374e-307 / 2.20104730041374e-307 = 1.0
field decay(t = 401.0): 7.01525638552965e-12 / 7.01525638552965e-12 = 1.0
field decay(t = 601.5): 3.55627941401709e-8 / 3.55627941401709e-8 = 1.0
field decay(t = 802.0): 3.31778269088121e-5 / 3.31778269088121e-5 = 1.0
field decay(t = 1002.5): 0.00439031330545536 / 0.00439031330545536 = 1.0
field decay(t = 1203.0): 0.0830320548934512 / 0.0830320548934512 = 1.0
field decay(t = 1403.5): 0.227387342157038 / 0.227387342157038 = 1.0
field decay(t = 1604.0): 0.216835560577182 / 0.227387342157038 =
0.953595563060986
field decay(t = 1804.5): 0.0553931896100671 / 0.227387342157038 =
0.243607181844852
field decay(t = 2005.0): 0.00232470418754533 / 0.227387342157038 =
0.0102235426365107
field decay(t = 2205.5): 1.98672676114093e-5 / 0.227387342157038 =
8.73719153535315e-5
on time step 4553 (time=2276.5), 0.000878568 s/step
field decay(t = 2406.0): 5.23174485458109e-8 / 0.227387342157038 =
2.30080742619699e-7
run 0 finished at t = 2406.0 (4812 timesteps)
flux1:, 0.0075, 1.45138920040302
flux1:, 0.00760204081632653, 3.26552555959972
flux1:, 0.0123979591836734, 8.66617482613586
flux1:, 0.0125, 3.94358407536073
...
...
...
Elapsed run time = 4.47439 s
mpirun (mpirm_setprocs): job done, sending "done" to API
-----------
3D case .ctl file:
---------------
(define-param thick1 5); Thickness of the first layer (define-param thick2 5);
Thickness of the last layer (define-param thick_diel 5); Thickness of the
sandwitched dielectrics (define-param thick_pml 100); Thickness of the PML
(define-param thick_stack 100); Thickness of the stack between the PML and
; the layers (define-param sx (+ thick_pml
(+ thick_stack
(+ thick1
(+ thick_diel
(+ thick2
(+ thick_stack thick_pml)
)))))
); The size along x-direction. Include the PML, the distance between the
; the PML and the layer, etc.
(define-param sy (* thick1 5));
(define-param sz (* thick1 5)); The size along x- and y- direction. 5
; times the thickness of slab
(set! geometry-lattice (make lattice (size sx sy sz))); (set! geometry (list
(make block
(center (- (- 0 (/ thick_diel 2)) (/ thick1 2)) 0 0)
(size thick1 sy sz)
(material (make dielectric (epsilon 3))))
(make block
(center 0 0 0)
(size thick_diel sy sz)
(material (make dielectric (epsilon 2.0))))
(make block
(center (+ (/ thick2 2) (/ thick_diel 2)) 0 0)
(size thick2 sy sz)
(material (make dielectric (epsilon 3))))
));
(set! sources (list (make source
(src (make gaussian-src (frequency 0.01) (fwidth
0.005)))
(component Ey)
(center (- 0 (+ (+ (/ thick_diel 2) (/ thick1 2))
thick_stack) 100) 0 0)
(size 0 sy sz)
)));
(set! pml-layers (list (make pml (thickness thick_pml) (direction X)))); (set!
resolution 1); (define trans (add-flux 0.01 0.005 25
(make flux-region (center (- (* 0.5 sx) thick_pml 10) 0 0) (size 0
sy sz))));
(meep-fields-set-boundary fields High Y Metallic); (meep-fields-set-boundary
fields Low Y Metallic); (meep-fields-set-boundary fields High Z Magnetic);
(meep-fields-set-boundary fields Low Z Magnetic);
(run-sources+
(stop-when-fields-decayed 200 Ey
(vector3 (- (* 0.5 sx) thick_pml 10) 0 0)
1e-6));
(display-fluxes trans);
--------
3D out file:
-------
....
-----------
Initializing structure...
Working in 3D dimensions.
block, center = (-5,0,0)
size (5,25,25)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 3
block, center = (0,0,0)
size (5,25,25)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 2
block, center = (5,0,0)
size (5,25,25)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon = 3
time for set_epsilon = 0.111727 s
-----------
on time step 102 (time=51), 0.0394899 s/step
on time step 223 (time=111.5), 0.0332617 s/step
on time step 343 (time=171.5), 0.0335118 s/step
field decay(t = 200.5): 7.86823638077021e-190 / 7.86823638077021e-190 = 1.0
on time step 462 (time=231), 0.0336207 s/step
on time step 583 (time=291.5), 0.0332875 s/step
on time step 700 (time=350), 0.0343074 s/step
field decay(t = 401.0): 2.17259080821291e-34 / 2.17259080821291e-34 = 1.0
on time step 820 (time=410), 0.0335085 s/step
on time step 941 (time=470.5), 0.0332634 s/step
on time step 1062 (time=531), 0.0333275 s/step
on time step 1183 (time=591.5), 0.0333089 s/step
field decay(t = 601.5): 2.68147525745871e-29 / 2.68147525745871e-29 = 1.0
on time step 1304 (time=652), 0.0332227 s/step
on time step 1425 (time=712.5), 0.033179 s/step
on time step 1545 (time=772.5), 0.0334744 s/step
field decay(t = 802.0): 3.35560267132528e-26 / 3.35560267132528e-26 = 1.0
on time step 1665 (time=832.5), 0.0333847 s/step
on time step 1786 (time=893), 0.0331464 s/step
on time step 1907 (time=953.5), 0.0332618 s/step
field decay(t = 1002.5): 7.38529083253747e-25 / 7.38529083253747e-25 = 1.0
on time step 2028 (time=1014), 0.0331695 s/step
on time step 2149 (time=1074.5), 0.0331732 s/step
on time step 2270 (time=1135), 0.0331627 s/step
on time step 2390 (time=1195), 0.0333693 s/step
on time step 7104 (time=3552), 0.0330666 s/step
field decay(t = 3609.0): 9.95970929690159e-22 / 9.95970929690159e-22 = 1.0
on time step 7225 (time=3612.5), 0.0330928 s/step
on time step 7346 (time=3673), 0.0331924 s/step
on time step 7467 (time=3733.5), 0.0331271 s/step
on time step 7588 (time=3794), 0.033172 s/step
field decay(t = 3809.5): 1.10355785645641e-21 / 1.10355785645641e-21 = 1.0
on time step 7709 (time=3854.5), 0.033161 s/step
on time step 7829 (time=3914.5), 0.0333563 s/step
on time step 7950 (time=3975), 0.0331618 s/step
field decay(t = 4010.0): 8.42143200428359e-22 / 1.10355785645641e-21 =
0.763116492263062
on time step 8071 (time=4035.5), 0.0331005 s/step
on time step 8192 (time=4096), 0.0330795 s/step
on time step 8313 (time=4156.5), 0.0330648 s/step
field decay(t = 4210.5): 7.28424121109024e-22 / 1.10355785645641e-21 =
0.660068810028717
on time step 8435 (time=4217.5), 0.0330389 s/step
on time step 8556 (time=4278), 0.0331027 s/step
on time step 8677 (time=4338.5), 0.0331827 s/step
on time step 8798 (time=4399), 0.0331503 s/step
....
....
....
field decay(t = 701950.5): 34472.8457895311 / 63036.3907622053 =
0.54687213802539
on time step 1403914 (time=701957), 0.0340987 s/step
on time step 1404032 (time=702016), 0.0340862 s/step
on time step 1404150 (time=702075), 0.0341135 s/step
on time step 1404268 (time=702134), 0.0340847 s/step
field decay(t = 702151.0): 31455.7227755882 / 63036.3907622053 =
0.499008943805966
on time step 1404386 (time=702193), 0.034169 s/step
on time step 1404504 (time=702252), 0.034137 s/step
on time step 1404622 (time=702311), 0.0341068 s/step
field decay(t = 702351.5): 37313.2340585566 / 63036.3907622053 =
0.591931638334352
on time step 1404740 (time=702370), 0.0341675 s/step
on time step 1404858 (time=702429), 0.0341454 s/step
----------
Jingjing Li, Ph.D
Postdoc, Quantumn Science Research
Hewlett-Packard Laboratories
1501 Page Mill Road, MS 1l23
Palo Alto, CA 94304
Phone: 650 857 5454
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
