Dear Meep's users,
I'm a PHD student at the Laboratory of Photonics and Nanostructure (CNRS, Franc
e), and I'm trying to simulate the photonic cavity L3. My simulations sometimes
blow up. When I turn on both eps-averaging and symmetries, the fields of the c
avity diverge exponentially. When I turn off the eps-averaging or the symmetry
along the X direction it goes right.
eps-averaging and symmetry seems to be partially incompatible?
I could run the calculations without the symmetry, but it greatly lengthens the
time to calculate them, and for what I am planning to do, time is critical (I
should calculate many of them). With eps-averaging turned off the center of the
holes are systematically centered on the calculations grid (I don't know why).
Therefore shift=0.06 and shift=0.07 give two strictly identical eps-fields whe
n eps-averaging is turned off (shift is a personnal parameter to shift some hol
es). Such problem doesn't appear with eps-averaging turned on which is why I wa
nt to use this option.
I enclosed a "light" version of my ctl file. For example, with:
meep run-number=1 subpixel=1 symX=true symY=true symZ=true instabilite.ctl > in
stabilite-1.log
meep run-number=2 subpixel=1 symX=false symY=true symZ=true instabilite.ctl > i
nstabilite-2.log
meep run-number=3 subpixel=0 symX=true symY=true symZ=true instabilite.ctl > in
stabilite-3.log
run number 2 and 3 go fine, but number 1 numerically diverges (subpixel=0 turns
off the eps-averaging)
I searched in forums if someone encountered something like that. Some instabili
ty can be due to the Courant, but in my case I think its value is good, because
if it wasn't, my simulations should allways diverge and not only when eps-aver
aging is on, isn't it?
Do somebody have any clue for my problem?
Thank you very much for whatever time you should spend for me,
Matthieu Larqué
LPN CNRS France
File instability.ctl:
; <<<<<<<<<<<<<< Some parameters to describe the geometry <<<<<<<<<<<<<<
; dielectric constant of waveguide, index squared
(define-param eps (* 3.5 3.5))
; dielectric constant for holes
(define-param special (* 1 1))
(define-param Ra 0.42) ; period, the prefix 'R' means 'real value'
(define-param Rr 0.12) ; radius of holes
(define-param Rhigh 0.25) ; membrane thickness
(define-param N 4) ; 2*N rows of holes around the center
(define-param bord 2) ; the equivalent in distance of 'bord*Ra' of
; holes around the cristal of membrane without holes
(define-param shift 0.18) ; shift of holes
(define-param dr 0.9) ; modified holes of radii Rr*dr
(define-param Rlambda 1.6) ; pulse center wavelength
(define-param noscil 10) ; roughly the number of temporal oscillations
; of the source inside its halwidth intensity
(define-param src-cutoff 5.0) ; when turning on and off gaussian source
(define-param temps 300) ; time after source turned off
(define-param resol 15) ; number of pixels per unit
(define-param distpml 3.0) ; distance between the membrane and pml, in
; normalized units
(define-param subpixel 1) ; sort of more global parameter to control
; eps-averaging and its accuracy
(define-param symX true) ; do we use x symetry?
(define-param symY true) ; idem y
(define-param symZ true) ; idem z
(define-param run-number 0) ; run number
(set! filename-prefix (string-append (get-filename-prefix) "-" (number->string
run-number)))
; <<<<<<<<<<<<<< Normalised values and meep's parameters <<<<<<<<<<<<<<
(set-param! resolution 15)
(set-param! Courant (* 0.03 resolution))
(define a 1) ; lattice
(define r (/ Rr Ra)) ; radius of holes
(define high (/ Rhigh Ra)) ; membrane thickness
(define r-petit (* r dr)) ; radius of little holes
; size of cell in x direction (pair integer number of pixels, PML excluded)
(define sx (* 2 (+ (* 2 N) bord)))
(set! sx (/ (+ (floor (* sx resolution 0.5)) 1) (* resolution 0.5)))
; size of cell in y direction (pair number of pixels, PML excluded)
(define sy (* 2 (+ (* (sqrt 3) N) bord)))
(set! sy (/ (+ (floor (* sy resolution 0.5)) 1) (* resolution 0.5)))
; size of cell in z direction (pair number of pixels, PML excluded)
(define sz (+ distpml high distpml))
(set! sz (/ (+ (floor (* sz resolution 0.5)) 1) (* resolution 0.5)))
(define dpml 1) ; PML thickness
(define pmlstrength 1.0) ; PML strength
(define fcen (/ Ra Rlambda)) ; central frequency of the source
(define df (/ fcen noscil)) ; frequency width of the source
(set! geometry-lattice (make lattice (size (+ dpml sx dpml) (+ dpml sy dpml)
(+ dpml sz dpml))))
(set! progress-interval 100)
(set! ensure-periodicity false)
(set! eps-averaging? (> subpixel 0))
(set-param! subpixel-tol (/ 0.0001 subpixel))
(set-param! subpixel-maxeval (inexact->exact (* 100000 subpixel)))
; <<<<<<<<<<<<<< Geometrical elements <<<<<<<<<<<<<<
; Convention: zero coordinates at center of the membrane
; basic material
(define mat (make dielectric (epsilon eps)))
; air
(define matspec (make dielectric (epsilon special)))
; membrane
(define membrane (make block (center 0 0 0) (size infinity infinity high)
(material mat)))
; the whole world
(define tout (make block (center 0 0 0) (size infinity infinity infinity)
(material matspec)))
; a hole in the membrane
(define trou (make cylinder (center 0 0 0) (radius r) (height high)
(material matspec)))
; a little hole in the membrane
(define trou-petit (make cylinder (center 0 0 0) (radius r-petit) (height high)
(material matspec)))
; refilling a hole
(define rempli (make cylinder (center 0 0 0) (radius r) (height high)
(material mat)))
(define repy (* a (sqrt 3)))
(define ligne-2N+1 (lambda (ypos) (geometric-object-duplicates (vector3 a 0 0)
(- 0 (* 2 N)) (* 2 N) (shift-geometric-object trou (vector3 0 ypos
0)))))
(define ligne-2N (lambda (ypos) (geometric-object-duplicates (vector3 a 0 0)
(- 0 (* 2 N)) (- (* 2 N) 1) (shift-geometric-object trou
(vector3 (/ a 2) ypos 0)))))
(define cristal (append (geometric-objects-duplicates (vector3 0 repy 0)
(- 0 N) N (ligne-2N+1 0)) (geometric-objects-duplicates (vector3 0 repy
0)
(- 0 N) (- N 1) (ligne-2N (/ repy 2)))))
(define shiftx-L3 (+ (* 2 a) shift))
(define L3m (list rempli
(shift-geometric-object rempli (vector3 a 0 0))
(shift-geometric-object rempli (vector3 (- 0 a) 0 0))
(shift-geometric-object rempli (vector3 (* 2 a) 0 0))
(shift-geometric-object rempli (vector3 (- 0 (* 2 a)) 0 0))
(shift-geometric-object trou-petit (vector3 shiftx-L3 0 0))
(shift-geometric-object trou-petit (vector3 (- 0 shiftx-L3) 0 0))
))
; <<<<<<<<<<<<<< Geometry <<<<<<<<<<<<<<
(set! geometry (append (list tout membrane) cristal L3m))
(set! pml-layers (list
(make pml (direction Y) (thickness dpml) (strength pmlstrength))
(make pml (direction X) (thickness dpml) (strength pmlstrength))
(make pml (direction Z) (thickness dpml) (strength pmlstrength))
))
(set! symmetries (list))
(if symX (set! symmetries (append symmetries (list (make mirror-sym
(direction X) (phase +1))))))
(if symY (set! symmetries (append symmetries (list (make mirror-sym
(direction Y) (phase -1))))))
(if symZ (set! symmetries (append symmetries (list (make mirror-sym
(direction Z) (phase +1))))))
; <<<<<<<<<<<<<< Sources <<<<<<<<<<<<<<
(define source-gauss (make gaussian-src (frequency fcen) (fwidth df)
(cutoff src-cutoff)))
(set! sources (list (make source (src source-gauss) (component Ey)
(center 0 0 0))))
; <<<<<<<<<<<<<< Timings <<<<<<<<<<<<<<
(define ts (/ Courant resolution)) ; time step
(define flag-src 1)
; sera tourne à 0 quand la source s'arrete
; <<<<<<<<<<<<<< Time field <<<<<<<<<<<<<<
; this section to output in a file the meep-time, the value of the evolution
; of the field in the center of the cavity, and if the source is on
(define cout-line (lambda (line port) (begin
(display line port)
(newline port)
(force-output port)
)))
(define cout-current (lambda (line) (cout-line line (current-output-port))))
(define cout-timefield (open-output-file (string-append (get-filename-prefix)
"-time.csv")))
(define display-timefield (lambda ()
(cout-line (string-append (number->string (meep-time))
", " (number->string (get-field-point Ey (vector3 0 0 0)))
", " (number->string flag-src)) cout-timefield)
))
; <<<<<<<<<<<<<< Run <<<<<<<<<<<<<<
(define src-end 0) ; meep-time when the source turned off
(define run-end 0) ; meep-time when the simulation should stop
(define (when-source-off) (begin
(set! flag-src 0)
(set! src-end (meep-time))
(set! run-end (+ src-end temps))
))
; maximum value of the field at the center of the cavity when source is on
(define src-max 0)
(define (max-val) (set! src-max (max src-max (get-field-point Ey
(vector3 0 0 0)))))
(define (test-end-evol) (or (< run-end (meep-time))
(< src-max (get-field-point Ey (vector3 0 0 0)))))
(define (when-end-evol) (begin
(if (< src-max (get-field-point Ey (vector3 0 0 0)))
(cout-current "Instability"))
(output-efield)
(output-hfield)
))
; <<<<<<<<<<<<<< Run <<<<<<<<<<<<<<
(cout-current "Run: source pulse" )
(run-sources
(at-beginning output-epsilon)
(at-every ts max-val)
(at-every ts display-timefield)
(at-end when-source-off)
)
(cout-current "Run: mode evolution" )
; will stop when the expected time was simulated,
; or when the field goes stronger than the maximum value reached
; when the source was on
; on the second case, "Instability" is printed on the log file
(run-until test-end-evol
(at-every ts display-timefield)
(at-end when-end-evol)
)
(close-output-port cout-timefield)
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
