Dear Tjeerd, Steven and other meep users,
I am calculating a transmission spectrum through 2D photonic crystal and
I expirienced with the same problem.
For some parameters and when the structure is present I cannot get the
fields to decay, while normalization run is OK.
If length L is set to 16 (see .ctl file below) and n to 90 ( 90th mode --
meaning that transverse momenta equals to 2*\pi*n/sy ) the fields decay
till some value and then blow up. I tried to experiment a bit, and I have
figured out that the process can be made stable by disabling subpixel
averaging (I hope the following table could be usefull):
-----------------------------------------------
parametrs ! minimal rate !
-----------------------------------------------
res=20,sx=30, ! 1.14*10^-7 !
ph_diff=2\pi ! then blow up !
-----------------------------------------------
res=20,sx=30, ! 1.10*10^-7 !
ph_diff=0 ! then blow up !
-----------------------------------------------
res=20,sx=50, ! 1.39*10^-7 !
ph_diff=2\pi ! then blow up !
-----------------------------------------------
res=20,sx=30, ! 1.15*10^-7 !
ph_diff=2\pi ! then blow up !
pml_st=2.0 ! !
-----------------------------------------------
res=22,sx=50, ! 5.28*10^-4 !
ph_diff=2\pi ! then blow up !
-----------------------------------------------
res=22,sx=50, ! 5.26*10^-4 !
ph_diff=0 ! then blow up !
-----------------------------------------------
res=24,sx=50, ! 2.29*10^-6 !
ph_diff=2\pi ! then blow up !
-----------------------------------------------
res=20,sx=30, ! < 10^-7 !
ph_diff=2\pi ! !
eps-averaging? ! !
false ! !
-----------------------------------------------
where "res" -- resolution, "ph_diff" -- phase difference between top and
bottom of the computational cell,
"pml_st" -- strenght of pml layers.
In this example the decay rate was set to 10^-7; however for some other
parameters and structures fields decay at most till 0.001. It could be that
for the same mode they decay for L=5 or 10 and blow up for L=9.
So, try to (set! eps-averaging? false) and see whether it blows up or not.
P.S. I also noticed that setting phase difference to get periodic boundary
conditions to 0 instead of 2\pi decreases remarkably memory usage.
ctl file:
--------------------------------------------------------------------------------------------------------------------------------------------------------------
(define-param pmlthick 5)
(define-param sx (+ 30 (* 2 pmlthick)) )
(define-param sy 300)
(set! geometry-lattice (make lattice (size sx sy no-size)))
;PARAMETERS
(define-param norm false)
(define-param f 0.431) ; filling fraction f~(4 Pi r^2)/(2 Sqrt[3] a^2)
(define-param eps 14) ; dielectric
(define-param L 16)
(define-param a 1) ; lattice constant (distance between neibouhring rods)
(define-param N (+ 1 (* 2 (round (/ (* 0.5 L) (* (sqrt 3) a)) ) )))
(define-param r (* (sqrt (* f (/ (* 2 (sqrt 3)) (* 4 pi)))) a)) ; radius of
rodes
(if (not norm)
(set! geometry
(append
(geometric-objects-duplicates (vector3 (* 1.73205 a) 0) (* -1 (/ (- N 1)
2)) (/ (- N 1) 2)
(geometric-object-duplicates (vector3 0 a) 0 (+ (/ sy a) 1)
(make cylinder (center 0 (* -0.5 sy)) (radius r) (height infinity)
(material (make dielectric (epsilon eps))))
)
)
(geometric-objects-duplicates (vector3 (* 1.73205 a) 0) (* -1 (/ (- N 1)
2)) (- (/ (- N 1) 2) 1)
(geometric-object-duplicates (vector3 0 a) 0 (/ sy a)
(make cylinder (center (* 0.5 (* 1.73205 a)) (+ (* -1 (* 0.5 sy)) (*
0.5 a))) (radius r) (height infinity)
(material (make dielectric (epsilon eps))))
)
)
)
)
)
(set! pml-layers (list (make pml (thickness pmlthick) (direction X) )))
(set-param! k-point (vector3 0 0) )
(set! resolution 20)
(define-param fcen 0.4) ; pulse center frequency
(define-param df 0.06) ; pulse width (in frequency)
(define-param n 0)
(define-param K (vector3 0(/ (* (* 2 pi) n) sy)) )
(define ((pw-amp K y0) y)
(exp (* 0+1i (vector3-dot K (vector3+ y y0)))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Hz)
(center (+ pmlthick (* -0.5 sx)) 0)
(size 0 sy)
(amp-func (pw-amp K (vector3 0 (* 0.5 sy))))
)))
(define-param nfreq 100) ; number of frequencies at which to compute flux
(define trans ; transmitted flux
(add-flux fcen df nfreq
(make flux-region
(center (- (* 0.5 sx) (+ 1 pmlthick)) 0) (size 0 sy))))
(define refl ; reflected flux
(add-flux fcen df nfreq
(make flux-region
(center (+ (+ 1 pmlthick) (* -0.5 sx)) 0) (size 0
sy))))
(if (not norm) (load-minus-flux "refl-flux" refl))
(run-sources+
(stop-when-fields-decayed 50 Hz
(vector3 (- (* 0.5 sx) (+ 1 pmlthick)) (-
(* 0.5sy) 2))
1e-7) )
(if norm (save-flux "refl-flux" refl))
(display-fluxes trans refl)
-----------------------------------------------------------------------------------------------------------------------------------------------------------------
Sincerely,
Ruslan
_______________________________________________
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
