diff options
| -rw-r--r-- | Makefile | 11 | ||||
| -rw-r--r-- | fourier.cpp | 77 | ||||
| -rw-r--r-- | fourier.hpp | 32 | ||||
| -rw-r--r-- | log-fourier_integrator.cpp | 114 |
4 files changed, 228 insertions, 6 deletions
@@ -1,4 +1,4 @@ -all: walk correlations integrator fourier_integrator get_energy +all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 @@ -9,13 +9,16 @@ correlations: correlations.cpp $(CC) correlations.cpp -o correlations integrator: integrator.cpp fourier.o fourier.hpp - $(CC) integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o integrator + $(CC) integrator.cpp fourier.o -lfftw3 -lgsl -o integrator fourier.o: fourier.cpp fourier.hpp $(CC) fourier.cpp -c -o fourier.o fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o - $(CC) fourier_integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o fourier_integrator + $(CC) fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o fourier_integrator + +log-fourier_integrator: log-fourier_integrator.cpp fourier.hpp fourier.o + $(CC) log-fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o log-fourier_integrator get_energy: get_energy.cpp fourier.hpp fourier.o - $(CC) get_energy.cpp fourier.o -lfftw3 -lfftw3_omp -o get_energy + $(CC) get_energy.cpp fourier.o -lfftw3 -lgsl -o get_energy diff --git a/fourier.cpp b/fourier.cpp index bc4b633..cf45ad1 100644 --- a/fourier.cpp +++ b/fourier.cpp @@ -1,4 +1,5 @@ #include "fourier.hpp" +#include <boost/multiprecision/fwd.hpp> inline Real fP(unsigned p, Real q) { return 0.5 * pow(q, p); @@ -27,8 +28,8 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q) { FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : n(n), Δω(Δω), Δτ(Δτ) { a = fftw_alloc_real(2 * n); â = reinterpret_cast<Complex*>(fftw_alloc_complex(n + 1)); - fftw_init_threads(); - fftw_plan_with_nthreads(FFTW_THREADS); +// fftw_init_threads(); +// fftw_plan_with_nthreads(FFTW_THREADS); fftw_import_wisdom_from_filename("fftw.wisdom"); plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a, reinterpret_cast<fftw_complex*>(â), flags); plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags); @@ -99,6 +100,78 @@ void FourierTransform::writeToA(unsigned i, Real ai) { a[i] = ai; } +LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs) : N(N), k(k), Δτ(Δτ), Δω(Δω), Δs(Δs) { + τₛ = -0.5 * N; + ωₛ = -0.5 * N; + sₛ = -0.5 * N; + a = reinterpret_cast<Complex*>(fftw_alloc_complex(N)); + â = reinterpret_cast<Complex*>(fftw_alloc_complex(N)); + fftw_import_wisdom_from_filename("fftw.wisdom"); + a_to_â = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), 1, 0); + â_to_a = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), 1, 0); + fftw_export_wisdom_to_filename("fftw.wisdom"); +} + +LogarithmicFourierTransform::~LogarithmicFourierTransform() { + fftw_destroy_plan(a_to_â); + fftw_destroy_plan(â_to_a); + fftw_free(a); + fftw_free(â); + fftw_cleanup(); +} + +Real LogarithmicFourierTransform::τ(unsigned n) const { + return Δτ * (n + τₛ); +} + +Real LogarithmicFourierTransform::ω(unsigned n) const { + return Δω * (n + ωₛ); +} + +Real LogarithmicFourierTransform::s(unsigned n) const { + return Δs * (n + sₛ); +} + +Real LogarithmicFourierTransform::t(unsigned n) const { + return exp(τ(n)); +} + +Real LogarithmicFourierTransform::ν(unsigned n) const { + return exp(ω(n)); +} + +Complex gamma(Complex z) { + gsl_sf_result logΓ; + gsl_sf_result argΓ; + + gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ); + + return exp(logΓ.val + 1i * argΓ.val); +} + +std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric){ + std::vector<Complex> ĉ(N); + std::vector<Real> σs = {1}; + if (symmetric){ + σs.push_back(-1); + } + for (Real σ : σs) { + for (unsigned n = 0; n < N; n++) { + a[n] = c[n] * exp((1 - k) * τ(n)); + } + fftw_execute(a_to_â); + for (unsigned n = 0; n < N; n++) { + â[(n + N / 2) % N] *= pow(-1i * σ, 1i * (s(n)) - k) * gamma(k - 1i * (s(n))); + } + fftw_execute(â_to_a); + for (unsigned n = 0; n < N; n++) { + ĉ[n] += exp(-k * ω(n) * 4 * M_PI) * a[n] / pow(2 * M_PI, 1); + } + } + + return ĉ; +} + std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) { return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat"; } diff --git a/fourier.hpp b/fourier.hpp index 791953b..22a57e7 100644 --- a/fourier.hpp +++ b/fourier.hpp @@ -4,6 +4,7 @@ #include <complex> #include <vector> #include <tuple> +#include <gsl/gsl_sf_gamma.h> #ifndef FFTW_THREADS #define FFTW_THREADS 1 @@ -36,6 +37,37 @@ public: std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R); }; +class LogarithmicFourierTransform { +private: + Complex* a; + Complex* â; + fftw_plan a_to_â; + fftw_plan â_to_a; + unsigned N; + Real k; + Real Δτ; + Real Δω; + Real Δs; + Real τₛ; + Real ωₛ; + Real sₛ; +public: + LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs); + ~LogarithmicFourierTransform(); + Real τ(unsigned n) const; + Real ω(unsigned n) const; + Real t(unsigned n) const; + Real ν(unsigned n) const; + Real s(unsigned n) const; + std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric); + std::vector<Complex> fourier(); + std::vector<Real> inverse(const std::vector<Complex>& ĉ); + void writeToA(unsigned i, Real ai); + std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R); +}; + +Complex gamma(Complex z); + std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ); Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ); std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ); diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp new file mode 100644 index 0000000..247270b --- /dev/null +++ b/log-fourier_integrator.cpp @@ -0,0 +1,114 @@ +#include "fourier.hpp" +#include <getopt.h> +#include <iostream> +#include <fstream> + +int main(int argc, char* argv[]) { + unsigned p = 2; + unsigned s = 2; + Real λ = 0.5; + Real τ₀ = 0; + Real β₀ = 0; + Real yₘₐₓ = 0.5; + Real Δy = 0.05; + + unsigned log2n = 8; + Real τₘₐₓ = 20; + + unsigned maxIterations = 1000; + Real ε = 1e-14; + Real γ = 1; + + int opt; + + while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) { + switch (opt) { + case 'p': + p = atoi(optarg); + break; + case 's': + s = atoi(optarg); + break; + case '2': + log2n = atoi(optarg); + break; + case 'T': + τₘₐₓ = atof(optarg); + break; + case 't': + τ₀ = atof(optarg); + break; + case '0': + β₀ = atof(optarg); + break; + case 'y': + yₘₐₓ = atof(optarg); + break; + case 'd': + Δy = atof(optarg); + break; + case 'I': + maxIterations = (unsigned)atof(optarg); + break; + case 'g': + γ = atof(optarg); + break; + default: + exit(1); + } + } + + unsigned N = pow(2, log2n); + + Real Δτ = 1e-2; + Real Δω = 1e-2; + Real Δs = 1e-2; + Real k = 0.1; + + Real Γ₀ = 1.0; + Real μ = Γ₀; + if (τ₀ > 0) { + μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); + } + + std::vector<Real> C(N); + std::vector<Real> R(N); + + LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs); + + for (unsigned n = 0; n < N; n++) { + if (τ₀ > 0) { + C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2)); + } else { + C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ; + } + R[n] = exp(-μ * fft.t(n)); + } + + std::vector<Complex> Ct = fft.fourier(C, true); + std::vector<Complex> Rt = fft.fourier(R, false); + + /* + for (unsigned n = 0; n < N; n++) { + std::cout << fft.t(n) << " " << C[n] << std::endl; + } + */ + + for (unsigned n = 0; n < N; n++) { + std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl; + } + + /* + // start from the exact solution for τ₀ = 0 + for (unsigned i = 0; i < N + 1; i++) { + Real ω = i * Δω; + Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2)); + Rt[i] = 1.0 / (μ + 1i * ω); + Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2)); + } + C = fft.inverse(Ct); + R = fft.inverse(Rt); + */ + + return 0; +} |