From b37d1b067174f85e86f0ac8fcde6796ec1f3eb5c Mon Sep 17 00:00:00 2001 From: smartalecH Date: Wed, 12 Jun 2019 09:45:35 -0600 Subject: [PATCH 1/8] init --- src/iir-susceptibility.cpp | 35 +++++++++++++++++++++++++++++ src/meep.hpp | 45 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) create mode 100644 src/iir-susceptibility.cpp diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp new file mode 100644 index 000000000..df1471425 --- /dev/null +++ b/src/iir-susceptibility.cpp @@ -0,0 +1,35 @@ +/* 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 + */ + +typedef struct { + size_t sz_data; + size_t ntot; + realnum *P; + realnum *P_prev; + realnum *W_prev; + realnum data[1]; +} iir_data; + +// constructor must take the s domain coefficients and transform +// them to the z domain and store them in the class. +iir_susceptibility::iir_susceptibility(double *num, int num_N, double *den, int den_N, double dt){ + numz_N = num_N; + denz_N = den_N; + tustins_method(num,num_N,den,den_N, dt); + numz = num; + denz = den; +} \ No newline at end of file diff --git a/src/meep.hpp b/src/meep.hpp index 94fd8f03b..a4ee290f2 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -320,6 +320,51 @@ 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: + double *numz, *denz; + int numz_N, denz_N; + + iir_susceptibility(double *num, int num_N, double *den, int den_N, double dt); + + virtual susceptibility *clone() const { return new iir_susceptibility(*this); } + + virtual ~iir_susceptibility() {} + + virtual void update_P(realnum *W[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; } +}; +// generalized method to convert from s domain to z domain. +void polynomial_transform(int N, double alpha, double beta, double delta, double gamma, double *r); +// tustins method (aka trapezoidal rule) to convert from s domain to z domain. +void tustins_method(double *num, int num_N, double *den, int den_N, double T); +// backward difference method used to convert from s domain to z domain. +void backward_difference(double *num, int num_N, double *den, int den_N, double T); + class grace; // h5file.cpp: HDF5 file I/O. Most users, if they use this From 33486dece72a449dd9a20e74e65991a8c23f6d8b Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 12 Jun 2019 16:57:59 -0600 Subject: [PATCH 2/8] more progress --- python/meep.i | 6 +++ src/Makefile.am | 2 +- src/iir-susceptibility.cpp | 107 +++++++++++++++++++++++++++++++++++++ src/meep.hpp | 11 ++-- 4 files changed, 121 insertions(+), 5 deletions(-) diff --git a/python/meep.i b/python/meep.i index 6a0e824a8..aa725265c 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 *num, int num_N), (double *den, int den_N)}; + // For some reason SWIG needs the namespaced version too %apply material_type_list { meep_geom::material_type_list }; 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 index df1471425..365f9ef52 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -15,6 +15,109 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include +#include +#include "meep.hpp" +#include "meep_internals.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 polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma) { + int c, k, j, n=N-1; + double sum, beta_powers[n], gamma_powers[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 Date: Thu, 13 Jun 2019 15:47:31 -0600 Subject: [PATCH 3/8] updatse --- src/iir-susceptibility.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index 365f9ef52..980639f3e 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -53,14 +53,6 @@ void polynomial_transform(double *r, int N, double alpha, double beta, double de for (j=0; j Date: Mon, 24 Jun 2019 10:18:08 -0600 Subject: [PATCH 4/8] more progess --- python/geom.py | 16 +++ python/meep.i | 3 +- python/typemap_utils.cpp | 30 ++++ src/iir-susceptibility.cpp | 274 +++++++++++++++++++++++++++---------- src/material_data.hpp | 3 + src/meep.hpp | 30 ++-- src/meepgeom.cpp | 5 + 7 files changed, 269 insertions(+), 92 deletions(-) diff --git a/python/geom.py b/python/geom.py index 4adb5290e..9a690762a 100755 --- a/python/geom.py +++ b/python/geom.py @@ -302,6 +302,22 @@ 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 = np.squeeze(num).tolist() + self.den = np.squeeze(den).tolist() + + def eval_susceptibility(self,freq): + freq = np.squeeze(freq) + eval_num = np.zeros(freq.shape,dtype=np.float64) + eval_den = np.zeros(freq.shape,dtype=np.float64) + N = self.num.size + D = self.den.size + 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 eval_num / eval_den class DrudeSusceptibility(Susceptibility): diff --git a/python/meep.i b/python/meep.i index a60293807..2f9338642 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1173,7 +1173,7 @@ 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 *num, int num_N), (double *den, int den_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 }; @@ -1395,6 +1395,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 6e3d7029b..5effcbaf7 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -335,6 +335,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/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index 980639f3e..a3a17aac4 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -19,6 +19,7 @@ #include #include "meep.hpp" #include "meep_internals.hpp" +#include "material_data.hpp" using namespace std; @@ -31,104 +32,227 @@ namespace meep { * 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 polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma) { - int c, k, j, n=N-1; - double sum, beta_powers[n], gamma_powers[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; i--) factorial *= i; + return factorial; } -// Uses Tustins method to discretize the transfer function. S = (2/T)*(1-z^-1)/(1+z^-1) -void tustins_method(double *num, int num_N, double *den, int den_N, double dt) { - double c = 2.0 / dt; - double alpha = 2.0, beta = -2.0, delta = dt, gamma = dt; - int k; - - if (den_N < num_N) abort("error: the numerator cannot be higher rank than the denominator"); - if (num_N < den_N) { - double * temp = new double[den_N]; - for (k=0; k < den_N; k++) { - temp[k] = (k <=num_N) ? num[k] : 0.0; - } - delete num; - num = temp; - num_N = den_N; - } - - // Process numerator - polynomial_transform(num,num_N,alpha,beta,delta,gamma); - - // Process denominator - polynomial_transform(den,den_N,alpha,beta,delta,gamma); - - // Normalize - double norm = den[0]; - for (k=0; k < num_N; k++) num[k] = num[k] / norm; - for (k=0; k < den_N; k++) den[k] = den[k] / norm; - for (k=0;k bprime(Np + 1, 0.0); + std::vector aprime(Dp + 1, 0.0); + + int j, i, k, l; + double 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 + iir_data *d = (iir_data *)data; + size_t sz_data = d->sz_data; + memset(d, 0, sz_data); + d->sz_data = sz_data; + size_t ntot = d->ntot = gv.ntot() * Dz; + realnum *P = d->data; + realnum *P_prev = d->data + ntot; + realnum *W_prev = P_prev + ntot; + FOR_COMPONENTS(c) DOCMP2 { + if (needs_P(c, cmp, W)) { + for (int dd = X; dd < R; dd++) { + d->P[c][cmp] = P; + d->P_prev[c][cmp] = P_prev; + d->W_prev[c][cmp] = W_prev; + P += 2 * ntot; + P_prev += 2 * ntot; + } + } + } +} + +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; + (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], *pp = d->P_prev[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; + + 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 + LOOP_OVER_VOL_OWNED(gv, c, i) { + // s[i] != 0 check is a bit of a hack to work around + // some instabilities that occur near the boundaries + // of materials; see PR #666 + if (s[i] != 0) { + realnum pcur = p[i]; + p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + + omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is) + + OFFDIAG(s2, w2, is2, is))); + pp[i] = pcur; + } + } + } else if (s1) { // 2x2 anisotropic + LOOP_OVER_VOL_OWNED(gv, c, i) { + if (s[i] != 0) { // see above + realnum pcur = p[i]; + p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + + omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is))); + pp[i] = pcur; + } + } + } else { // isotropic + LOOP_OVER_VOL_OWNED(gv, c, i) { + realnum pcur = p[i]; + // Assign current polarization by looping through first + // the numerator coefficients and then the denominator + // coefficients. + data->p[i] = numZ[0] * W[c2][cmp][i]; + for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; + for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; + + // shift the W field + // shift the P field + p[i] = gamma1inv * + (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + omega0dtsqr * (s[i] * w[i])); + pp[i] = pcur; + } + } + } + } + } +} + // constructor must take the s domain coefficients and transform // them to the z domain and store them in the class. -iir_susceptibility::iir_susceptibility(double *num, int num_N, double *den, int den_N, double dt){ - numz_N = num_N; - denz_N = den_N; - tustins_method(num,num_N,den,den_N, dt); - numz = num; - denz = den; +iir_susceptibility::iir_susceptibility(std::vector vec_numS, std::vector vec_denS, double dt) { + if (vec_numS.size() > vec_denS.size()) abort("error: the numerator rank cannot be higher than the denominator rank.\n"); + + // copy over the S domain numerator and denominator + N = vec_numS.size(); + D = vec_denS.size(); + N_z = D; + D_z = D; // both numerator and denominator of Z domain are same rank. + numS = new double[N]; + std::copy(vec_numS.begin(), vec_numS.end(), numS); + denS = new double[D]; + std::copy(vec_denS.begin(), vec_denS.end(), denS); + numZ = new double[N_z]; + denZ = new double[D_z]; + + T = dt; + + // calculate the Z domain numerator and denominator + tustins_method(numS, N, denS, D, numZ, N_z, denZ,D_z, T); } -*/ +// 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 45a08da93..a98a3de2d 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -364,18 +364,15 @@ class multilevel_susceptibility : public susceptibility { 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: - double *numz, *denz; - int numz_N, denz_N; + iir_susceptibility(std::vector numS, std::vector denS, double d); - iir_susceptibility(double *num, int num_N, double *den, int den_N, double dt); - - virtual susceptibility *clone() const { return new iir_susceptibility(*this); } - - virtual ~iir_susceptibility() {} + //virtual susceptibility *clone() const { return new iir_susceptibility(*this); } + virtual ~iir_susceptibility(); + virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *P_internal_data) const; @@ -396,16 +393,17 @@ class iir_susceptibility : public susceptibility { virtual void dump_params(h5file *h5f, size_t *start); virtual int get_num_params() { return 4; } -}; -*/ -// generalized method to convert from s domain to z domain. -void polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma); + -// tustins method (aka trapezoidal rule) to convert from s domain to z domain. -void tustins_method(double *num, int num_N, double *den, int den_N, double dt); +protected: +double *numS, *denS, *numZ, *denZ, T; +int N, D, N_z, D_z; -// backward difference method used to convert from s domain to z domain. -void backward_difference(double *num, int num_N, double *den, int den_N, double dt); +}; +// tustins method (aka trapezoidal rule) to convert from s domain to z domain. +void tustins_method(double* vec_numS, int N, double* vec_denS, int D, double* vec_numZ, int N_z, double* vec_denZ, int D_z, double T); +double factorial(int i); +double comb(int n, int k); class grace; diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index a85a1250d..cfdaf0743 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1385,18 +1385,23 @@ 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 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); } master_printf("%s%s susceptibility: frequency=%g, gamma=%g", noisy ? "noisy " : gyrotropic ? "gyrotropic " : "", ss->saturated_gyrotropy ? "Landau-Lifshitz-Gilbert-type" + : 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"); } From 924e95a50545579c9e8a573752d7bc72316a845f Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 25 Jun 2019 15:50:00 -0600 Subject: [PATCH 5/8] more progress --- src/iir-susceptibility.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index a3a17aac4..bac75af4d 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -204,12 +204,13 @@ void iir_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], } } else { // isotropic LOOP_OVER_VOL_OWNED(gv, c, i) { - realnum pcur = p[i]; + realnum pcur = data->p[i]; // Assign current polarization by looping through first - // the numerator coefficients and then the denominator - // coefficients. + // the numerator coefficients... data->p[i] = numZ[0] * W[c2][cmp][i]; for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; + //... and now throught the denominator coefficients. + data->p[i] = denZ[0] * for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; // shift the W field From f71a9a83dd4c9a35e46aaacadd461b8600dd69f8 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Thu, 25 Jul 2019 16:01:03 -0600 Subject: [PATCH 6/8] progress --- python/geom.py | 8 +- src/iir-susceptibility.cpp | 196 +++++++++++++++++++++++-------------- src/meep.hpp | 5 +- 3 files changed, 131 insertions(+), 78 deletions(-) diff --git a/python/geom.py b/python/geom.py index 9a690762a..d4d60698b 100755 --- a/python/geom.py +++ b/python/geom.py @@ -309,15 +309,17 @@ def __init__(self, num, den, **kwargs): self.den = np.squeeze(den).tolist() 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.float64) eval_den = np.zeros(freq.shape,dtype=np.float64) N = self.num.size D = self.den.size - 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] + for k in range(N): eval_num += pow(1j*freq,k)*self.num[k] + for k in range(D): eval_den += pow(1j*freq,k)*self.den[k] - return eval_num / eval_den + return sigma * eval_num / eval_den class DrudeSusceptibility(Susceptibility): diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index bac75af4d..049f6d9a6 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -42,7 +42,6 @@ double 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( @@ -103,21 +102,28 @@ void tustins_method( // IIR transfer function susceptibility // --------------------------------------------------------------- // +// The internal data needs to contain the current polarizations, +// the previous polarizations, and the previous fields. typedef struct { size_t sz_data; - size_t ntot; - realnum *P; - realnum *W; + realnum *P[NUM_FIELD_COMPONENTS][2]; // current polarization value used in step equations + size_t P_tot; + realnum *P_prev[NUM_FIELD_COMPONENTS][2]; // previous polarization values used to update + size_t P_prev_tot; + realnum *W_prev[NUM_FIELD_COMPONENTS][2]; // previous field values used to update + size_t W_prev_tot; realnum data[1]; } iir_data; -// The internal data is just a backup of P AND W from the previous timestep(s). void *iir_susceptibility::new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const { + // Count to see how large we need our data structure to be int num = 0; FOR_COMPONENTS(c) DOCMP2 { - if (needs_P(c, cmp, W)) num += 2 * gv.ntot() * 2 * Dz; // for each component, we need 2 directions, ntot index points, and 2*Dz filter taps (1 for P, 1 for W) + // for each component, we need 2 directions, ntot index points, + // Dz denominator filter taps, and Nz-1 numerator filter taps + if (needs_P(c, cmp, W)) num += 2 * gv.ntot() * D_z * (N_z-1); } - size_t sz = sizeof(iir_susceptibility) + sizeof(realnum) * (num - 1); + size_t sz = sizeof(iir_data) + sizeof(realnum) * (num); iir_data *d = (iir_data *)malloc(sz); d->sz_data = sz; return (void *)d; @@ -126,98 +132,91 @@ void *iir_susceptibility::new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], void iir_susceptibility::init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *data) const { (void)dt; // unused 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() * Dz; + + // Take that contiguous block of data passed by the user and + // distribute it to the appropriate array locations + d->P_tot = gv.ntot(); + d->P_prev_tot = gv.ntot() * D_z; + d->W_prev_tot = gv.ntot() * (N_z -1); realnum *P = d->data; - realnum *P_prev = d->data + ntot; - realnum *W_prev = P_prev + ntot; + realnum *P_prev = d->data + d->P_tot; + realnum *W_prev = P_prev + d->P_prev_tot; + FOR_COMPONENTS(c) DOCMP2 { if (needs_P(c, cmp, W)) { for (int dd = X; dd < R; dd++) { d->P[c][cmp] = P; d->P_prev[c][cmp] = P_prev; d->W_prev[c][cmp] = W_prev; - P += 2 * ntot; - P_prev += 2 * ntot; + P += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + P_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + W_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); } } } } +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->P_tot = d->P_tot; + dnew->P_prev_tot = d->P_prev_tot; + dnew->W_prev_tot = d->W_prev_tot; + realnum *P = d->data; + realnum *P_prev = d->data + d->P_tot; + realnum *W_prev = P_prev + d->P_prev_tot; + FOR_COMPONENTS(c) DOCMP2 { + if (d->P[c][cmp]) { + dnew->P[c][cmp] = P; + dnew->P_prev[c][cmp] = P_prev; + dnew->W_prev[c][cmp] = W_prev; + P += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + P_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + W_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + } + } + return (void *)dnew; +} + 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; (void)W_prev; // unused; - + (void)dt; 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], *pp = d->P_prev[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; + if (s != 0) { + realnum *p = d->P[c][cmp], *pp = d->P_prev[c][cmp], *wp = d->W_prev[c][cmp]; - 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 - LOOP_OVER_VOL_OWNED(gv, c, i) { - // s[i] != 0 check is a bit of a hack to work around - // some instabilities that occur near the boundaries - // of materials; see PR #666 - if (s[i] != 0) { - realnum pcur = p[i]; - p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + - omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is) + - OFFDIAG(s2, w2, is2, is))); - pp[i] = pcur; - } - } - } else if (s1) { // 2x2 anisotropic LOOP_OVER_VOL_OWNED(gv, c, i) { - if (s[i] != 0) { // see above - realnum pcur = p[i]; - p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + - omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is))); - pp[i] = pcur; - } - } - } else { // isotropic - LOOP_OVER_VOL_OWNED(gv, c, i) { - realnum pcur = data->p[i]; - // Assign current polarization by looping through first - // the numerator coefficients... - data->p[i] = numZ[0] * W[c2][cmp][i]; - for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; - //... and now throught the denominator coefficients. - data->p[i] = denZ[0] * - for (int num_c=0; num_cp[i] += numZ[num_c+1] * data->W_prev[num_c]; - - // shift the W field - // shift the P field - p[i] = gamma1inv * - (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + omega0dtsqr * (s[i] * w[i])); - pp[i] = pcur; + // Initialize current polarization by weighting the + // current field value with the first numerator tap ... + p[i] = numZ[0] * W[c][cmp][i]; + // ... then continue by looping through the rest of + // the numerator taps... + for (int num_c=0; num_c0; num_c--) wp[num_c][i] = wp[num_c-1][i]; + wp[0][i] = W[c][cmp][i]; + // shift the P field + for (int den_c=D_z-1; den_c>0; den_c--) pp[den_c][i] = pp[den_c-1][i]; + pp[0][i] = p[i]; } } } @@ -225,10 +224,29 @@ void iir_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], } } +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) { - if (vec_numS.size() > vec_denS.size()) abort("error: the numerator rank cannot be higher than the denominator rank.\n"); + 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 N = vec_numS.size(); @@ -248,6 +266,38 @@ iir_susceptibility::iir_susceptibility(std::vector vec_numS, std::vector tustins_method(numS, N, denS, D, numZ, N_z, denZ,D_z, T); } +int iir_susceptibility::num_cinternal_notowned_needed(component c, + void *P_internal_data) const { + iir_data *d = (iir_data *)P_internal_data; + return d->P[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),i) * numS[i]; + for (int i=0;i(0,freq),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; diff --git a/src/meep.hpp b/src/meep.hpp index a98a3de2d..d930fa868 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -369,11 +369,12 @@ class iir_susceptibility : public susceptibility { public: iir_susceptibility(std::vector numS, std::vector denS, double d); - //virtual susceptibility *clone() const { return new iir_susceptibility(*this); } + virtual susceptibility *clone() const { return new iir_susceptibility(*this); } virtual ~iir_susceptibility(); - virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, + 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], From d0eb4f7a2647770ed93d4a9bab44d7bdd5f15a55 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Wed, 31 Jul 2019 13:43:04 -0600 Subject: [PATCH 7/8] checkpoint --- python/geom.py | 20 ++-- src/iir-susceptibility.cpp | 234 ++++++++++++++++++++++++------------- src/meep.hpp | 17 ++- src/susceptibility.cpp | 18 +++ 4 files changed, 191 insertions(+), 98 deletions(-) diff --git a/python/geom.py b/python/geom.py index 00b77f310..5070f04cd 100755 --- a/python/geom.py +++ b/python/geom.py @@ -309,21 +309,21 @@ def eval_susceptibility(self,freq): class IIR_Susceptibility(Susceptibility): def __init__(self, num, den, **kwargs): super(IIR_Susceptibility, self).__init__(**kwargs) - self.num = np.squeeze(num).tolist() - self.den = np.squeeze(den).tolist() + 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.float64) - eval_den = np.zeros(freq.shape,dtype=np.float64) - N = self.num.size - D = self.den.size - for k in range(N): eval_num += pow(1j*freq,k)*self.num[k] - for k in range(D): eval_den += pow(1j*freq,k)*self.den[k] - - return sigma * eval_num / eval_den + 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/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index 049f6d9a6..2c9deb2ff 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -32,39 +32,89 @@ namespace meep { * in the s domain to the z domain. */ -double factorial(int i) { - double factorial = 1; - for (i; i > 0; i--) factorial *= i; + // 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 polynomial_transform(double *r, int N, double alpha, double beta, double delta, double gamma) { + 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; } -double comb(int n, int k) { +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( - double* vec_numS, + realnum* vec_numS, int N, - double* vec_denS, + realnum* vec_denS, int D, - double* vec_numZ, + realnum* vec_numZ, int N_z, - double* vec_denZ, + realnum* vec_denZ, int D_z, - double T + 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); + std::vector bprime(Np + 1, 0.0); + std::vector aprime(Dp + 1, 0.0); int j, i, k, l; - double val; + realnum val; for (j=0;j<=Np;j++) { val = 0.0; for(i=0;isz_data = sz; + d->ntot = ntot; + 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; + size_t ntot = d->ntot; memset(d, 0, sz_data); d->sz_data = sz_data; + d->ntot = ntot; // Take that contiguous block of data passed by the user and // distribute it to the appropriate array locations - d->P_tot = gv.ntot(); - d->P_prev_tot = gv.ntot() * D_z; - d->W_prev_tot = gv.ntot() * (N_z -1); realnum *P = d->data; - realnum *P_prev = d->data + d->P_tot; - realnum *W_prev = P_prev + d->P_prev_tot; FOR_COMPONENTS(c) DOCMP2 { if (needs_P(c, cmp, W)) { - for (int dd = X; dd < R; dd++) { d->P[c][cmp] = P; - d->P_prev[c][cmp] = P_prev; - d->W_prev[c][cmp] = W_prev; - P += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); - P_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); - W_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); - } + P += ntot; + d->P_prev[c][cmp] = P; + P += ntot * (D_z-1); + d->W_prev[c][cmp] = P; + P += ntot * (N_z-1); } } } @@ -167,20 +213,18 @@ void *iir_susceptibility::copy_internal_data(void *data) const { iir_data *dnew = (iir_data *)malloc(d->sz_data); memcpy(dnew, d, d->sz_data); - dnew->P_tot = d->P_tot; - dnew->P_prev_tot = d->P_prev_tot; - dnew->W_prev_tot = d->W_prev_tot; - realnum *P = d->data; - realnum *P_prev = d->data + d->P_tot; - realnum *W_prev = P_prev + d->P_prev_tot; + 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; - dnew->P_prev[c][cmp] = P_prev; - dnew->W_prev[c][cmp] = W_prev; - P += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); - P_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); - W_prev += 2 * (d->P_tot + d->P_prev_tot + d->W_prev_tot); + P += ntot; + dnew->P_prev[c][cmp] = P; + P += ntot * (D_z-1); + dnew->W_prev[c][cmp] = P; + P += ntot * (N_z-1); } } return (void *)dnew; @@ -190,34 +234,42 @@ 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; - (void)W_prev; // unused; + 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) { - if (s != 0) { - realnum *p = d->P[c][cmp], *pp = d->P_prev[c][cmp], *wp = d->W_prev[c][cmp]; - - LOOP_OVER_VOL_OWNED(gv, c, i) { - // Initialize current polarization by weighting the - // current field value with the first numerator tap ... - p[i] = numZ[0] * W[c][cmp][i]; - // ... then continue by looping through the rest of - // the numerator taps... - for (int num_c=0; num_c0; num_c--) wp[num_c][i] = wp[num_c-1][i]; - wp[0][i] = W[c][cmp][i]; - // shift the P field - for (int den_c=D_z-1; den_c>0; den_c--) pp[den_c][i] = pp[den_c-1][i]; - pp[0][i] = p[i]; + realnum *p = d->P[c][cmp], *pp = d->P_prev[c][cmp], *wp = d->W_prev[c][cmp]; + LOOP_OVER_VOL_OWNED(gv, c, i) { + //p[i] = 0; + if (s[i] != 0) { + for (int den_c=0; den_c0; num_c--) wp[num_c * ntot + i] = wp[(num_c-1) * ntot + i]; + wp[0 * ntot+i] = w[i]; + // shift the P field + for (int den_c=D_z-2; den_c>0; den_c--) pp[den_c * ntot + i] = pp[(den_c-1) * ntot + i]; + pp[0 * ntot+i] = p[i]; + //master_printf("W[%i]: %3.1e\n",i,w[i]); } } } @@ -245,25 +297,42 @@ void iir_susceptibility::subtract_P(field_type ft, // 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) { - if (vec_numS.size() > vec_denS.size()) abort("error: for iir-filter susceptibilities, the numerator rank cannot be higher than the denominator rank.\n"); +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 - N = vec_numS.size(); - D = vec_denS.size(); - N_z = D; - D_z = D; // both numerator and denominator of Z domain are same rank. - numS = new double[N]; - std::copy(vec_numS.begin(), vec_numS.end(), numS); - denS = new double[D]; - std::copy(vec_denS.begin(), vec_denS.end(), denS); - numZ = new double[N_z]; - denZ = new double[D_z]; - - T = dt; + // copy over the S domain numerator and denominator vectors to arrays + for (int i=0; i iir_susceptibility::chi1(double freq, double sigma) { std::complex num = 0; std::complex den = 0; - for (int i=0;i(0,freq),i) * numS[i]; - for (int i=0;i(0,freq),i) * denS[i]; + 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; } diff --git a/src/meep.hpp b/src/meep.hpp index fbcb6b6d7..dd4f07ca0 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -378,7 +378,8 @@ class multilevel_susceptibility : public susceptibility { class iir_susceptibility : public susceptibility { public: - iir_susceptibility(std::vector numS, std::vector denS, double d); + 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); } @@ -405,17 +406,21 @@ class iir_susceptibility : public susceptibility { 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: -double *numS, *denS, *numZ, *denZ, T; 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 tustins_method(double* vec_numS, int N, double* vec_denS, int D, double* vec_numZ, int N_z, double* vec_denZ, int D_z, double T); -double factorial(int i); -double comb(int n, int k); +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; diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index 80a88543f..78b815e7a 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; @@ -253,6 +261,8 @@ void lorentzian_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] + omega0dtsqr * (s[i] * w[i])); pp[i] = pcur; + //master_printf("w_0: %3.3e, p_1: %3.3e, p_2: %3.3e\n",gamma1inv * omega0dtsqr * s[i] ,(2 - omega0dtsqr_denom) * gamma1inv, -gamma1 * gamma1inv); + //master_printf("pcure: %3.6e\n",pcur); } } } @@ -313,6 +323,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 +360,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) { From 02c693da71a182e967939740a2f5dd513ba55274 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 2 Aug 2019 16:06:24 -0600 Subject: [PATCH 8/8] stuck --- src/iir-susceptibility.cpp | 118 ++++++++++++++++++++++++------------- src/susceptibility.cpp | 2 - 2 files changed, 78 insertions(+), 42 deletions(-) diff --git a/src/iir-susceptibility.cpp b/src/iir-susceptibility.cpp index 2c9deb2ff..3638c63dd 100644 --- a/src/iir-susceptibility.cpp +++ b/src/iir-susceptibility.cpp @@ -36,7 +36,18 @@ namespace meep { // 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"); @@ -54,7 +65,10 @@ void polynomial_transform(double *r, int N, double alpha, double beta, double de for (j=0; jsz_data = sz; - d->ntot = ntot; return (void *)d; } @@ -186,10 +203,9 @@ void iir_susceptibility::init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], // initialize everything to 0 size_t sz_data = d->sz_data; - size_t ntot = d->ntot; memset(d, 0, sz_data); d->sz_data = sz_data; - d->ntot = ntot; + 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 @@ -199,10 +215,8 @@ void iir_susceptibility::init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], if (needs_P(c, cmp, W)) { d->P[c][cmp] = P; P += ntot; - d->P_prev[c][cmp] = P; + d->storage[c][cmp] = P; P += ntot * (D_z-1); - d->W_prev[c][cmp] = P; - P += ntot * (N_z-1); } } } @@ -221,15 +235,24 @@ void *iir_susceptibility::copy_internal_data(void *data) const { if (d->P[c][cmp]) { dnew->P[c][cmp] = P; P += ntot; - dnew->P_prev[c][cmp] = P; + dnew->storage[c][cmp] = P; P += ntot * (D_z-1); - dnew->W_prev[c][cmp] = P; - P += ntot * (N_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 { @@ -242,34 +265,48 @@ void iir_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2], 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], *pp = d->P_prev[c][cmp], *wp = d->W_prev[c][cmp]; + 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) { - //p[i] = 0; if (s[i] != 0) { - for (int den_c=0; den_c0; num_c--) wp[num_c * ntot + i] = wp[(num_c-1) * ntot + i]; - wp[0 * ntot+i] = w[i]; - // shift the P field - for (int den_c=D_z-2; den_c>0; den_c--) pp[den_c * ntot + i] = pp[(den_c-1) * ntot + i]; - pp[0 * ntot+i] = p[i]; - //master_printf("W[%i]: %3.1e\n",i,w[i]); } } } @@ -311,11 +348,12 @@ iir_susceptibility::iir_susceptibility(std::vector vec_numS, std::vector // calculate the Z domain numerator and denominator for (int i=0; i