mpb.scm.in

麻省理工的计算光子晶体的程序

; variables holding the band range data and current list of gaps, in
; the format returned by update-band-range-data and output-gaps, above:
(define band-range-data '())
(define gap-list '())

; Return the frequency gap from the band #lower-band to the band
; #(lower-band+1), as a percentage of mid-gap frequency.  The "gap"
; may be negative if the maximum of the lower band is higher than the
; minimum of the upper band.  (The gap is computed from the
; band-range-data of the previous run.)
(define (retrieve-gap lower-band)
  (if (> (+ lower-band 1) (length band-range-data))
      (error "retrieve-gap called for higher band than was calculated")
      (let ((f1 (cadr (list-ref band-range-data (- lower-band 1))))
	    (f2 (caar (list-ref band-range-data lower-band))))
	(/ (- f2 f1) (* 0.005 (+ f1 f2))))))

; ****************************************************************

; stuff to keep statistics on the eigensolver performance, for tuning:
(define eigensolver-iters '()) ; the iterations used, updated by (run)
(define total-run-time 0.0) ; the total time used by (run) functions (seconds)

(define (display-eigensolver-stats)
  (let ((num-runs (length eigensolver-iters)))
    (if (> num-runs 0)
	(let ((min-iters (apply min eigensolver-iters))
	      (max-iters (apply max eigensolver-iters))
	      (mean-iters (/ (fold-right + 0 eigensolver-iters) num-runs)))
	  (print "eigensolver iterations for " num-runs " k-points: "
		 min-iters "-" max-iters ", mean = "  mean-iters)
	  (if (defined? 'sort)  ; sort was added in Guile 1.3.x
	      (let ((sorted-iters (sort eigensolver-iters <)))
		(let ((median-iters (* 0.5 (+ (list-ref sorted-iters
							(quotient num-runs 2))
					      (list-ref sorted-iters
							(- (quotient 
							    (+ num-runs 1) 2)
							   1))))))
		  (print ", median = " median-iters))))
	  (print "\nmean flops per iteration = "
		 (/ eigensolver-flops (* num-runs mean-iters)) "\n")
	  (print "mean time per iteration = "
		 (/ total-run-time (* mean-iters num-runs)) " s\n")))))

; ****************************************************************

; Define an easy way for the user to split the k-points list over multiple
; processes.  k-split-num is the number of chunks to split the k-points into,
; and k-split-index is the index of the current chunk (0 to k-split-num - 1).
(define-param k-split-num 1)
(define-param k-split-index 0)

; Split a list L into num more-or-less equal pieces, returning the piece
; given by index (in 0..num-1), along with the index in L of the first
; element of the piece, as a car pair: (first-index . piece-of-L).
(define (list-split L num index)
  (define (list-sub L start len index rest)
    (if (null? L)
	(reverse rest)
	(if (and (>= index start) (< index (+ start len)))
	    (list-sub (cdr L) start len (+ index 1) (cons (car L) rest))
	    (list-sub (cdr L) start len (+ index 1) rest))))
  (if (or (>= index num) (negative? index))
      (cons (length L) '())
      (let ((block-size (quotient (+ (length L) num -1) num)))
	(let ((start (* index block-size))
	      (len (min block-size (- (length L) (* index block-size)))))
	  (cons start (list-sub L start len 0 '()))))))

; ****************************************************************

(define current-k (vector3 0)) ; current k point in the run function
(define all-freqs '()) ; list of all freqs computed in a run

; (run) functions, to do vanilla calculations.  They all take zero or
; more "band functions."  Each function should take a single
; parameter, the band index, and is called for each band index at
; every k point.  These are typically used to output the bands.

(define (run-parity p reset-fields . band-functions)
 (set! total-run-time (+ total-run-time
  (begin-time "total elapsed time for run: "
   (set! all-freqs '())
   (set! band-range-data '())
   (set! interactive? false)  ; don't be interactive if we call (run)
   (begin-time "elapsed time for initialization: "
	       (init-params p (if reset-fields true false))
	       (if (string? reset-fields) (load-eigenvectors reset-fields)))
   (let ((k-split (list-split k-points k-split-num k-split-index)))
     (set-kpoint-index (car k-split))
     (if (zero? (car k-split))
	 (output-epsilon)) ; output epsilon immediately for 1st k block
     (if (> num-bands 0)
	 (begin
	   (map (lambda (k)
		  (set! current-k k)
		  (begin-time "elapsed time for k point: " (solve-kpoint k))
		  (set! all-freqs (cons freqs all-freqs))
		  (set! band-range-data 
			(update-band-range-data band-range-data freqs k))
		  (set! eigensolver-iters
			(append eigensolver-iters
				(list (/ iterations num-bands))))
		  (map (lambda (f)
			 (if (zero? (procedure-num-args f))
			     (f) ; f is a thunk: evaluate once per k-point
			     (do ((band 1 (+ band 1))) ((> band num-bands))
			       (f band))))
		       band-functions))
		(cdr k-split))
	   (if (> (length (cdr k-split)) 1)
	       (begin
		 (output-band-range-data band-range-data)
		 (set! gap-list (output-gaps band-range-data)))
	       (set! gap-list '()))))))))
 (set! all-freqs (reverse all-freqs)) ; put them in the right order
 (print "done.\n"))

(define run-polarization run-parity) ; backwards compatibility

; a macro to create a run function with a given name and parity
(defmacro-public define-run (name parity)
  `(define (,name . band-functions)
     (apply run-parity (append (list ,parity true) band-functions))))

(define-run run NO-PARITY)
(define-run run-zeven EVEN-Z)
(define-run run-zodd ODD-Z)
(define-run run-yeven EVEN-Y)
(define-run run-yodd ODD-Y)
(define-run run-yeven-zeven (+ EVEN-Y EVEN-Z))
(define-run run-yeven-zodd (+ EVEN-Y ODD-Z))
(define-run run-yodd-zeven (+ ODD-Y EVEN-Z))
(define-run run-yodd-zodd (+ ODD-Y ODD-Z))

(define run-even run-zeven) ; backwards compatibility
(define run-odd run-zodd) ; backwards compatibility
(define run-te run-zeven)
(define run-tm run-zodd)
(define run-te-yeven run-yeven-zeven)
(define run-te-yodd run-yodd-zeven)
(define run-tm-yeven run-yeven-zodd)
(define run-tm-yodd run-yodd-zodd)

; ****************************************************************

; Some predefined output functions (functions of the band index),
; for passing to (run).

(define (output-hfield which-band)
  (get-hfield which-band)
  (output-field))
(define (output-hfield-x which-band)
  (get-hfield which-band)
  (output-field-x))
(define (output-hfield-y which-band)
  (get-hfield which-band)
  (output-field-y))
(define (output-hfield-z which-band)
  (get-hfield which-band)
  (output-field-z))

(define (output-dfield which-band)
  (get-dfield which-band)
  (output-field))
(define (output-dfield-x which-band)
  (get-dfield which-band)
  (output-field-x))
(define (output-dfield-y which-band)
  (get-dfield which-band)
  (output-field-y))
(define (output-dfield-z which-band)
  (get-dfield which-band)
  (output-field-z))

(define (output-efield which-band)
  (get-efield which-band)
  (output-field))
(define (output-efield-x which-band)
  (get-efield which-band)
  (output-field-x))
(define (output-efield-y which-band)
  (get-efield which-band)
  (output-field-y))
(define (output-efield-z which-band)
  (get-efield which-band)
  (output-field-z))

(define (output-hpwr which-band)
  (get-hfield which-band)
  (compute-field-energy)
  (output-field))

(define (output-dpwr which-band)
  (get-dfield which-band)
  (compute-field-energy)
  (output-field))

(define (get-poynting which-band)
  (get-efield which-band)                      ; put E in cur-field
  (let ((e (field-copy cur-field)))            ; ... and copy to local var.
    (get-hfield which-band)                    ; put H in cur-field
    (field-map! cur-field                      ; write ExH to cur-field
		(lambda (e h) (vector3-cross (vector3-conj e) h))
		e cur-field)
    (cvector-field-nonbloch! cur-field)))
(define (output-poynting which-band)
  (get-poynting which-band)
  (output-field-to-file -1 (string-append filename-prefix "flux.")))
(define (output-poynting-x which-band)
  (get-poynting which-band)
  (output-field-to-file 0 (string-append filename-prefix "flux.")))
(define (output-poynting-y which-band)
  (get-poynting which-band)
  (output-field-to-file 1 (string-append filename-prefix "flux.")))
(define (output-poynting-z which-band)
  (get-poynting which-band)
  (output-field-to-file 2 (string-append filename-prefix "flux.")))

(define (get-tot-pwr which-band)
  (get-dfield which-band)
  (compute-field-energy)
  (let ((epwr (field-copy cur-field))
        (tot-pwr (rscalar-field-make cur-field)))
    (get-hfield which-band)
    (compute-field-energy)
    (field-map! tot-pwr
                (lambda (epwr hpwr) (+ epwr hpwr))
                epwr cur-field)
    (field-load tot-pwr)))
(define (output-tot-pwr which-band)
  (get-tot-pwr which-band)
  (output-field-to-file -1 (string-append filename-prefix "tot.")))

; We need a special function to evaluate band functions, since
; band functions can either be a function of the band number or
; a thunk (function of no arguments, evaluated once per k-point).
(define (apply-band-func-thunk band-func which-band eval-thunk?)
  (if (zero? (procedure-num-args band-func))
      (if eval-thunk? (band-func)) ; evaluate thunks once per k-point
      (band-func which-band)))
(define (apply-band-func band-func which-band)
  (apply-band-func-thunk band-func which-band (= which-band 1)))

; The following function returns an output function that calls
; output-func for bands with D energy in objects > min-energy.
; For example, (output-dpwr-in-objects output-dfield 0.20 some-object)
; would return an output function that would spit out the D field
; for bands with at least %20 of their D energy in some-object.
(define (output-dpwr-in-objects output-func min-energy . objects)
  (lambda (which-band)
    (get-dfield which-band)
    (compute-field-energy)
    (let ((energy (compute-energy-in-object-list objects)))
        ; output the computed energy for grepping:
	(print "dpwr:, " which-band ", "
	       (list-ref freqs (- which-band 1)) ", " energy "\n")
	(if (>= energy min-energy)
	    (apply-band-func output-func which-band)))))

; Combines zero or more band functions into one:
(define (combine-band-functions . band-funcs)
  (lambda (which-band)
    (map (lambda (f) (apply-band-func f which-band)) band-funcs)))

; Only invoke the given band functions for the specified k-point:
(define (output-at-kpoint kpoint . band-funcs)
  (let ((band-func (apply combine-band-functions band-funcs)))
    (lambda (which-band)
      (if (vector3= current-k kpoint)
	  (band-func which-band)))))

; Band functions to pick a canonical phase for the eigenstate of the
; given band based upon the spatial representation of the given field:
(define (fix-hfield-phase which-band)
  (get-hfield which-band)
  (fix-field-phase))
(define (fix-dfield-phase which-band)
  (get-dfield which-band)
  (fix-field-phase))
(define (fix-efield-phase which-band)
  (get-efield which-band)
  (fix-field-phase))

; ****************************************************************
; Here, we solve the inverse problem, that of solving for the
; wavevectors for a set of bands at a given frequency.  To do
; this, we use the fact that we can compute the group velocities
; cheaply, and thus can employ find-root-deriv (Newton's method).
; Moreover, we save information gathered while finding the k's of
; higher bands to speed the computation for lower bands.

(define (find-k p omega band-min band-max kdir
		tol kmag-guess kmag-min kmag-max . band-funcs)
  (define (ncdr n lst) (if (> n 0) (ncdr (- n 1) (cdr lst)) lst))
  (let ((num-bands-save num-bands) (k-points-save k-points)
	(nb (- band-max band-min -1)) 
	(kdir1 (cartesian->reciprocal (unit-vector3 (reciprocal->cartesian kdir))))
	; k0s is an array caching the best k value found for each band:
	(k0s (if (list? kmag-guess) (list->vector kmag-guess)
		 (make-vector (- band-max band-min -1) kmag-guess)))
	; bktab is a table (assoc. list) to memoize all (band . k) results:
	(bktab '()))
    (define ((rootfun b) k)
      (let ((tab-val (assoc (cons b k) bktab))) ; first, look in cached table
	(if tab-val 
	    (begin ; use cached result if available
	      (print "find-k " b " at " k ": " (cadr tab-val) " (cached)\n")
	      (cdr tab-val))
	    (begin ; otherwise, compute bands and cache results
	      (set! num-bands b)
	      (set! k-points (list (vector3-scale k kdir1)))
	      (run-parity p false)
	      (let ((v (compute-group-velocity-component kdir1)))
		; cache computed values:
		(map (lambda (b f v)
		       (let ((tabval (assoc 
				      (cons b (vector-ref k0s (- b band-min)))
				      bktab)))
			 (if (or (not tabval)
				 (< (abs (- f omega)) (abs (cadr tabval))))
			     (vector-set! k0s (- b band-min) k))) ; cache k0
		       (set! bktab (cons (cons (cons b k) (cons (- f omega) v))
					 bktab)))
		     (arith-sequence band-min 1 (- b band-min -1)) 
		     (ncdr (- band-min 1) freqs)
		     (ncdr (- band-min 1) v))
		; finally return (frequency - omega . derivative):
		(let ((fun (- (car (reverse freqs)) omega)))
		  (print "find-k " b " at " k ": " fun "\n")
		  (cons fun (car (reverse v)))))))))
    (randomize-fields) ; don't let previous computations interfere
    (let ((ks (reverse (map
			(lambda (b)
			  (find-root-deriv (rootfun b) tol kmag-min kmag-max
					   (vector-ref k0s (- b band-min))))
			(arith-sequence band-max -1 nb)))))
      (if band-funcs
	  (map (lambda (b k)
		 (set! num-bands b)
		 (set! k-points (list (vector3-scale k kdir1)))
		 (run-parity p false
			     (lambda (b')
			       (if (= b' b)
				   (map (lambda (f) 
					  (apply-band-func-thunk f b true))
					band-funcs)))))
	       (arith-sequence band-max -1 nb) (reverse ks)))
      (set! num-bands num-bands-save)
      (set! k-points k-points-save)
      (print parity "kvals:, " omega ", " band-min ", " band-max)
      (vector-map (lambda (k) (print ", " k)) kdir1)
      (map (lambda (k) (print ", " k)) ks) 
      (print "\n")
      ks)))

; ****************************************************************

(define (sqmatrix-diag m)
  (map (lambda (i) (sqmatrix-ref m i i))
       (arith-sequence 0 1 (sqmatrix-size m))))

(define (fix-phase-consistency old-eigs first-band)
  (let ((dots (dot-eigenvectors old-eigs first-band)))
    (let ((phases (map (lambda (d) (conj (make-polar 1 (angle d))))
		       (sqmatrix-diag dots))))
      (map (lambda (i phase) 
	     (scale-eigenvector i phase)
	     (conj phase))
	   (arith-sequence first-band 1 (length phases)) phases))))

; ****************************************************************

; Load GNU Readline support, for easier command-line editing support.
; This is not loaded in by default in Guile 1.3.2+ because readline
; is licensed under the GPL, which would have caused Guile to effectively
; be under the GPL itself.  However, since the MIT Photonic Bands package
; is under the GPL too, we can load Readline by default with no problems.

@ACTIVATE_READLINE@  ; command to activate readline is determined by configure

(set! scm-repl-prompt "mpb> ")

; ****************************************************************