Hello all, I'm simulating PhC resonators by launching pulses into the system. Like in the bend waveguide example I put a field measurement at the output to make sure that the fields at the output are decayed enough. I put the decay value to 1e-3.
The fields generally tend to decay first, but after a while they start to ramp up to infinity. This cannot be a valid outcome, the pml's should absorb all energy in the end (or am I wrong here?), so there must be a maximum of the field values. Has anyone seen this kind of behavior before? What can be the solution? I have already started trying several things including removing all holes from the PML, using a normal gaussian pulse point source without the envelope, increasing the resolution (which seems to worsen the problem). I'm currently experimenting with the PML strength. My question, is it just the PML that is set wrong? Or do I do unhandy things with the launched pulse like ramp it on too fast? Or can this weird behavior of the simulator be explained otherwise? I attach one typical simulation (sorry for the dirty programming in the PhC functions) and some (un)nice field plots. the simulation can be run with: ./runsimulation.sh resonator005 "fcen=0.285 df=0.08 nfreq=1000 resolution=10" to obtain a reflection transmissionplot (octave ../retrplot.oct in output dir) and with: ./runsimulation.sh resonator005 to obtain the (un)nice plots below, these were obtained by running: h5topng -t 0:143 -m -1 -M 1 -Zc dkbluered -a yarg -A measure-eps-000000.00.h5 measure-hz.h5 in the output dir (I scaled them) (note: look first how much datasets are there after killing the simulation because it wont stop automatically). I hope these things can also help you to improve meep, its a nice simulator. Yours, Tjeerd Pinkert P.S. better versions of the field plots (and a few more) can be found at: http://www.scintilla.utwente.nl/~tjeerdp/study/fieldplots/
;
; PhC funtions for MEEP.
;
;;mylattice: (vector3 basis1) (vector3 basis2) (vector3 basis3) (vector3
basis-size) (vector3 size) ->
;; (#(x y z) #(x y z) #(x y z))
;;take latice parameters and return a list of lattice coordinates in cartesian
(vector3's)
(define (mylattice basis1 basis2 basis3 basis-size size)
(let ; set locale variabelen
(
(thelattice (crystal basis1 basis2 basis3 basis-size (vector3- size
(vector3 1 1 1))))
)
; transpose everything half-a-lattice down
(shiftlattice thelattice (vector3* (latticesize basis1 basis2 basis3
basis-size size) (/ -2)))
)
)
;;mysqlattice: (vector3 basis1) (vector3 basis2) (vector3 basis3) (vector3
basis-size) (vector3 size)
;; mylattice
;;takes lattice parameters but makes the lattice square with regard to vector
basis1.
(define (mysqlattice basis1 basis2 basis3 basis-size size)
(let* ; set locale variabelen
(
; maak het crystal
(thelattice (crystal basis1 basis2 basis3 basis-size (vector3- size
(vector3 1 1 1))))
; maak het vierkant
(thelattice (makelatticesquare thelattice basis1 basis2 basis3 basis-size
size))
; shift half rooster naar beneden
(thelattice (shiftlattice thelattice (vector3* (latticesize basis1 basis2
basis3 basis-size size) (/ -2))))
)
thelattice
)
)
;;shiftlattice: (mylattice lattice) (vector3 shiftvector)
;; returns: mylattice
;;shifts lattice by shiftvector
(define (shiftlattice lattice shiftvector)
(if (not (null? lattice))
(begin
(append
(list
(vector3+ (car lattice) shiftvector)
)
(shiftlattice (cdr lattice) shiftvector)
)
)
(list)
)
)
;;latticesize: (vector3 basis1) (vector3 basis2) (vector3 basis3) (vector3
basis-size) (vector3 size)
;; returns: (vector3 x-size y-size z-size)
;;calculate the size of a general lattice with input parameters
(define (latticesize basis1 basis2 basis3 basis-size size)
(vector3
(max
(* (vector3-x basis1) (vector3-x basis-size) (- (vector3-x size) 1))
(* (vector3-x basis2) (vector3-x basis-size) (- (vector3-x size) 1))
(* (vector3-x basis3) (vector3-x basis-size) (- (vector3-x size) 1))
)
(max
(* (vector3-y basis1) (vector3-y basis-size) (- (vector3-y size) 1))
(* (vector3-y basis2) (vector3-y basis-size) (- (vector3-y size) 1))
(* (vector3-y basis3) (vector3-y basis-size) (- (vector3-y size) 1))
)
(max
(* (vector3-z basis1) (vector3-z basis-size) (- (vector3-z size) 1))
(* (vector3-z basis2) (vector3-z basis-size) (- (vector3-z size) 1))
(* (vector3-z basis3) (vector3-z basis-size) (- (vector3-z size) 1))
)
)
)
;;mirror-x: (mylattice lattice)
;; returns: (mylattice lattice)
;;mirrors a list of points in the x-axis
(define (mirror-x lattice)
(if (not (null? lattice))
(begin
(append
(list (vector3 (- (vector3-x (car lattice))) (vector3-y (car lattice))
(vector3-z (car lattice))))
(mirror-x (cdr lattice))
)
)
(list)
)
)
;;mirror-y: (mylattice lattice)
;; returns: (mylattice lattice)
;;mirrors a list of points in the y-axis
(define (mirror-y lattice)
(if (not (null? lattice))
(begin
(append
(list (vector3 (vector3-x (car lattice)) (- (vector3-y (car lattice)))
(vector3-z (car lattice))))
(mirror-y (cdr lattice))
)
)
(list)
)
)
;;mirror-z: (mylattice lattice)
;; returns: (mylattice lattice)
;;mirrors a list of points in the z-axis
(define (mirror-z lattice)
(if (not (null? lattice))
(begin
(append
(list (vector3 (vector3-x (car lattice)) (vector3-y (car lattice)) (-
(vector3-z (car lattice)))))
(mirror-z (cdr lattice))
)
)
(list)
)
)
;;putonlattice: (mylattice lattice) (geometry unitcell)
;; returns: geometry
;;shifts a centered geometry "unitcell" on each latticepoint and returns the
resulting structure as a geometry list
(define (putonlattice lattice unitcell)
(if (not (null? lattice))
(begin
(append
(shift-geometric-objects unitcell (car lattice))
(putonlattice (cdr lattice) unitcell)
)
)
(list)
)
)
;;gaussian-gaussian-source: frequency fwidth component center (size gwidth)
;; returns: the (make source ....) output.
;;makes a spatially gaussian and temporally gaussian shaped source
;;note that the variable gwidth^ needs to exist globally
(define (gaussian-gaussian-source freq df comp cntr gwidth)
(define tl (vector3* gwidth 7)); the size of the source must be bigger than
; the width of the gaussian function.
(set! gwidth^ (vector3-norm gwidth)) ; the extra parameter for
gaussian-beam...
(make source
(src
(make gaussian-src
(frequency freq)
(fwidth df)
)
)
(component comp)
(center cntr)
(size tl)
(amp-func gaussian-beam)
)
)
;;gaussian-continuous-source: frequency component center gwidth
;; returns: the (make source ....) output.
;;makes a spatially gaussian shaped and temporally continues wave source
;;note that the variable gwidth^ needs to exist globally
(define (gaussian-continuous-source freq comp cntr gwidth)
(define tl (vector3* gwidth 7)); the size of the source must be bigger than
; the width of the gaussian function.
(set! gwidth^ (vector3-norm gwidth)) ; the extra parameter for
gaussian-beam...
(make source
(src
(make continuous-src
(frequency freq)
(component comp)
)
)
(center cntr)
(size tl)
(amp-func gaussian-beam)
)
)
;;
;; Functions below here are private to the above functions.
;;
;;
;; shift a list of geometric-objects, a unitcell is a list of geometric objects
(define (shift-geometric-objects objectlist shiftvector)
(if (not (null? objectlist))
(begin
(append
(list (shift-geometric-object (car objectlist) shiftvector))
(shift-geometric-objects (cdr objectlist) shiftvector)
)
)
(list)
)
)
;; make a no square lattice square with regard to basis1 vector
(define (makelatticesquare lattice basis1 basis2 basis3 basis-size size)
(let
(
(basis1-length (* (sqrt (+ (* (vector3-x basis1) (vector3-x basis1)) (*
(vector3-y basis1) (vector3-y basis1)) (* (vector3-z basis1) (vector3-z
basis1)))) (vector3-x basis-size) (vector3-x size)))
(shiftvector (vector3- (vector3 0 0 0) (vector3* basis1 (* (vector3-x
basis-size) (vector3-x size)))))
(basis1-phi (angle-phi basis1))
(basis1-theta (angle-theta basis1))
)
(begin
; loop through lattice list
(if (not (null? lattice))
(begin
(let*
(
; make a projection of the latticepoint on the basisvector.
(latpoint (car lattice))
; calculate the length of the latticepointvector (from 0 0 0)
(latpointlength (sqrt (+ (* (vector3-x latpoint) (vector3-x
latpoint)) (* (vector3-y latpoint) (vector3-y latpoint)) (* (vector3-z
latpoint) (vector3-z latpoint)))))
; calculate the difference of the phi angles
(diffangle-phi (abs (- basis1-phi (angle-phi latpoint))))
; calculate the difference of the theta angles
(diffangle-theta (abs (- basis1-theta (angle-theta latpoint))))
; calculate projectionlength on the z-basis1 plane
(z1length=latpointlength*cos(phi1-phi2))
(projlength1 (* latpointlength (cos diffangle-phi)))
; calculate projectionlength on the basis1 vector
(projectionlength=z1length*cos(theta1-theta2))
(projlength2 (* projlength1 (cos diffangle-theta)))
; if the length of this projection is longer than the basisvector
)
(begin
(if (or (not (< projlength2 basis1-length)) (are-close?
projlength2 basis1-length (* 1e-3 (vector3-x basis-size))))
; transpose crystalpoint one crystallength along basis1 vector.
(set! latpoint (car (shiftlattice (list latpoint) shiftvector)))
)
(append
(list latpoint)
(makelatticesquare (cdr lattice) basis1 basis2 basis3
basis-size size)
)
)
)
)
(list)
)
)
)
)
;two numbers are closer than tolerance
(define (are-close? num1 num2 tolerance)
(if (and (< num2 (+ num1 tolerance)) (> num2 (- num1 tolerance)))
#t
#f
)
)
;the angle of the x-axis in the x-y plane
(define (angle-phi vector)
(if (= 0 (vector3-x vector)) ;to avoid division by 0
0
(atan (/ (vector3-y vector) (vector3-x vector)))
)
)
;the angle with the projection of 'vector' in the x-y plane
(define (angle-theta vector)
(if (vector3= vector (vector3 0 0 0))
0
; note: only because the (sqrt 0) is 0.0 this works, else we would get
overflows
; see angle-phi how to solve this if it happens
(atan (/ (vector3-z vector) (sqrt (+ (* (vector3-y vector) (vector3-y
vector)) (* (vector3-x vector) (vector3-x vector))))))
)
)
;; calculate two vectors that define the perpendicular plane of the vector1
input
(define (perpendicular-plane vector1 vector2)
(let
(
(vlakvector1 (vector3-cross vector1 vector2))
)
(list vlakvector1 (vector3-cross vector1 vlakvector1))
)
)
;; define most extended position in this crystal...
(define (crystalsitecoords basis1 basis2 basis3 basis-size position)
(vector3
(+
(* (vector3-x basis1) (vector3-x basis-size) (vector3-x position))
(* (vector3-x basis2) (vector3-y basis-size) (vector3-y position))
(* (vector3-x basis3) (vector3-z basis-size) (vector3-z position))
)
(+
(* (vector3-y basis1) (vector3-x basis-size) (vector3-x position))
(* (vector3-y basis2) (vector3-y basis-size) (vector3-y position))
(* (vector3-y basis3) (vector3-z basis-size) (vector3-z position))
)
(+
(* (vector3-z basis1) (vector3-x basis-size) (vector3-x position))
(* (vector3-z basis2) (vector3-y basis-size) (vector3-y position))
(* (vector3-z basis3) (vector3-z basis-size) (vector3-z position))
)
)
)
;crystalcolumn(n)=crystalcolumn(n-1)+latticesite(n)
(define (crystalcolumn basis1 basis2 basis3 basis-size size)
(if (< (vector3-x size) 1)
(list
(crystalsitecoords basis1 basis2 basis3 basis-size size)
)
(append
(list
(crystalsitecoords basis1 basis2 basis3 basis-size size)
)
(let
((size (vector3 (- (vector3-x size) 1) (vector3-y size) (vector3-z
size))))
(crystalcolumn basis1 basis2 basis3 basis-size size)
)
)
)
)
;crystalslab(n)=crystalslab(n-1)+crystalcolumn(n)
(define (crystalslab basis1 basis2 basis3 basis-size size)
(if (< (vector3-y size) 1)
(crystalcolumn basis1 basis2 basis3 basis-size size)
(append
(crystalcolumn basis1 basis2 basis3 basis-size size)
(let
((size (vector3 (vector3-x size) (- (vector3-y size) 1) (vector3-z
size))))
(crystalslab basis1 basis2 basis3 basis-size size)
)
)
)
)
;crystal(n)=crystal(n-1)+crystalslab(n)
(define (crystal basis1 basis2 basis3 basis-size size)
(if (< (vector3-z size) 1)
(crystalslab basis1 basis2 basis3 basis-size size)
(append
(crystalslab basis1 basis2 basis3 basis-size size)
(let
((size (vector3 (vector3-x size) (vector3-y size) (- (vector3-z size)
1))))
(crystal basis1 basis2 basis3 basis-size size)
)
)
)
)
; unitcell op lattice
; (geometry-lattice geometry lattice)
; lees een coordinaat uit lattice
; maak een unitcell
; (shift-geometric-object unitcell coordinaat)
; stop het verschoven object in de geometry
; define aplitude functions for the sources for the simulator
; NOTE: to use this function the variable gwidth^ must be defined in the
environment.
(define (gaussian-beam point)
(define x (vector3-x point))
(define y (vector3-y point))
(define z (vector3-z point))
(exp (/ (+ (* x x -1) (* y y -1) (* z z -1)) (* gwidth^ gwidth^)))
)
(load "phc-funcs.ctl")
; this is ugly but fixes our stack overflow problem, we should rewrite the
phc-funcs.ctl file
; to use proper tail recursive functions or to avoid excessive stack usage at
all by head
; recursion...
(debug-set! stack 2000000)
;;
;; commandline params
;;
; define if we run our "gauge" run or not
(define-param r 135)
(define-param a 440)
(define-param gaugerun #f)
(define-param testrun #f)
(define-param si-n 2.9)
(define-param fcen 0.261)
(define-param df 0.02)
(define-param nfreq 500)
(define-param fielddecay 1e-3) ; set higher to avoid unstable situation
(define-param pmlstrength 1.0) ; set higher to avoid unstable situation
;resolution 20 (this is set after the structure is defined)
;normalize r/a ratio
(set! r (/ r a))
; define the lattice parameters for our triangular lattice
(define mybasis1 (vector3 1 0 0))
(define mybasis2 (vector3 (/ 2) (/ (sqrt 3) 2) 0))
(define mybasis3 (vector3 0 0 1))
(define mybasis-size (vector3 1 1 1))
(define mysize (vector3 41 41 1))
(define mylat (mysqlattice mybasis1 mybasis2 mybasis3 mybasis-size mysize))
(define guide1
(shiftlattice (mylattice mybasis1 mybasis2 mybasis3 mybasis-size (vector3 18
1 1)) (vector3 -11.5 0 0))
)
(define guide2
(shiftlattice (mylattice mybasis1 mybasis2 mybasis3 mybasis-size (vector3 18
1 1)) (vector3 12.5 0 0))
)
(define defect
(list (vector3 0 0 0) (vector3 1 0 0) (vector3 0.5 (/ (sqrt 3) -2)) (vector3
0.5 (/ (sqrt 3) 2)))
)
(define guidegaugerun
(mylattice mybasis1 mybasis2 mybasis3 mybasis-size (vector3 41 1 1))
)
; define the structure
(set! geometry-lattice (make lattice (size 40 30 no-size)))
(set! default-material (make dielectric (epsilon (* si-n si-n))))
(define wgmaterial (make dielectric (epsilon (* si-n si-n))))
(define testmaterial (make dielectric (epsilon 4)))
(define unitcell
(list
(make cylinder (center 0 0 0) (radius r) (height infinity) (material air))
)
)
(define fillholes
(list
(make cylinder (center 0 0 0) (radius r) (height infinity) (material
wgmaterial))
)
)
(define mylatticegeom (putonlattice mylat unitcell))
(define guide1geom (putonlattice guide1 fillholes))
(define guide2geom (putonlattice guide2 fillholes))
(define defectgeom (putonlattice defect fillholes))
(define gaugerungeom (putonlattice guidegaugerun fillholes))
; now the simulator parameters...
(if gaugerun
(set! geometry
(append
mylatticegeom
gaugerungeom
)
)
(set! geometry
(append
mylatticegeom
guide1geom
guide2geom
defectgeom
; this was to check the position of the fluxplane and field measurement
spot.
;(list
; (make cylinder (center 8 0 0) (radius r) (height infinity) (material
testmaterial))
;)
)
)
)
(define gwidth^ 0) ; we need this variable for our own gaussian source
; the sources
(set! sources
(list
(gaussian-gaussian-source fcen df Hz (vector3 -17 0 0) (vector3 0 .707 0))
)
)
; the rest...
(set! pml-layers (list (make pml (thickness 1.0) (strength pmlstrength))))
(set-param! resolution 20)
;set output directory name
(define output-dir (string-append (number->string si-n) "-" (number->string
fcen) "-" (number->string df) "-" (number->string nfreq) "-" (number->string
resolution)))
(use-output-directory output-dir)
;prepend gauge/measure to filenames to distinguish both runs...
(if gaugerun
(set! filename-prefix "gauge")
(set! filename-prefix "measure")
)
; write the output directory to a file for our bash script...
(define dirfile (open-output-file "simdir.txt"))
(display output-dir dirfile)
(close-output-port dirfile)
(define in-fwd-center (vector3 -6 0 0))
(define in-bwd-center (vector3 -6 0 0))
(define out-fwd-center (vector3 8 0 0))
; define flux regions after the rest of the simulation is completely set up
(define in-fwd ; the forward input flux
(add-flux fcen df nfreq
(make flux-region (center in-fwd-center) (size 0 15))
)
)
(define in-bwd ; the reflected input flux
(add-flux fcen df nfreq
(make flux-region (center in-bwd-center) (size 0 15))
)
)
; the output port fluxes
(define out-fwd
(add-flux fcen df nfreq
(make flux-region (center out-fwd-center) (size 0 15))
)
)
; Als we de reflectie willen meten, moeten we de flux van de gaugerun
; (die wordt opgeslagen) eraf trekken
(if (not testrun)
(if (not gaugerun)
(begin
(set! filename-prefix #f)
(load-minus-flux "gauge-in-bwd" in-bwd)
(set! filename-prefix "measure")
)
)
)
(run-sources+
(stop-when-fields-decayed 50 Hz
(if gaugerun
(vector3 17 0 0)
out-fwd-center
)
fielddecay)
(at-beginning output-epsilon)
(to-appended "hz" (at-every (* fcen 500.1) output-hfield-z))
)
; sla de flux van de gaugerun op om een reflectiemeting te doen.
(if gaugerun (save-flux "in-bwd" in-bwd))
(display-fluxes in-fwd in-bwd out-fwd)
runsimulation.sh
Description: application/shellscript
gaugerun = dlmread "fluxdata.gaugerun"; measurerun = dlmread "fluxdata.measurerun"; axis([0,1,0,2],"autox"); plot(gaugerun(:,2),(measurerun(:,5)./gaugerun(:,3)), ";Transmission;"); hold on; plot(gaugerun(:,2),-1.*(measurerun(:,4)./gaugerun(:,3)), ";Reflection;"); pause;
<<inline: measure-hz.t006.jpg>>
<<inline: measure-hz.t048.jpg>>
_______________________________________________ meep-discuss mailing list meep-discuss@ab-initio.mit.edu http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
