summaryrefslogtreecommitdiff
path: root/fourier_integrator.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'fourier_integrator.cpp')
-rw-r--r--fourier_integrator.cpp102
1 files changed, 77 insertions, 25 deletions
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 4f4eff5..9f8db97 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -8,9 +8,9 @@ int main(int argc, char* argv[]) {
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
- Real y₀ = 0;
- Real yₘₐₓ = 0.5;
- Real Δy = 0.05;
+ Real β₀ = 0;
+ Real βₘₐₓ = 0.5;
+ Real Δβ = 0.05;
unsigned log2n = 8;
Real τₘₐₓ = 20;
@@ -41,13 +41,13 @@ int main(int argc, char* argv[]) {
τ₀ = atof(optarg);
break;
case '0':
- y₀ = atof(optarg);
+ β₀ = atof(optarg);
break;
case 'y':
- yₘₐₓ = atof(optarg);
+ βₘₐₓ = atof(optarg);
break;
case 'd':
- Δy = atof(optarg);
+ Δβ = atof(optarg);
break;
case 'I':
maxIterations = (unsigned)atof(optarg);
@@ -68,8 +68,11 @@ int main(int argc, char* argv[]) {
Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
- Real z = 0.5;
- Real Γ₀ = 1 + τ₀ / 2;
+ Real Γ₀ = 1.0;
+ Real μ = Γ₀;
+ if (τ₀ > 0) {
+ μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+ }
std::vector<Real> C(2 * n);
std::vector<Real> R(2 * n);
@@ -78,42 +81,52 @@ int main(int argc, char* argv[]) {
std::vector<Complex> Ct;
std::vector<Complex> Rt;
- Real y = y₀;
+ Real β = β₀;
if (!loadData) {
// start from the exact solution for τ₀ = 0
for (unsigned i = 0; i < n; i++) {
Real τ = i * Δτ * M_PI;
if (τ₀ > 0) {
- C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+ C[i] = Γ₀ * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
} else {
- C[i] = Γ₀ / 2 * exp(-z * τ) / z;
+ C[i] = Γ₀ * exp(-μ * τ) / μ;
}
if (i > 0) {
C[2 * n - i] = C[i];
}
- R[i] = exp(-z * τ);
+ R[i] = exp(-μ * τ);
}
Ct = fft.fourier(C);
Rt = fft.fourier(R);
} else {
- std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+ std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
cfile.close();
for (unsigned i = 1; i < n; i++) {
C[2 * n - i] = C[i];
}
- std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+ std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
rfile.close();
Ct = fft.fourier(C);
Rt = fft.fourier(R);
- z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
+ μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);
}
- while (y += Δy, y <= yₘₐₓ) {
+ std::vector<Real> Cb = C;
+ std::vector<Real> Rb = R;
+ std::vector<Complex> Ctb = Ct;
+ std::vector<Complex> Rtb = Rt;
+
+ Real fac = 1;
+ while (β += Δβ, β <= βₘₐₓ) {
+ Real Δμ = 1e-2;
+ Real μ₁ = 0;
+ Real μ₂ = 0;
+ while (true) {
Real ΔC = 1;
Real ΔCprev = 1000;
unsigned it = 0;
@@ -123,12 +136,12 @@ int main(int argc, char* argv[]) {
for (unsigned i = 0; i < Rt.size(); i++) {
Real ω = i * Δω;
- Rt[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
+ Rt[i] = (1.0 + pow(β, 2) * RddfCt[i] * Rt[i]) / (μ + 1i * ω);
}
for (unsigned i = 0; i < Ct.size(); i++) {
Real ω = i * Δω;
- Ct[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
+ Ct[i] = (2 * Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(β, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (μ + 1i * ω);
}
std::vector<Real> Cnew = fft.inverse(Ct);
@@ -151,8 +164,7 @@ int main(int argc, char* argv[]) {
R[i] += γ * (Rnew[i] - R[i]);
}
- z *= Cnew[0];
-
+ /*
if (it % maxIterations == 0) {
if (ΔC < ΔCprev) {
γ = std::min(1.1 * γ, 1.0);
@@ -162,20 +174,60 @@ int main(int argc, char* argv[]) {
ΔCprev = ΔC;
}
+ */
- std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
+ std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];
+ std::cerr << "\r";
}
+ if (std::isnan(C[0])) {
+ C = Cb;
+ R = Rb;
+ Ct = Ctb;
+ Rt = Rtb;
+ μ /= sqrt(sqrt(fac*std::tanh(Cb[0]-1)+1));
+ fac /= 2;
+ μ₁ = 0;
+ μ₂ = 0;
+ } else {
+ Cb = C;
+ Rb = R;
+ Ctb = Ct;
+ Rtb = Rt;
+ if (pow(C[0] - 1, 2) < ε) {
+ break;
+ }
+ if (μ₁ == 0 || μ₂ == 0) {
+ if (C[0] > 1 && μ₁ == 0) {
+ /* We found a lower bound */
+ μ₁ = μ;
+ }
+ if (C[0] < 1 && μ₂ == 0) {
+ /* We found an upper bound */
+ μ₂ = μ;
+ }
+ μ *= sqrt(sqrt(fac*std::tanh(C[0]-1)+1));
+ } else {
+ /* Once the bounds are set, we can use bisection */
+ if (C[0] > 1) {
+ μ₁ = μ;
+ } else {
+ μ₂ = μ;
+ }
+ μ = (μ₁ + μ₂) / 2;
+ }
+ }
+ }
- Real e = energy(C, R, p, s, λ, y, Δτ);
+ Real e = energy(C, R, p, s, λ, β, Δτ);
- std::cerr << "y " << y << " " << e << " " << z << std::endl;
+ std::cerr << "y " << β << " " << e << " " << μ << std::endl;
- std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+ std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real));
outfile.close();
- std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+ std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
outfileR.close();
}