diff --git a/python/geom.py b/python/geom.py index a77c31db6..5070f04cd 100755 --- a/python/geom.py +++ b/python/geom.py @@ -306,6 +306,24 @@ def eval_susceptibility(self,freq): else: return self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) * sigma +class IIR_Susceptibility(Susceptibility): + def __init__(self, num, den, **kwargs): + super(IIR_Susceptibility, self).__init__(**kwargs) + self.num = num + self.den = den + + def eval_susceptibility(self,freq): + sigma = np.expand_dims(Matrix(diag=self.sigma_diag,offdiag=self.sigma_offdiag),axis=0) + + freq = np.squeeze(freq) + eval_num = np.zeros(freq.shape,dtype=np.complex128) + eval_den = np.zeros(freq.shape,dtype=np.complex128) + N = len(self.num) + D = len(self.den) + for k in range(N): eval_num += pow(1j*freq,N-k)*self.num[k] + for k in range(D): eval_den += pow(1j*freq,D-k)*self.den[k] + + return sigma * eval_num[:,np.newaxis,np.newaxis] / eval_den[:,np.newaxis,np.newaxis] class DrudeSusceptibility(Susceptibility): diff --git a/python/meep.i b/python/meep.i index 6222c0406..ba8a10530 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1169,6 +1169,12 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, } } +//-------------------------------------------------- +// IIR transfer function typemaps +//-------------------------------------------------- +%apply (double *INPLACE_ARRAY1, int DIM1) {(double *r, int N)}; +%apply (double *INPLACE_ARRAY1, int DIM1) {(double *vec_numS, int N), (double *vec_denS, int D), (double *vec_numZ, int N_z), (double *vec_denZ, int D_z)}; + // For some reason SWIG needs the namespaced version too %apply material_type_list { meep_geom::material_type_list }; @@ -1398,6 +1404,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where GyrotropicDrudeSusceptibility, GyrotropicLorentzianSusceptibility, GyrotropicSaturatedSusceptibility, + IIR_Susceptibility, Lattice, LorentzianSusceptibility, Matrix, diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 00e9e07df..f67d5bafb 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -336,6 +336,9 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru s->saturated_gyrotropy = false; s->transitions.resize(0); s->initial_populations.resize(0); + s->iir = false; + s->numS.resize(0); + s->denS.resize(0); if (PyObject_HasAttrString(po, "frequency")) { if (!get_attr_dbl(po, &s->frequency, "frequency")) { return 0; } @@ -390,6 +393,33 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru s->drude = false; } + // ---------- IIR filter sus. typemaps --------------------// + if (PyObject_HasAttrString(po, "num")) { + s->iir = true; + + // numerator python list to C++ vector + PyObject *py_pop = PyObject_GetAttrString(po, "num"); + if (!py_pop) { return 0; } + int length = PyList_Size(py_pop); + s->numS.resize(length); + + for (int i = 0; i < length; ++i) { + s->numS[i] = PyFloat_AsDouble(PyList_GetItem(py_pop, i)); + } + + // denominator python list to C++ vector + py_pop = PyObject_GetAttrString(po, "den"); + if (!py_pop) { return 0; } + length = PyList_Size(py_pop); + s->denS.resize(length); + + for (int i = 0; i < length; ++i) { + s->denS[i] = PyFloat_AsDouble(PyList_GetItem(py_pop, i)); + } + Py_DECREF(py_pop); + + } + s->is_file = false; return 1; diff --git a/src/Makefile.am b/src/Makefile.am index ad571b2d2..afbbb4ef3 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -12,7 +12,7 @@ bicgstab.hpp meepgeom.hpp material_data.hpp libmeep_la_SOURCES = array_slice.cpp anisotropic_averaging.cpp \ bands.cpp boundaries.cpp bicgstab.cpp casimir.cpp \ control_c.cpp cw_fields.cpp dft.cpp dft_ldos.cpp energy_and_flux.cpp \ -fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp \ +fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp iir-susceptibility.cpp \ initialize.cpp integrate.cpp integrate2.cpp monitor.cpp mympi.cpp \ multilevel-atom.cpp near2far.cpp output_directory.cpp random.cpp \ sources.cpp step.cpp step_db.cpp stress.cpp structure.cpp structure_dump.cpp \ diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp new file mode 100644 index 000000000..3638c63dd --- /dev/null +++ b/src/iir-susceptibility.cpp @@ -0,0 +1,417 @@ +/* Copyright (C) 2005-2019 Massachusetts Institute of Technology. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +#include +#include +#include "meep.hpp" +#include "meep_internals.hpp" +#include "material_data.hpp" + +using namespace std; + +namespace meep { + +// --------------------------------------------------------------- // +// S to Z transform routines +// --------------------------------------------------------------- // +/* The following routines are use to convert the IIR transfer function + * in the s domain to the z domain. + */ + + // generalized way to convert a polynomial in the s domain to the z +// domain. Based on the paper: +// Scott, Dan M. "A simplified method for the bilinear sz transformation." IEEE Transactions on Education 37.3 (1994): 289-292. +// https://digitalcommons.unl.edu/cgi/viewcontent.cgi?referer=https://www.google.com/&httpsredir=1&article=1085&context=imsefacpub +void flip_array(realnum* arr, int N){ + realnum temp; + for (int i = 0; i < N/2; ++i) { + temp = arr[N-i-1]; + arr[N-i-1] = arr[i]; + arr[i] = temp; + } +} + +void polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma) { + flip_array(r,N); + + int k, j, n=N-1; + double sum, beta_powers[n], gamma_powers[n]; + master_printf("entered PT -------------- \n"); + beta_powers[0] = gamma_powers[0] = 1; + for (j=1;j<=n;j++) + { + beta_powers[j] = beta*beta_powers[j- 1]; + gamma_powers[j] = gamma*gamma_powers[j- 1]; + } + + for (k=0; k 0; k--) factorial *= k; + return factorial; +} + +realnum comb(int n, int k) { + return factorial(n) / (factorial(k) * factorial(n-k)); +} + +// Uses Tustins method to discretize the transfer function. S = (2/T)*(1-z^-1)/(1+z^-1) +// https://github.com/scipy/scipy/blob/v1.3.0/scipy/signal/filter_design.py#L1965-L2043 +void tustins_method( + realnum* vec_numS, + int N, + realnum* vec_denS, + int D, + realnum* vec_numZ, + int N_z, + realnum* vec_denZ, + int D_z, + realnum T + ) { + (void)N_z; + (void)D_z; + + int M = D-1; //Denominator should always have highest rank + int Np = M; + int Dp = M; + + std::vector bprime(Np + 1, 0.0); + std::vector aprime(Dp + 1, 0.0); + + int j, i, k, l; + realnum val; + for (j=0;j<=Np;j++) { + val = 0.0; + for(i=0;isz_data = sz; + + return (void *)d; +} + +void iir_susceptibility::init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *data) const { + (void)dt; // unused + (void)gv; + iir_data *d = (iir_data *)data; + + // initialize everything to 0 + size_t sz_data = d->sz_data; + memset(d, 0, sz_data); + d->sz_data = sz_data; + size_t ntot = d->ntot = gv.ntot(); + + // Take that contiguous block of data passed by the user and + // distribute it to the appropriate array locations + realnum *P = d->data; + + FOR_COMPONENTS(c) DOCMP2 { + if (needs_P(c, cmp, W)) { + d->P[c][cmp] = P; + P += ntot; + d->storage[c][cmp] = P; + P += ntot * (D_z-1); + } + } +} + +void *iir_susceptibility::copy_internal_data(void *data) const { + iir_data *d = (iir_data *)data; + if (!d) return 0; + iir_data *dnew = (iir_data *)malloc(d->sz_data); + memcpy(dnew, d, d->sz_data); + + dnew->sz_data = d->sz_data; + dnew->ntot = d->ntot; + realnum *P = dnew->data; + + FOR_COMPONENTS(c) DOCMP2 { + if (d->P[c][cmp]) { + dnew->P[c][cmp] = P; + P += ntot; + dnew->storage[c][cmp] = P; + P += ntot * (D_z-1); + } + } + return (void *)dnew; +} + +#define SWAP(t, a, b) \ + { \ + t SWAP_temp = a; \ + a = b; \ + b = SWAP_temp; \ + } + +// stable averaging of offdiagonal components +#define OFFDIAG(u, g, sx, s) \ + (0.25 * ((g[i] + g[i - sx]) * u[i] + (g[i + s] + g[(i + s) - sx]) * u[i + s])) + +void iir_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], + realnum *W_prev[NUM_FIELD_COMPONENTS][2], double dt, + const grid_volume &gv, void *P_internal_data) const { + iir_data *d = (iir_data *)P_internal_data; + size_t ntot = d->ntot; + (void)dt; + (void)W_prev; // unused; + + FOR_COMPONENTS(c) DOCMP2 { + if (d->P[c][cmp]) { + const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)]; + if (w && s) { + realnum *p = d->P[c][cmp], *storage = d->storage[c][cmp]; + + // directions/strides for offdiagonal terms, similar to update_eh + const direction d = component_direction(c); + const ptrdiff_t is = gv.stride(d) * (is_magnetic(c) ? -1 : +1); + direction d1 = cycle_direction(gv.dim, d, 1); + component c1 = direction_component(c, d1); + ptrdiff_t is1 = gv.stride(d1) * (is_magnetic(c) ? -1 : +1); + const realnum *w1 = W[c1][cmp]; + const realnum *s1 = w1 ? sigma[c][d1] : NULL; + direction d2 = cycle_direction(gv.dim, d, 2); + component c2 = direction_component(c, d2); + ptrdiff_t is2 = gv.stride(d2) * (is_magnetic(c) ? -1 : +1); + const realnum *w2 = W[c2][cmp]; + const realnum *s2 = w2 ? sigma[c][d2] : NULL; + + LOOP_OVER_VOL_OWNED(gv, c, i) { + if (s[i] != 0) { + realnum w_cur; + if (s2 && !s1) { // make s1 the non-NULL one if possible + SWAP(direction, d1, d2); + SWAP(component, c1, c2); + SWAP(ptrdiff_t, is1, is2); + SWAP(const realnum *, w1, w2); + SWAP(const realnum *, s1, s2); + } + if (s1 && s2) { // 3x3 anisotropic + w_cur = s[i] * w[i] + OFFDIAG(s1, w1, is1, is) + OFFDIAG(s2, w2, is2, is); + }else if (s1) { // 2x2 + w_cur = s[i] * w[i] + OFFDIAG(s1, w1, is1, is); + }else{ // isotropic + w_cur = s[i] * w[i]; + } + + // Implement a "Direct Form II Transposed" difference equation + //https://ocw.mit.edu/courses/mechanical-engineering/2-161-signal-processing-continuous-and-discrete-fall-2008/lecture-notes/lecture_20.pdf + p[i] = storage[0*ntot + i] + numZ[0] * w_cur; + for (int J = 0; J < D_z-2; J++) { + storage[(J)*ntot + i] = storage[(J+1)*ntot + i] + numZ[J+1] * w_cur - denZ[J+1] * p[i]; + } + storage[(D_z-2)*ntot + i] = numZ[D_z-1] * w_cur - denZ[N_z-1] * p[i]; + } + } + } + } + } +} + +void iir_susceptibility::subtract_P(field_type ft, + realnum *f_minus_p[NUM_FIELD_COMPONENTS][2], + void *P_internal_data) const { + iir_data *d = (iir_data *)P_internal_data; + field_type ft2 = ft == E_stuff ? D_stuff : B_stuff; // for sources etc. + size_t ntot = d->ntot; + FOR_FT_COMPONENTS(ft, ec) DOCMP2 { + if (d->P[ec][cmp]) { + component dc = field_type_component(ft2, ec); + if (f_minus_p[dc][cmp]) { + realnum *p = d->P[ec][cmp]; + realnum *fmp = f_minus_p[dc][cmp]; + for (size_t i = 0; i < ntot; ++i) + fmp[i] -= p[i]; + } + } + } +} + +// constructor must take the s domain coefficients and transform +// them to the z domain and store them in the class. +iir_susceptibility::iir_susceptibility(std::vector vec_numS, std::vector vec_denS, double dt): + N (vec_numS.size()), D (vec_denS.size()), N_z (vec_denS.size()), D_z (vec_denS.size()), T (dt), + numS (new realnum[N]), denS (new realnum[D]), + numZ (new realnum[N_z]), denZ (new realnum[D_z]) + { + //if (vec_numS.size() > vec_denS.size()) abort("error: for iir-filter susceptibilities, the numerator rank cannot be higher than the denominator rank.\n"); + + // copy over the S domain numerator and denominator vectors to arrays + for (int i=0; iP[c][0][0] ? 1 : 0; +} + +realnum *iir_susceptibility::cinternal_notowned_ptr(int inotowned, component c, int cmp, + int n, void *P_internal_data) const { + iir_data *d = (iir_data *)P_internal_data; + (void)inotowned; // always = 0 + if (!d || !d->P[c][cmp]) return NULL; + return d->P[c][cmp] + n; +} + +std::complex iir_susceptibility::chi1(double freq, double sigma) { + std::complex num = 0; + std::complex den = 0; + for (int i=0;i(0,freq),N-i) * numS[i]; + for (int i=0;i(0,freq),D-i) * denS[i]; + + return sigma * num / den; +} + +void iir_susceptibility::dump_params(h5file *h5f, size_t *start) { + size_t num_params = N_z + D_z; + size_t params_dims[1] = {num_params}; + double params_data[num_params]; + for (int i=0; iwrite_chunk(1, start, params_dims, params_data); + *start += num_params; +} + +// be sure to clean up the numerator and denominator arrays +iir_susceptibility::~iir_susceptibility() { + delete[] numZ; + delete[] denZ; + delete[] numS; + delete[] denS; +} + +} \ No newline at end of file diff --git a/src/material_data.hpp b/src/material_data.hpp index 20c6f9823..8f456f0bc 100644 --- a/src/material_data.hpp +++ b/src/material_data.hpp @@ -69,6 +69,9 @@ typedef struct susceptibility_struct { bool is_file; std::vector transitions; std::vector initial_populations; + // IIR filter + bool iir; + std::vector numS, denS; } susceptibility; struct susceptibility_list { diff --git a/src/meep.hpp b/src/meep.hpp index 9fa946848..dd4f07ca0 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -370,6 +370,58 @@ class multilevel_susceptibility : public susceptibility { realnum *sigmat; // 5*T transition-specific sigma-diagonal factors }; +/* a IIR susceptibility + Susceptibility modeled by a general rational function (ratios of p(w)/q(w)). + The user provides the numerator and denominator coefficients in the continuous (s) + domain. One of the conversion routines is used to transform to the discrete (z) + domain. */ + +class iir_susceptibility : public susceptibility { +public: + iir_susceptibility(std::vector vec_numS, std::vector vec_denS, double dt); + iir_susceptibility(const iir_susceptibility &from); + + virtual susceptibility *clone() const { return new iir_susceptibility(*this); } + + virtual ~iir_susceptibility(); + + virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2], + realnum *W_prev[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, + void *P_internal_data) const; + + virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2], + void *P_internal_data) const; + + virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const; + + virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], double dt, + const grid_volume &gv, void *data) const; + + virtual void *copy_internal_data(void *data) const; + + virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const; + + virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp, int n, + void *P_internal_data) const; + + virtual void dump_params(h5file *h5f, size_t *start); + virtual int get_num_params() { return 4; } + + virtual std::complex chi1(double freq, double sigma = 1); + +protected: +int N, D, N_z, D_z; +realnum T, *numS, *denS, *numZ, *denZ; + + +}; +// tustins method (aka trapezoidal rule) to convert from s domain to z domain. +void polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma); +void backward_difference(realnum* vec_numS, int N, realnum* vec_denS, int D, realnum* vec_numZ, int N_z, realnum* vec_denZ, int D_z, realnum T); +void tustins_method(realnum* vec_numS, int N, realnum* vec_denS, int D, realnum* vec_numZ, int N_z, realnum* vec_denZ, int D_z, realnum T); +realnum factorial(int i); +realnum comb(int n, int k); + class grace; // h5file.cpp: HDF5 file I/O. Most users, if they use this diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index ad4b41ce0..25b78c1e6 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1397,21 +1397,24 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s) : ss->drude ? meep::GYROTROPIC_DRUDE : meep::GYROTROPIC_LORENTZIAN; sus = new meep::gyrotropic_susceptibility(meep::vec(ss->bias.x, ss->bias.y, ss->bias.z), ss->frequency, ss->gamma, ss->alpha, model); - } - else { + } else if (ss->iir) { + // Initialize IIR filter sus + sus = new meep::iir_susceptibility(ss->numS,ss->denS,s->dt); + } else { sus = new meep::lorentzian_susceptibility(ss->frequency, ss->gamma, ss->drude); } if (!meep::quiet) { master_printf("%s%s susceptibility: frequency=%g, gamma=%g", noisy ? "noisy " : gyrotropic ? "gyrotropic " : "", ss->saturated_gyrotropy ? "Landau-Lifshitz-Gilbert-type" - : ss->drude ? "drude" : "lorentzian", - ss->frequency, ss->gamma); + : ss->iir ? "IIR filter" + : ss->drude ? "drude" : "lorentzian", ss->frequency, ss->gamma); if (noisy) master_printf(", amp=%g ", ss->noise_amp); if (gyrotropic) { if (ss->saturated_gyrotropy) master_printf(", alpha=%g", ss->alpha); master_printf(", bias=(%g,%g,%g)", ss->bias.x, ss->bias.y, ss->bias.z); } + if (ss->iir) master_printf(", numerator=%lu, denominator=%lu", ss->numS.size(), ss->denS.size()); master_printf("\n"); } } diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index 80a88543f..6033a8600 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -35,6 +35,10 @@ using namespace std; namespace meep { +// ---------------------------------------------------------- // +// Base susceptibility class +// ---------------------------------------------------------- // + int susceptibility::cur_id = 0; susceptibility *susceptibility::clone() const { @@ -95,6 +99,10 @@ bool susceptibility::needs_W_notowned(component c, realnum *W[NUM_FIELD_COMPONEN return false; } +// ---------------------------------------------------------- // +// Lorentzian susceptibility +// ---------------------------------------------------------- // + typedef struct { size_t sz_data; size_t ntot; @@ -313,6 +321,10 @@ void lorentzian_susceptibility::dump_params(h5file *h5f, size_t *start) { *start += num_params; } +// ---------------------------------------------------------- // +// Noisy lorentzian +// ---------------------------------------------------------- // + void noisy_lorentzian_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], realnum *W_prev[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *P_internal_data) const { @@ -346,6 +358,10 @@ void noisy_lorentzian_susceptibility::dump_params(h5file *h5f, size_t *start) { *start += num_params; } +// ---------------------------------------------------------- // +// Gyrotropic susceptibility +// ---------------------------------------------------------- // + gyrotropic_susceptibility::gyrotropic_susceptibility(const vec &bias, double omega_0, double gamma, double alpha, gyrotropy_model model) : omega_0(omega_0), gamma(gamma), alpha(alpha), model(model) {