From 56812253790e50a519105cae37763047ecd10856 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 11:22:43 +0200 Subject: [PATCH 01/14] Enabled --prefix= option passed by autogen.sh - important for experimental compiling --- autogen.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/autogen.sh b/autogen.sh index 143d469e1..b7e0e4ebd 100755 --- a/autogen.sh +++ b/autogen.sh @@ -9,6 +9,7 @@ while test $# -ge 1; do --disable-*) configure_args="$configure_args $1" ;; --with-*) configure_args="$configure_args $1" ;; --without-*) configure_args="$configure_args $1" ;; + --prefix=*) configure_args="$configure_args $1" ;; *) echo "unknown argument $1"; exit 1 ;; esac shift From 281658610b056502de106b91d2f6358eea75c5cf Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 11:26:59 +0200 Subject: [PATCH 02/14] Fixed bicgstab according to Wu http://permalink.gmane.org/gmane.comp.science.electromagnetism.meep.general/4787 --- src/bicgstab.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/bicgstab.cpp b/src/bicgstab.cpp index 3004bfcd7..d7aa67579 100644 --- a/src/bicgstab.cpp +++ b/src/bicgstab.cpp @@ -79,7 +79,8 @@ static double dot(int n, const realnum *x, const realnum *y) } static double norm2(int n, const realnum *x) { - int i; + return sqrt(dot(n, x, x)); + /* int i; // did not work properly double xmax = 0, scale; long double sum = 0; for (i = 0; i < n; ++i) { @@ -93,7 +94,7 @@ static double norm2(int n, const realnum *x) { double xs = scale * x[i]; sum += xs * xs; } - return xmax * sqrt(sum_to_all(sum)); + return xmax * sqrt(sum_to_all(sum)); */ } static void xpay(int n, realnum *x, double a, const realnum *y) { From 746ca7011b3f0aa18452aa85900529166cc0d75e Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 13:37:35 +0200 Subject: [PATCH 03/14] Added support for a source with rectangular spectrum (see http://fzu.cz/~dominecf/misc/meep/index.html#bandsource) --- src/meep.hpp | 18 +++++++++++++ src/sources.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 70204b509..982ef3f98 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -737,6 +737,24 @@ class gaussian_src_time : public src_time { double freq, width, peak_time, cutoff; }; + +// band source with given minimum and maximum frequency and duration +class band_src_time : public src_time { + public: + band_src_time(double f, double fwidth, double user_cutoff); + virtual ~band_src_time() {} + + virtual complex dipole(double time) const; + virtual double last_time() const { return float(peak_time + cutoff); }; + virtual src_time *clone() const { return new band_src_time(*this); } + virtual bool is_equal(const src_time &t) const; + virtual complex frequency() const { return freq; } + virtual void set_frequency(complex f) { freq = real(f); } + + private: + double freq, width, peak_time, cutoff; + //complex test; +}; // Continuous (CW) source with (optional) slow turn-on and/or turn-off. class continuous_src_time : public src_time { public: diff --git a/src/sources.cpp b/src/sources.cpp index c3c7b619d..aabe9c9da 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -23,6 +23,12 @@ #include "meep.hpp" #include "meep_internals.hpp" +// For sine-integral function used in the band-source +#define HAVE_LIBGSL_EXPINT + +#ifdef HAVE_LIBGSL_EXPINT +# include +#endif namespace meep { /*********************************************************************/ @@ -107,6 +113,70 @@ bool gaussian_src_time::is_equal(const src_time &t) const else return 0; } +/// Band source +band_src_time::band_src_time(double f, double fwidth, double user_cutoff) +{ + freq = f; + width = 1.0 / fwidth; + peak_time = user_cutoff/2; + cutoff = user_cutoff; +#ifdef HAVE_LIBGSL_EXPINT + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tExperimental: using GSL si() function for flat-top spectrum\n"); +#else + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tWarning: not compiled with GSL, the source spectrum will not be flat-top\n", peak_time, cutoff); +#endif + cutoff = float(cutoff); // don't make cutoff sensitive to roundoff error +} + +complex band_src_time::dipole(double time) const +{ + double tt = time - peak_time; + + // The function that introduces a rectangular band in the spectrum (centered around zero frequency) +#ifdef HAVE_LIBGSL_EXPINT + // The emitted field is the derivative of the dipole, so if possible, we use the sine integral: + // Si(x) = \int_0^x dt sin(t)/t + complex func = gsl_sf_Si(tt*2*pi*(freq-.5/width)) - gsl_sf_Si(tt*2*pi*(freq+.5/width)); +#else + // If GSL not available, a reasonable approximation is the sinc(t) function (but has not flat top) + complex func = sin(tt*2*pi / width/2)/(tt*2*pi * width) * polar(1.0, -2*pi*freq*tt); +#endif + + // The envelope that suppresses ringing and side lobes of the rectangle + double wnd_BlackmanN= (0.3635819 + 0.4891775*cos(tt/(cutoff)*pi*2) + + 0.1365995*cos(tt/(cutoff)*pi*4)+ 0.0106411*cos(tt/(cutoff)*pi*6)); + //double wnd_Hann = (.5 + .5*cos(tt/(cutoff)*pi*2)); + //double wnd_Hamming = (0.53836+0.46164*cos(tt/(cutoff)*pi*2)) ; + //double wnd_Blackman = (0.42659 + 0.49656*cos(tt/(cutoff)*pi*2) + 0.076849*cos(tt/(cutoff)*pi*4)) ; + //double wnd_Gauss = exp(-((tt)/(cutoff)*3)**2) // Gaussian window, char. width = 1/3; + //double wnd_ContRect = // should crop sinc when passing zero: rect(tt, 0, (trunc((cutoff)/8/width))*8*width) ; + + // correction factor so that current amplitude (= d(dipole)/dt) is + // ~ 1 near the peak of the band. + complex amp = 1.0 / complex(0,-2*pi*freq); // TODO + + // The needed time-domain source amplitude is the sinc function constrained by window and shifted + // in frequency by complex exponential + if (abs(tt) < cutoff/2) + //return sinc; + return func * wnd_BlackmanN; + //return func * wnd_BlackmanN * amp; + else + return 0; + +} + +bool band_src_time::is_equal(const src_time &t) const +{ + const band_src_time *tp = dynamic_cast(&t); + if (tp) + return(tp->freq == freq && tp->width == width && + tp->peak_time == peak_time && tp->cutoff == cutoff); + else + return 0; +} complex continuous_src_time::dipole(double time) const { From 39ddb6850e9519d567be5b1d98ca2c8140177095 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 5 Feb 2015 20:41:01 +0100 Subject: [PATCH 04/14] Fixing autogen.sh to upstream --- autogen.sh | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/autogen.sh b/autogen.sh index b7e0e4ebd..bac4c3f0f 100755 --- a/autogen.sh +++ b/autogen.sh @@ -4,13 +4,13 @@ configure_args="" while test $# -ge 1; do case $1 in - --verbose) verbose=yes ;; - --enable-*) configure_args="$configure_args $1" ;; - --disable-*) configure_args="$configure_args $1" ;; - --with-*) configure_args="$configure_args $1" ;; - --without-*) configure_args="$configure_args $1" ;; - --prefix=*) configure_args="$configure_args $1" ;; - *) echo "unknown argument $1"; exit 1 ;; + --verbose) verbose=yes ;; + --enable-*) configure_args="$configure_args $1" ;; + --disable-*) configure_args="$configure_args $1" ;; + --with-*) configure_args="$configure_args $1" ;; + --without-*) configure_args="$configure_args $1" ;; + --prefix=*) configure_args="$configure_args $1" ;; + *) echo "unknown argument $1"; exit 1 ;; esac shift done From cb974764c4aedc02390f0125a1ec5507905674b6 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 11:26:59 +0200 Subject: [PATCH 05/14] Fixed bicgstab according to Wu http://permalink.gmane.org/gmane.comp.science.electromagnetism.meep.general/4787 --- src/bicgstab.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bicgstab.cpp b/src/bicgstab.cpp index c66047fd0..6fa651bd1 100644 --- a/src/bicgstab.cpp +++ b/src/bicgstab.cpp @@ -96,7 +96,7 @@ static double norm2(int n, const realnum *x) { double xs = scale * x[i]; sum += xs * xs; } - return xmax * sqrt(sum_to_all(sum)); + return xmax * sqrt(sum_to_all(sum)); } static void xpay(int n, realnum *x, double a, const realnum *y) { From f9685b97e50172f93064491b12eb124a3ea0c52a Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 13:37:35 +0200 Subject: [PATCH 06/14] Added support for a source with rectangular spectrum (see http://fzu.cz/~dominecf/misc/meep/index.html#bandsource) --- src/meep.hpp | 19 +++++++++++++++++++ src/sources.cpp | 1 + 2 files changed, 20 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 62ee6421b..ae7d89f1c 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -755,6 +755,25 @@ class band_src_time : public src_time { double freq, width, peak_time, cutoff; //complex test; }; + +// band source with given minimum and maximum frequency and duration +class band_src_time : public src_time { + public: + band_src_time(double f, double fwidth, double user_cutoff); + virtual ~band_src_time() {} + + virtual complex dipole(double time) const; + virtual double last_time() const { return float(peak_time + cutoff); }; + virtual src_time *clone() const { return new band_src_time(*this); } + virtual bool is_equal(const src_time &t) const; + virtual complex frequency() const { return freq; } + virtual void set_frequency(complex f) { freq = real(f); } + + private: + double freq, width, peak_time, cutoff; + //complex test; +}; + // Continuous (CW) source with (optional) slow turn-on and/or turn-off. class continuous_src_time : public src_time { public: diff --git a/src/sources.cpp b/src/sources.cpp index 20c6f7142..ec7acace7 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -29,6 +29,7 @@ #ifdef HAVE_LIBGSL_EXPINT # include #endif + using namespace std; namespace meep { From 72c308516bbd463c9e21e655e33e617ef6dfb979 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Thu, 10 Apr 2014 13:37:35 +0200 Subject: [PATCH 07/14] Added support for a source with rectangular spectrum (see http://fzu.cz/~dominecf/misc/meep/index.html#bandsource) --- src/meep.hpp | 18 +++++++++++++ src/sources.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 70204b509..982ef3f98 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -737,6 +737,24 @@ class gaussian_src_time : public src_time { double freq, width, peak_time, cutoff; }; + +// band source with given minimum and maximum frequency and duration +class band_src_time : public src_time { + public: + band_src_time(double f, double fwidth, double user_cutoff); + virtual ~band_src_time() {} + + virtual complex dipole(double time) const; + virtual double last_time() const { return float(peak_time + cutoff); }; + virtual src_time *clone() const { return new band_src_time(*this); } + virtual bool is_equal(const src_time &t) const; + virtual complex frequency() const { return freq; } + virtual void set_frequency(complex f) { freq = real(f); } + + private: + double freq, width, peak_time, cutoff; + //complex test; +}; // Continuous (CW) source with (optional) slow turn-on and/or turn-off. class continuous_src_time : public src_time { public: diff --git a/src/sources.cpp b/src/sources.cpp index c3c7b619d..aabe9c9da 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -23,6 +23,12 @@ #include "meep.hpp" #include "meep_internals.hpp" +// For sine-integral function used in the band-source +#define HAVE_LIBGSL_EXPINT + +#ifdef HAVE_LIBGSL_EXPINT +# include +#endif namespace meep { /*********************************************************************/ @@ -107,6 +113,70 @@ bool gaussian_src_time::is_equal(const src_time &t) const else return 0; } +/// Band source +band_src_time::band_src_time(double f, double fwidth, double user_cutoff) +{ + freq = f; + width = 1.0 / fwidth; + peak_time = user_cutoff/2; + cutoff = user_cutoff; +#ifdef HAVE_LIBGSL_EXPINT + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tExperimental: using GSL si() function for flat-top spectrum\n"); +#else + master_printf("Initializing band source for time_domain (peak_time = %g, cutoff = %g)\n", peak_time, cutoff); + master_printf("\tWarning: not compiled with GSL, the source spectrum will not be flat-top\n", peak_time, cutoff); +#endif + cutoff = float(cutoff); // don't make cutoff sensitive to roundoff error +} + +complex band_src_time::dipole(double time) const +{ + double tt = time - peak_time; + + // The function that introduces a rectangular band in the spectrum (centered around zero frequency) +#ifdef HAVE_LIBGSL_EXPINT + // The emitted field is the derivative of the dipole, so if possible, we use the sine integral: + // Si(x) = \int_0^x dt sin(t)/t + complex func = gsl_sf_Si(tt*2*pi*(freq-.5/width)) - gsl_sf_Si(tt*2*pi*(freq+.5/width)); +#else + // If GSL not available, a reasonable approximation is the sinc(t) function (but has not flat top) + complex func = sin(tt*2*pi / width/2)/(tt*2*pi * width) * polar(1.0, -2*pi*freq*tt); +#endif + + // The envelope that suppresses ringing and side lobes of the rectangle + double wnd_BlackmanN= (0.3635819 + 0.4891775*cos(tt/(cutoff)*pi*2) + + 0.1365995*cos(tt/(cutoff)*pi*4)+ 0.0106411*cos(tt/(cutoff)*pi*6)); + //double wnd_Hann = (.5 + .5*cos(tt/(cutoff)*pi*2)); + //double wnd_Hamming = (0.53836+0.46164*cos(tt/(cutoff)*pi*2)) ; + //double wnd_Blackman = (0.42659 + 0.49656*cos(tt/(cutoff)*pi*2) + 0.076849*cos(tt/(cutoff)*pi*4)) ; + //double wnd_Gauss = exp(-((tt)/(cutoff)*3)**2) // Gaussian window, char. width = 1/3; + //double wnd_ContRect = // should crop sinc when passing zero: rect(tt, 0, (trunc((cutoff)/8/width))*8*width) ; + + // correction factor so that current amplitude (= d(dipole)/dt) is + // ~ 1 near the peak of the band. + complex amp = 1.0 / complex(0,-2*pi*freq); // TODO + + // The needed time-domain source amplitude is the sinc function constrained by window and shifted + // in frequency by complex exponential + if (abs(tt) < cutoff/2) + //return sinc; + return func * wnd_BlackmanN; + //return func * wnd_BlackmanN * amp; + else + return 0; + +} + +bool band_src_time::is_equal(const src_time &t) const +{ + const band_src_time *tp = dynamic_cast(&t); + if (tp) + return(tp->freq == freq && tp->width == width && + tp->peak_time == peak_time && tp->cutoff == cutoff); + else + return 0; +} complex continuous_src_time::dipole(double time) const { From 115f45c8103b8e98b317df349103a5c5027b4334 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Fri, 6 Feb 2015 16:46:34 +0100 Subject: [PATCH 08/14] update namespace std:: in meep header --- src/meep.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 62ee6421b..ad31eb594 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -744,12 +744,12 @@ class band_src_time : public src_time { band_src_time(double f, double fwidth, double user_cutoff); virtual ~band_src_time() {} - virtual complex dipole(double time) const; + virtual std::complex dipole(double time) const; virtual double last_time() const { return float(peak_time + cutoff); }; virtual src_time *clone() const { return new band_src_time(*this); } virtual bool is_equal(const src_time &t) const; - virtual complex frequency() const { return freq; } - virtual void set_frequency(complex f) { freq = real(f); } + virtual std::complex frequency() const { return freq; } + virtual void set_frequency(std::complex f) { freq = real(f); } private: double freq, width, peak_time, cutoff; From c5043fd6fe5965dc98a69e2d01691aeee93f0a9c Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Fri, 6 Feb 2015 16:46:34 +0100 Subject: [PATCH 09/14] update namespace std:: in meep header (+whitespace fix) --- src/bicgstab.cpp | 2 +- src/meep.hpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/bicgstab.cpp b/src/bicgstab.cpp index 6fa651bd1..c66047fd0 100644 --- a/src/bicgstab.cpp +++ b/src/bicgstab.cpp @@ -96,7 +96,7 @@ static double norm2(int n, const realnum *x) { double xs = scale * x[i]; sum += xs * xs; } - return xmax * sqrt(sum_to_all(sum)); + return xmax * sqrt(sum_to_all(sum)); } static void xpay(int n, realnum *x, double a, const realnum *y) { diff --git a/src/meep.hpp b/src/meep.hpp index 62ee6421b..ad31eb594 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -744,12 +744,12 @@ class band_src_time : public src_time { band_src_time(double f, double fwidth, double user_cutoff); virtual ~band_src_time() {} - virtual complex dipole(double time) const; + virtual std::complex dipole(double time) const; virtual double last_time() const { return float(peak_time + cutoff); }; virtual src_time *clone() const { return new band_src_time(*this); } virtual bool is_equal(const src_time &t) const; - virtual complex frequency() const { return freq; } - virtual void set_frequency(complex f) { freq = real(f); } + virtual std::complex frequency() const { return freq; } + virtual void set_frequency(std::complex f) { freq = real(f); } private: double freq, width, peak_time, cutoff; From 3f58eee1a768d4fe2963a7774a0c7170c9ec893d Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Fri, 6 Feb 2015 19:03:36 +0100 Subject: [PATCH 10/14] Lorentzian stability check should write a warning only (suspected to report false positives/negatives) --- src/susceptibility.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index 74ef0afbb..e2ad7f687 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -179,7 +179,7 @@ void lorentzian_susceptibility::update_P if (!no_omega_0_denominator && gamma >= 0 && lorentzian_unstable(omega_0, gamma, dt)) - abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); + master_printf("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); FOR_COMPONENTS(c) DOCMP2 if (d->P[c][cmp]) { const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)]; From 33429cea4dcc40634ffc4d61fc0d847e5f9ef36b Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Mon, 23 Feb 2015 20:15:23 +0100 Subject: [PATCH 11/14] removing the Lorentzian stability check (it is wrong) --- src/susceptibility.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index 74ef0afbb..c65dd67a5 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -177,9 +177,9 @@ void lorentzian_susceptibility::update_P const double omega0dtsqr_denom = no_omega_0_denominator ? 0 : omega0dtsqr; (void) W_prev; // unused; - if (!no_omega_0_denominator && gamma >= 0 - && lorentzian_unstable(omega_0, gamma, dt)) - abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); + //if (!no_omega_0_denominator && gamma >= 0 + //&& lorentzian_unstable(omega_0, gamma, dt)) + //abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); FOR_COMPONENTS(c) DOCMP2 if (d->P[c][cmp]) { const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)]; From e12f3fb7f90a4a070e92dee1424e92da63741004 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Mon, 23 Feb 2015 20:50:04 +0100 Subject: [PATCH 12/14] lorentzian_unstable() changed to a start-up warning only --- src/susceptibility.cpp | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index c65dd67a5..899ad87d9 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -111,6 +111,21 @@ void *lorentzian_susceptibility::new_internal_data( return (void*) d; } + +/* Return true if the discretized Lorentzian ODE is intrinsically unstable, + i.e. if it corresponds to a filter with a pole z outside the unit circle. + Note that the pole satisfies the quadratic equation: + (z + 1/z - 2)/dt^2 + g*(z - 1/z)/(2*dt) + w^2 = 0 + where w = 2*pi*omega_0 and g = 2*pi*gamma. It is just a little + algebra from this to get the condition for a root with |z| > 1. */ +static bool lorentzian_unstable(double omega_0, double gamma, double dt) { + double w = 2*pi*omega_0, g = 2*pi*gamma; + double g2 = g*dt/2, w2 = (w*dt)*(w*dt); + double b = (1 - w2/2) / (1 + g2), c = (1 - g2) / (1 + g2); + return b*b > c && 2*b*b - c + 2*fabs(b)*sqrt(b*b - c) > 1; +} + + void lorentzian_susceptibility::init_internal_data( realnum *W[NUM_FIELD_COMPONENTS][2], double dt, const grid_volume &gv, void *data) const { @@ -128,6 +143,11 @@ void lorentzian_susceptibility::init_internal_data( P += 2*ntot; P_prev += 2*ntot; } + master_printf("CHECKING LORENTZIAN STABILITY\n"); + if (!no_omega_0_denominator && gamma >= 0 + && lorentzian_unstable(omega_0, gamma, dt)) + master_printf("Lorentzian expected unstable by the Von Neumann stability analysis (omega_0=%g, gamma=%g, dt=%g)\n", omega_0, gamma, dt); + master_printf("END CHECKING LORENTZIAN STABILITY\n"); } void *lorentzian_susceptibility::copy_internal_data(void *data) const { @@ -147,18 +167,6 @@ void *lorentzian_susceptibility::copy_internal_data(void *data) const { return (void*) dnew; } -/* Return true if the discretized Lorentzian ODE is intrinsically unstable, - i.e. if it corresponds to a filter with a pole z outside the unit circle. - Note that the pole satisfies the quadratic equation: - (z + 1/z - 2)/dt^2 + g*(z - 1/z)/(2*dt) + w^2 = 0 - where w = 2*pi*omega_0 and g = 2*pi*gamma. It is just a little - algebra from this to get the condition for a root with |z| > 1. */ -static bool lorentzian_unstable(double omega_0, double gamma, double dt) { - double w = 2*pi*omega_0, g = 2*pi*gamma; - double g2 = g*dt/2, w2 = (w*dt)*(w*dt); - double b = (1 - w2/2) / (1 + g2), c = (1 - g2) / (1 + g2); - return b*b > c && 2*b*b - c + 2*fabs(b)*sqrt(b*b - c) > 1; -} #define SWAP(t,a,b) { t SWAP_temp = a; a = b; b = SWAP_temp; } @@ -177,10 +185,6 @@ void lorentzian_susceptibility::update_P const double omega0dtsqr_denom = no_omega_0_denominator ? 0 : omega0dtsqr; (void) W_prev; // unused; - //if (!no_omega_0_denominator && gamma >= 0 - //&& lorentzian_unstable(omega_0, gamma, dt)) - //abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, 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) { From 6ae68bd967e66b8cd29c9b0651684d9c258c877b Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Mon, 23 Feb 2015 21:03:21 +0100 Subject: [PATCH 13/14] cleaned the instability printout --- src/susceptibility.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index 899ad87d9..6822ea95d 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -143,11 +143,11 @@ void lorentzian_susceptibility::init_internal_data( P += 2*ntot; P_prev += 2*ntot; } - master_printf("CHECKING LORENTZIAN STABILITY\n"); + + // TODO - ensure this is printed only once for each oscillator, not for every chunk if (!no_omega_0_denominator && gamma >= 0 && lorentzian_unstable(omega_0, gamma, dt)) - master_printf("Lorentzian expected unstable by the Von Neumann stability analysis (omega_0=%g, gamma=%g, dt=%g)\n", omega_0, gamma, dt); - master_printf("END CHECKING LORENTZIAN STABILITY\n"); + master_printf("Warning: Lorentzian may be unstable according to the Von Neumann stability analysis (omega_0=%g, gamma=%g, dt=%g). Proceed with caution.\n", omega_0, gamma, dt); } void *lorentzian_susceptibility::copy_internal_data(void *data) const { From d04b4719e8f69792356bdd969508a11a47214b11 Mon Sep 17 00:00:00 2001 From: Filip Dominec Date: Tue, 23 Feb 2016 15:36:53 +0100 Subject: [PATCH 14/14] lorentzian pole printout --- src/susceptibility.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp index e2ad7f687..102fcb9ae 100644 --- a/src/susceptibility.cpp +++ b/src/susceptibility.cpp @@ -177,9 +177,9 @@ void lorentzian_susceptibility::update_P const double omega0dtsqr_denom = no_omega_0_denominator ? 0 : omega0dtsqr; (void) W_prev; // unused; - if (!no_omega_0_denominator && gamma >= 0 - && lorentzian_unstable(omega_0, gamma, dt)) - master_printf("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); + //if (!no_omega_0_denominator && gamma >= 0 + //&& lorentzian_unstable(omega_0, gamma, dt)) + //master_printf("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt); FOR_COMPONENTS(c) DOCMP2 if (d->P[c][cmp]) { const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)];