diff options
| -rw-r--r-- | Makefile | 8 | ||||
| -rw-r--r-- | log-fourier.cpp | 28 | ||||
| -rw-r--r-- | log-fourier.hpp | 1 | ||||
| -rw-r--r-- | log-fourier_test.cpp | 50 | ||||
| -rw-r--r-- | log_integrator.cpp | 54 |
5 files changed, 101 insertions, 40 deletions
@@ -1,4 +1,4 @@ -all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy +all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy log_integrator log-fourier_test CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 @@ -32,9 +32,15 @@ fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o $(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator +log-fourier_test: log-fourier_test.cpp log-fourier.hpp log-fourier.o + $(CC) log-fourier_test.cpp log-fourier.o -lfftw3 -lgsl -o log-fourier_test + log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o $(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long +log_integrator: log_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o + $(CC) log_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log_integrator + get_energy: get_energy.cpp fourier.hpp fourier.o p-spin.o $(CC) get_energy.cpp fourier.o p-spin.o -lfftw3 -lgsl -o get_energy diff --git a/log-fourier.cpp b/log-fourier.cpp index 9d1f2cb..bc2dd87 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -64,15 +64,13 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real for (unsigned n = 0; n < pad*N; n++) { if (n < N) { a[n] = c[n] * exp((1 - k) * τ(n)); - } else if (n >= (pad - 1) * N) { - a[n] = c[pad*N-n-1] * exp((1 - k) * τ(pad*N-n-1)); } else { a[n] = 0; } } FFTW_EXECUTE(a_to_â); for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= std::exp(II * s(n)) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { @@ -87,30 +85,16 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex std::vector<Real> c(N); std::vector<Real> σs = {1, -1}; for (Real σ : σs) { - for (unsigned n = 0; n < pad * N; n++) { - if (n < N) { - if (σ < 0) { - a[n] = std::conj(ĉ[n]) * std::exp((1 - k) * ω(n)); - } else { - a[n] = ĉ[n] * std::exp((1 - k) * ω(n)); - } - } else if (n >= (pad - 1) * N) { - if (σ < 0) { - a[n] = ĉ[pad*N-n-1] * std::exp((1 - k) * ω(pad*N-n-1)); - } else { - a[n] = std::conj(ĉ[pad*N-n-1]) * std::exp((1 - k) * ω(pad*N-n-1)); - } - } else { - a[n] = 0; - } + for (unsigned n = 0; n < N; n++) { + a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n)); } FFTW_EXECUTE(a_to_â); - for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n)); + for (unsigned n = 0; n < N; n++) { + â[(N / 2 + n) % N] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { - c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI); + c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N); } } diff --git a/log-fourier.hpp b/log-fourier.hpp index f0b53ef..675b1b0 100644 --- a/log-fourier.hpp +++ b/log-fourier.hpp @@ -17,6 +17,7 @@ private: unsigned pad; Real k; Real Δτ; + Real Δs; Real τₛ; Real ωₛ; Real sₛ; diff --git a/log-fourier_test.cpp b/log-fourier_test.cpp new file mode 100644 index 0000000..678e4b1 --- /dev/null +++ b/log-fourier_test.cpp @@ -0,0 +1,50 @@ +#include "log-fourier.hpp" +#include <getopt.h> +#include <iostream> + +int main(int argc, char* argv[]) { + /* Log-Fourier parameters */ + unsigned log2n = 8; + Real Δτ = 0.1; + Real k = 0.1; + unsigned pad = 4; + + int opt; + + while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:e:0:lS:x:P:")) != -1) { + switch (opt) { + case '2': + log2n = atoi(optarg); + break; + case 'k': + k = atof(optarg); + break; + case 'D': + Δτ = atof(optarg); + break; + case 'P': + pad = atoi(optarg); + break; + default: + exit(1); + } + } + + unsigned N =log2n; + + LogarithmicFourierTransform fft(N, k, Δτ, pad); + + std::vector<Complex> a(N); + + for (unsigned i = 0; i < N; i++) { + a[i] = 1 / (1 + pow(fft.ν(i), 2)); + } + + std::vector<Real> â = fft.inverse(a); + + for (unsigned i = 0; i < N; i++) { + std::cout << fft.t(i) << " " << â[i] << std::endl; + } + + return 0; +} diff --git a/log_integrator.cpp b/log_integrator.cpp index a8d9778..426af40 100644 --- a/log_integrator.cpp +++ b/log_integrator.cpp @@ -92,6 +92,7 @@ int main(int argc, char* argv[]) { std::vector<Real> Cₜ₋₁(N); std::vector<Real> Rₜ₋₁(N); + std::vector<Real> Γ(N); std::vector<Complex> Ĉₜ₋₁(N); std::vector<Complex> Ȓₜ₋₁(N); @@ -108,6 +109,7 @@ int main(int argc, char* argv[]) { Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁; } Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n)); + Γ[n] = (Γ₀ / τ₀) * exp(-fft.t(n) / τ₀); Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n)); @@ -142,23 +144,47 @@ int main(int argc, char* argv[]) { for (unsigned i = 0; i < N; i++) { dC[i] += -μₜ * Cₜ[i]; - Real ΓR; - for (unsigned j = 0; j < N; j++) { - + dR[i] += -μₜ * Rₜ[i]; + for (unsigned j = 1; j < N - i; j++) { + Real Γᵢ₊ⱼ = (Γ[i + j] + Γ[i + j - 1]) / 2; + signed jR = round(log(fft.t(i+j) - fft.t(i)) / Δτ) + N / 2; + Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR]; + Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR]; + Real ΓR2 = Γᵢ₊ⱼ * Rⱼ; + Real dfCR = (dfC[i + j] * Rⱼ + dfC[i + j - 1] * Rⱼ) / 2; + Real RddfCC = (RddfC[i + j] * Cⱼ + RddfC[i + j - 1] * Cⱼ) / 2; + dC[i] += (fft.t(i+j) - fft.t(i+j - 1)) * (ΓR2 + pow(β, 2) * (dfCR + RddfCC)); } + for (unsigned j = 1; j < i; j++) { + signed jR = round(log(fft.t(i) - fft.t(i-j)) / Δτ) + N / 2; + Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR]; + Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR]; + dC[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Cⱼ + RddfC[i - 1 - j+1] * Cⱼ) / 2; + dR[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Rⱼ + RddfC[i - 1 - j+1] * Rⱼ) / 2; + } + + if (dC[i] > 0) dC[i] = 0; + if (dR[i] > 0) dR[i] = 0; } + std::vector<Real> Rₜ₊₁(N); + std::vector<Real> Cₜ₊₁(N); - std::vector<Complex> Ĉₜ₊₁(N); - std::vector<Complex> Ȓₜ₊₁(N); - for (unsigned n = 0; n < N; n++) { - Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n)); - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n)); + Rₜ₊₁[0] = 1; + Cₜ₊₁[N - 1] = 0; + + for (unsigned i = 1; i < N; i++) { + Cₜ₊₁[N - 1 - i] = Cₜ₊₁[N - i] - (fft.t(N - i) - fft.t(N - i - 1)) * (dC[N - i] + dC[N - 1 - i]) / 2; + Rₜ₊₁[i] = Rₜ₊₁[i - 1] + (fft.t(i) - fft.t(i - 1)) * (dR[i - 1] + dR[i]) / 2; } - std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); - std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); - μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, x); + μₜ /= pow(tanh(Cₜ₊₁[0]-1)+1, x); + + Real C₀ = Cₜ₊₁[0]; + + for (unsigned i = 0; i < N; i++) { + Rₜ₊₁[i] = Rₜ₊₁[i] < 0 ? 0 : Rₜ₊₁[i]; + } ΔCₜ = 0; for (unsigned i = 0; i < N; i++) { @@ -183,8 +209,6 @@ int main(int argc, char* argv[]) { for (unsigned i = 0; i < N; i++) { Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]); Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]); - Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]); - Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]); } std::cerr << "\x1b[2K" << "\r"; @@ -195,8 +219,6 @@ int main(int argc, char* argv[]) { γ₀ /= 2; Cₜ = Cₜ₋₁; Rₜ = Rₜ₋₁; - Ĉₜ = Ĉₜ₋₁; - Ȓₜ = Ȓₜ₋₁; μₜ = μₜ₋₁; } else { Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β); @@ -209,8 +231,6 @@ int main(int argc, char* argv[]) { β += Δβ; Cₜ₋₁ = Cₜ; Rₜ₋₁ = Rₜ; - Ĉₜ₋₁ = Ĉₜ; - Ȓₜ₋₁ = Ȓₜ; μₜ₋₁ = μₜ; } } |