summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile8
-rw-r--r--log-fourier.cpp28
-rw-r--r--log-fourier.hpp1
-rw-r--r--log-fourier_test.cpp50
-rw-r--r--log_integrator.cpp54
5 files changed, 40 insertions, 101 deletions
diff --git a/Makefile b/Makefile
index 0497df5..c5cac94 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy log_integrator log-fourier_test
+all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy
CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
@@ -32,15 +32,9 @@ 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 bc2dd87..9d1f2cb 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -64,13 +64,15 @@ 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::exp(II * s(n)) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
+ â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
}
FFTW_EXECUTE(â_to_a);
for (unsigned n = 0; n < N; n++) {
@@ -85,16 +87,30 @@ 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 < N; n++) {
- a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n));
+ 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;
+ }
}
FFTW_EXECUTE(a_to_â);
- for (unsigned n = 0; n < N; n++) {
- â[(N / 2 + n) % N] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+ 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));
}
FFTW_EXECUTE(â_to_a);
for (unsigned n = 0; n < N; n++) {
- c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N);
+ c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
}
}
diff --git a/log-fourier.hpp b/log-fourier.hpp
index 675b1b0..f0b53ef 100644
--- a/log-fourier.hpp
+++ b/log-fourier.hpp
@@ -17,7 +17,6 @@ private:
unsigned pad;
Real k;
Real Δτ;
- Real Δs;
Real τₛ;
Real ωₛ;
Real sₛ;
diff --git a/log-fourier_test.cpp b/log-fourier_test.cpp
deleted file mode 100644
index 678e4b1..0000000
--- a/log-fourier_test.cpp
+++ /dev/null
@@ -1,50 +0,0 @@
-#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 426af40..a8d9778 100644
--- a/log_integrator.cpp
+++ b/log_integrator.cpp
@@ -92,7 +92,6 @@ 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);
@@ -109,7 +108,6 @@ 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));
@@ -144,47 +142,23 @@ int main(int argc, char* argv[]) {
for (unsigned i = 0; i < N; i++) {
dC[i] += -μₜ * Cₜ[i];
- 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));
+ Real ΓR;
+ for (unsigned j = 0; j < N; j++) {
+
}
- 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);
- 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<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));
}
+ std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+ std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
- μₜ /= 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];
- }
+ μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, x);
ΔCₜ = 0;
for (unsigned i = 0; i < N; i++) {
@@ -209,6 +183,8 @@ 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";
@@ -219,6 +195,8 @@ int main(int argc, char* argv[]) {
γ₀ /= 2;
Cₜ = Cₜ₋₁;
Rₜ = Rₜ₋₁;
+ Ĉₜ = Ĉₜ₋₁;
+ Ȓₜ = Ȓₜ₋₁;
μₜ = μₜ₋₁;
} else {
Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);
@@ -231,6 +209,8 @@ int main(int argc, char* argv[]) {
β += Δβ;
Cₜ₋₁ = Cₜ;
Rₜ₋₁ = Rₜ;
+ Ĉₜ₋₁ = Ĉₜ;
+ Ȓₜ₋₁ = Ȓₜ;
μₜ₋₁ = μₜ;
}
}