From d8cf22d223fa5df9530a4ac38efeb4edfb4164b7 Mon Sep 17 00:00:00 2001 From: Christopher Hogan Date: Tue, 28 Aug 2018 10:15:50 -0500 Subject: [PATCH] Add guile API for get-eigenmode-coefficients (#477) * Add get-eigenmode-coefficients to scheme API * Add support for kpoint functions * Return kpoints from get-eigenmode-coefficients * Remove extra define --- scheme/examples/mode-coeffs.ctl | 99 +++++++++++++++++++++++++++++++++ scheme/meep-ctl-swig.hpp | 12 +++- scheme/meep.cpp | 18 +++++- scheme/meep.i | 86 ++++++++++++++++++++++++---- scheme/meep.scm.in | 46 +++++++++++++++ 5 files changed, 249 insertions(+), 12 deletions(-) create mode 100644 scheme/examples/mode-coeffs.ctl diff --git a/scheme/examples/mode-coeffs.ctl b/scheme/examples/mode-coeffs.ctl new file mode 100644 index 000000000..faf30ea41 --- /dev/null +++ b/scheme/examples/mode-coeffs.ctl @@ -0,0 +1,99 @@ + +(set-param! resolution 15) + +(define-param w 1) ; width of waveguide +(define-param L 10) ; length of waveguide + +(define-param dair 3.0) +(define-param dpml 3.0) + +;; mode frequency +(define-param fcen 0.20) ; > 0.5/sqrt(11) to have at least 2 modes + +(define (run-test mode-num kpoint-func) + (reset-meep) + + (let* ((sx (+ dpml L dpml)) + (sy (+ dpml dair w dair dpml)) + (prism_x (+ sx 1)) + (prism_y (/ w 2)) + (verts + (list + (vector3 (* prism_x -1) prism_y) + (vector3 prism_x prism_y) + (vector3 prism_x (* prism_y -1)) + (vector3 (* prism_x -1) (* prism_y -1)))) + (modes-to-check '(1 2)) ; indices of modes for which to compute expansion coefficients + (xm (- (* 0.5 sx) dpml))) ; x-coordinate of monitor + + (set! geometry-lattice (make lattice (size sx sy no-size))) + (set! geometry + (list + (make prism + (center (vector3 0 0)) + (vertices verts) + (height infinity) + (material (make medium (epsilon 12.0)))))) + + (set! pml-layers (list (make pml (thickness dpml)))) + + + (set! sources (list + (make eigenmode-source + (src (make gaussian-src (frequency fcen) (fwidth (* 0.5 fcen)))) + (center (vector3 (+ (* -0.5 sx) dpml) 0)) + (component ALL-COMPONENTS) + (size (vector3 0 (- sy (* 2 dpml)))) + (eig-match-freq? true) + (eig-band mode-num) + (eig-resolution 32)))) + + (set! symmetries + (list + (make mirror-sym (direction Y) (phase (if (odd? mode-num) 1 -1))))) + + (let ((mflux (add-mode-monitor fcen 0 1 + (make mode-region (center (vector3 xm 0)) (size (vector3 0 (- sy (* 2 dpml))))))) + (mode-flux (add-flux fcen 0 1 + (make flux-region (center (vector3 xm 0)) (size (vector3 0 (- sy (* 2 dpml)))))))) + + (run-sources+ 100) + + (let* ((alpha-vgrp-kpoints (get-eigenmode-coefficients mflux modes-to-check #:kpoint-func kpoint-func)) + (alpha (first alpha-vgrp-kpoints)) + (vgrp (second alpha-vgrp-kpoints)) + (kpoints (third alpha-vgrp-kpoints)) + (mode-power (car (get-fluxes mode-flux))) + (test-passed #t) + (tolerance 5.0e-3) + (c0 (array-ref alpha (- mode-num 1) 0 0))) ; coefficient of forward-traveling wave for mode # mode_num + + (if (or (not (vector3-close? (first kpoints) (vector3 0.604301 0 0) 1e-7)) + (not (vector3-close? (second kpoints) (vector3 0.494353 0 0) 1e-2))) + (begin + (print "Test failed") + (exit 1))) + + (map (lambda (nm) + (if (not (= mode-num nm)) + (let ((cfrel (/ (magnitude (array-ref alpha (- nm 1) 0 0)) (magnitude c0))) + (cbrel (/ (magnitude (array-ref alpha (- nm 1) 0 1)) (magnitude c0)))) + (if (or (> cfrel tolerance) (> cbrel tolerance)) + (set! test-passed #f))))) + modes-to-check) + + ;; test 1: coefficient of excited mode >> coeffs of all other modes + (if (not test-passed) + (begin + (print "Test failed") + (exit 1))) + + ;; test 2: |mode coeff|^2 = power + (if (> (abs (- 1 (/ mode-power (* (magnitude c0) (magnitude c0))))) 0.1) + (begin + (print "Test failed") + (exit 1))))))) + +(run-test 1 '()) +(run-test 2 '()) +(run-test 1 (lambda (freq mode) (vector3 0 0))) diff --git a/scheme/meep-ctl-swig.hpp b/scheme/meep-ctl-swig.hpp index 7248621c7..a6e3f486e 100644 --- a/scheme/meep-ctl-swig.hpp +++ b/scheme/meep-ctl-swig.hpp @@ -5,6 +5,11 @@ #ifndef MEEP_CTL_SWIG_HPP #define MEEP_CTL_SWIG_HPP 1 +struct kpoint_list { + meep::vec *kpoints; + size_t n; +}; + vector3 vec_to_vector3(const meep::vec &); meep::vec vector3_to_vec(const vector3 v3); void set_dimensions(int dims); @@ -23,12 +28,17 @@ meep::structure *make_structure(int dims, vector3 size, vector3 center, double global_D_conductivity_diag_, double global_B_conductivity_diag_); -ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, +ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, double fmin, double fmax, int maxbands, double spectral_density, double Q_thresh, double rel_err_thresh, double err_thresh, double rel_amp_thresh, double amp_thresh); +kpoint_list do_get_eigenmode_coefficients(meep::fields *f, meep::dft_flux flux, const meep::volume &eig_vol, + int *bands, int num_bands, int parity, std::complex *coeffs, + double *vgrp, double eig_resolution, double eigensolver_tol, + meep::kpoint_func user_kpoint_func, void *user_kpoint_data); + ctlio::number_list dft_flux_flux(meep::dft_flux *f); ctlio::number_list dft_force_force(meep::dft_force *f); ctlio::number_list dft_ldos_ldos(meep::dft_ldos *f); diff --git a/scheme/meep.cpp b/scheme/meep.cpp index f3af44fad..e9a0ad8a6 100644 --- a/scheme/meep.cpp +++ b/scheme/meep.cpp @@ -33,7 +33,7 @@ void ctl_export_hook(void) /**************************************************************************/ -ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, +ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, double fmin, double fmax, int maxbands, double spectral_density, double Q_thresh, double rel_err_thresh, double err_thresh, @@ -66,6 +66,22 @@ ctlio::cvector3_list do_harminv(ctlio::cnumber_list vals, double dt, return res; } +kpoint_list do_get_eigenmode_coefficients(fields *f, dft_flux flux, const volume &eig_vol, int *bands, int num_bands, + int parity, std::complex *coeffs, double *vgrp, double eig_resolution, + double eigensolver_tol, meep::kpoint_func user_kpoint_func, + void *user_kpoint_data) +{ + size_t num_kpoints = num_bands * flux.Nfreq; + meep::vec *kpoints = new meep::vec[num_kpoints]; + + f->get_eigenmode_coefficients(flux, eig_vol, bands, num_bands, parity, eig_resolution, eigensolver_tol, + coeffs, vgrp, user_kpoint_func, user_kpoint_data, kpoints); + + kpoint_list res = {kpoints, num_kpoints}; + + return res; +} + /**************************************************************************/ /* This is a wrapper function to fool SWIG...since our list constructor diff --git a/scheme/meep.i b/scheme/meep.i index de41b1bd4..c0d8c5c0b 100644 --- a/scheme/meep.i +++ b/scheme/meep.i @@ -60,22 +60,30 @@ static inline std::complex my_complex_func3(std::complex x) { return std::complex(cret.re, cret.im); } +static meep::vec my_kpoint_func(double freq, int mode, void *user_data) { + SCM scm_res = gh_call2((SCM)user_data, ctl_convert_number_to_scm(freq), + ctl_convert_number_to_scm(mode)); + vector3 v3 = ctl_convert_vector3_to_c(scm_res); + meep::vec result(v3.x, v3.y, v3.z); + return result; +} + %} %typecheck(SWIG_TYPECHECK_COMPLEX) std::complex { $1 = SwigComplex_Check($input); } -%typemap(guile,out) complex, std::complex, std::complex { +%typemap(out) complex, std::complex, std::complex { $result = scm_make_rectangular(ctl_convert_number_to_scm($1.real()), ctl_convert_number_to_scm($1.imag())); } -%typemap(guile,in) complex, std::complex, std::complex { +%typemap(in) complex, std::complex, std::complex { cnumber cnum = ctl_convert_cnumber_to_c($input); $1 = std::complex(cnum.re, cnum.im); } -%typemap(guile,in) std::complex(*)(meep::vec const &) { +%typemap(in) std::complex(*)(meep::vec const &) { my_complex_func_scm = $input; $1 = my_complex_func; } @@ -83,7 +91,7 @@ static inline std::complex my_complex_func3(std::complex x) { $1 = SCM_NFALSEP(scm_procedure_p($input)); } -%typemap(guile,in) std::complex(*)(std::complex) { +%typemap(in) std::complex(*)(std::complex) { my_complex_func3_scm = $input; $1 = my_complex_func3; } @@ -91,7 +99,7 @@ static inline std::complex my_complex_func3(std::complex x) { $1 = SCM_NFALSEP(scm_procedure_p($input)); } -%typemap(guile,in) (std::complex (*func)(double t, void *), void *data) { +%typemap(in) (std::complex (*func)(double t, void *), void *data) { $1 = my_complex_func2; $2 = (void *) $input; // input is SCM pointer to Scheme function } @@ -99,13 +107,13 @@ static inline std::complex my_complex_func3(std::complex x) { $1 = SCM_NFALSEP(scm_procedure_p($input)); } -%typemap(guile,in) meep::vec { +%typemap(in) meep::vec { $1 = vector3_to_vec(ctl_convert_vector3_to_c($input)); } -%typemap(guile,out) meep::vec { +%typemap(out) meep::vec { $result = ctl_convert_vector3_to_scm(vec_to_vector3($1)); } -%typemap(guile,in) meep::vec const & %{ +%typemap(in) meep::vec const & %{ meep::vec vec__$1 = vector3_to_vec(ctl_convert_vector3_to_c($input)); $1 = &vec__$1; %} @@ -115,7 +123,7 @@ static inline std::complex my_complex_func3(std::complex x) { /* field_function arguments are passed as a cons pair of (components . func) in order to set all four arguments at once. */ -%typemap(guile,in) (int num_fields, const meep::component *components, meep::field_function fun, void *fun_data_) (my_field_func_data data) { +%typemap(in) (int num_fields, const meep::component *components, meep::field_function fun, void *fun_data_) (my_field_func_data data) { $1 = list_length(gh_car($input)); $2 = new meep::component[$1]; for (int i = 0; i < $1; ++i) @@ -137,7 +145,7 @@ static inline std::complex my_complex_func3(std::complex x) { /* integrate2 arguments are passed as a cons pair of ((components1 . components2) . func) in order to set all six arguments at once. */ -%typemap(guile,in) (int num_fields1, const meep::component *components1, int num_fields2, const meep::component *components2, meep::field_function integrand, void *integrand_data_) (my_field_func_data data) { +%typemap(in) (int num_fields1, const meep::component *components1, int num_fields2, const meep::component *components2, meep::field_function integrand, void *integrand_data_) (my_field_func_data data) { $1 = list_length(gh_car(gh_car($input))); $2 = new meep::component[$1]; for (int i = 0; i < $1; ++i) @@ -182,6 +190,64 @@ static inline std::complex my_complex_func3(std::complex x) { } } +// do_get_eigenmode_coefficients + +%typemap(in) (int *bands, int num_bands) { + $2 = list_length($input); + $1 = new int[$2]; + + for (int i = 0; i < $2; ++i) { + $1[i] = integer_list_ref($input, i); + } +} + +%typemap(freearg) (int *bands, int num_bands) { + if ($1) { + delete[] $1; + } +} + +%typemap(in) (meep::kpoint_func user_kpoint_func, void *user_kpoint_data) { + if (SCM_NFALSEP(scm_procedure_p($input))) { + $1 = my_kpoint_func; + $2 = (void*)$input; + } + else { + $1 = NULL; + $2 = NULL; + } +} + +%typemap(in, noblock=1) std::complex *coeffs { + scm_t_array_handle coeffs_handle; + scm_array_get_handle($input, &coeffs_handle); + $1 = (std::complex*)scm_array_handle_uniform_writable_elements(&coeffs_handle); +} + +%typemap(in, noblock=1) double *vgrp { + scm_t_array_handle vgrp_handle; + scm_array_get_handle($input, &vgrp_handle); + $1 = (double*)scm_array_handle_uniform_writable_elements(&vgrp_handle); +} + +%typemap(freearg) std::complex *coeffs { + scm_array_handle_release(&coeffs_handle); +} + +%typemap(freearg) double *vgrp { + scm_array_handle_release(&vgrp_handle); +} + +%typemap(out) kpoint_list { + int i; + list cur_list = SCM_EOL; + for (i = $1.n - 1; i >= 0; --i) { + cur_list = gh_cons(vector32scm(vec_to_vector3($1.kpoints[i])), cur_list); + } + $result = cur_list; + delete[] $1.kpoints; +} + // Need to tell SWIG about any method that returns a new object // which needs to be garbage-collected. %newobject meep::fields::open_h5file; diff --git a/scheme/meep.scm.in b/scheme/meep.scm.in index aef640a01..56268e803 100644 --- a/scheme/meep.scm.in +++ b/scheme/meep.scm.in @@ -579,6 +579,26 @@ (load-flux fname flux) (meep-dft-flux-scale-dfts flux -1.0)) +; **************************************************************** +; Mode monitor + +(define mode-region flux-region) + +(define (fields-add-mode-monitor fields fcen df nfreq . fluxes) + (if (not (= (length fluxes) 1)) (error "add-mode-monitor expected just one mode-region.")) + (let* ((region (car fluxes)) + (v (volume (center (object-property-value region 'center)) + (size (object-property-value region 'size)))) + (d0 (object-property-value region 'direction)) + (d (if (< d0 0) + (meep-fields-normal-direction fields v) + d0))) + (meep-fields-add-mode-monitor fields d v (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq))) + +(define (add-mode-monitor fcen df nfreq . fluxes) + (if (null? fields) (init-fields)) + (apply fields-add-mode-monitor (append (list fields fcen df nfreq) fluxes))) + ; **************************************************************** ; Force spectra (from stress tensor) - very similar interface to flux spectra @@ -1084,6 +1104,32 @@ (define (harminv c pt fcen df . mxbands) (harminv! harminv-data harminv-data-dt harminv-results c pt fcen df mxbands)) + +; **************************************************************** +; get-eigenmode-coefficients + +(define (get-keyword-value args keyword default) + (let ((kv (memq keyword args))) + (if (and kv (>= (length kv) 2)) + (cadr kv) + default))) + +(define (get-eigenmode-coefficients flux bands . args) + (if (null? fields) + (error "init-fields is required before using get-eigenmode-coefficients")) + (let ((eig-parity (get-keyword-value args #:eig-parity NO-PARITY)) + (eig-vol (get-keyword-value args #:eig-vol (meep-dft-flux-where-get flux))) + (eig-resolution (get-keyword-value args #:eig-resolution 0)) + (eig-tolerance (get-keyword-value args #:eig-tolerance 1e-7)) + (num-bands (length bands)) + (kpoint-func (get-keyword-value args #:kpoint-func '()))) + (begin + (define coeffs (make-typed-array 'c64 '0 num-bands (meep-dft-flux-Nfreq-get flux) 2)) + (define vgrp (make-typed-array 'f64 '0 num-bands (meep-dft-flux-Nfreq-get flux))) + (define kpoints (do-get-eigenmode-coefficients fields flux eig-vol bands eig-parity coeffs vgrp + eig-resolution eig-tolerance kpoint-func)) + (list coeffs vgrp kpoints)))) + ; **************************************************************** ; dft-ldos step function