From 70101b6e893f2260216a8e5a1ea41ec22d118110 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 19 Apr 2025 23:35:26 -0300 Subject: Removed external μ control MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- fourier_integrator.cpp | 72 +++++++------------- log-fourier_integrator.cpp | 159 ++++++++++++++++++--------------------------- 2 files changed, 86 insertions(+), 145 deletions(-) diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp index 9ce47b6..96fbffd 100644 --- a/fourier_integrator.cpp +++ b/fourier_integrator.cpp @@ -120,13 +120,10 @@ int main(int argc, char* argv[]) { std::vector Rb = R; std::vector Ctb = Ct; std::vector Rtb = Rt; + Real μb = μ; Real fac = 1; - while (β += Δβ, β <= βₘₐₓ) { - Real Δμ = 1e-2; - Real μ₁ = 0; - Real μ₂ = 0; - while (true) { + while (β <= βₘₐₓ) { Real ΔC = 1; Real ΔCprev = 1000; unsigned it = 0; @@ -150,6 +147,8 @@ int main(int argc, char* argv[]) { Rnew[i] = 0; } + μ *= pow(tanh(C[0]-1)+1, 0.05); + ΔC = 0; for (unsigned i = 0; i < Cnew.size() / 2; i++) { ΔC += pow(Cnew[i] - C[i], 2); @@ -164,7 +163,6 @@ int main(int argc, char* argv[]) { R[i] += γ * (Rnew[i] - R[i]); } - /* if (it % maxIterations == 0) { if (ΔC < ΔCprev) { γ = std::min(1.1 * γ, 1.0); @@ -174,52 +172,23 @@ int main(int argc, char* argv[]) { ΔCprev = ΔC; } - */ + std::cerr << "\x1b[2K" << "\r"; 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)); - μ *= 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; - } - } } + if (std::isnan(C[0])) { + Δβ /= 2; + β -= Δβ; + C = Cb; + R = Rb; + Ct = Ctb; + Rt = Rtb; + μ = μb; + } else { + + std::cerr << "\x1b[2K" << "\r"; + Real e = energy(C, R, p, s, λ, β, Δτ); std::cerr << "y " << β << " " << e << " " << μ << std::endl; @@ -231,6 +200,13 @@ int main(int argc, char* argv[]) { 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(); + β += Δβ; + Cb = C; + Rb = R; + Rtb = Rt; + Ctb = Ct; + μb = μ; + } } return 0; diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index b15808b..531565c 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -93,108 +93,73 @@ int main(int argc, char* argv[]) { Real fac = 1; Real β = 0; while (β < βₘₐₓ) { - Real μ₁ = 0; - Real μ₂ = 0; - while (true) { - Real ΔC = 100; - Real ΔC₀ = 100; - unsigned it = 0; - while (ΔC > ε) { - std::vector dfC(N); - std::vector RddfC(N); - for (unsigned n = 0; n < N; n++) { - RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]); - dfC[n] = df(λ, p, s, Cₜ[n]); - } - std::vector RddfCt = fft.fourier(RddfC, false); - std::vector dfCt = fft.fourier(dfC, true); - - std::vector Ȓₜ₊₁(N); - std::vector Ĉₜ₊₁(N); - - for (unsigned n = 0; n < N; n++) { - Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n)); + Real ΔC = 100; + Real ΔC₀ = 100; + unsigned it = 0; + while (ΔC > ε) { + std::vector dfC(N); + std::vector RddfC(N); + for (unsigned n = 0; n < N; n++) { + RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]); + dfC[n] = df(λ, p, s, Cₜ[n]); + } + std::vector RddfCt = fft.fourier(RddfC, false); + std::vector dfCt = fft.fourier(dfC, true); + + std::vector Ȓₜ₊₁(N); + std::vector Ĉₜ₊₁(N); + + for (unsigned n = 0; n < N; n++) { + Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n)); // Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n); - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n)); - } - - std::vector Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); - - for (unsigned n = 0; n < N; n++) { - RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]); - } - RddfCt = fft.fourier(RddfC, false); - for (unsigned n = 0; n < N; n++) { - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n)); - } - std::vector Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); - - ΔC = 0; - for (unsigned i = 0; i < N; i++) { - ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]); - ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]); - } - ΔC = sqrt(ΔC) / (2*N); - - 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]); - } - - /* - if (ΔC < ΔC₀) { - ΔC₀ = ΔC; - it = 0; - γ = std::min(1.001 * γ, 1.0); - } else { - it++; - } - - if (it > 100) { - γ = std::max(0.5 * γ, 1e-3); - it = 0; - ΔC₀ = ΔC; - } - */ - - std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0]; - std::cerr << "\r"; + Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n)); } - if (std::isnan(Cₜ[0])) { - Cₜ = Cₜ₋₁; - Rₜ = Rₜ₋₁; - Ĉₜ = Ĉₜ₋₁; - Ȓₜ = Ȓₜ₋₁; - μ *= 1.1; - fac /= 2; - μ₁ = 0; - μ₂ = 0; - } else { - if (pow(Cₜ[0] - 1, 2) < ε) { - break; + + std::vector Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); + + for (unsigned n = 0; n < N; n++) { + RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]); } - if (μ₁ == 0 || μ₂ == 0) { - if (Cₜ[0] > 1 && μ₁ == 0) { - /* We found a lower bound */ - μ₁ = μ; - } - if (Cₜ[0] < 1 && μ₂ == 0) { - /* We found an upper bound */ - μ₂ = μ; - } - μ *= fac*std::tanh(Cₜ[0]-1)+1; + RddfCt = fft.fourier(RddfC, false); + for (unsigned n = 0; n < N; n++) { + Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n)); + } + std::vector Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); + + μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05); + + ΔC = 0; + for (unsigned i = 0; i < N; i++) { + ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]); + ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]); + } + ΔC = sqrt(ΔC) / (2*N); + + 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]); + } + + /* + if (ΔC < ΔC₀) { + ΔC₀ = ΔC; + it = 0; + γ = std::min(1.001 * γ, 1.0); } else { - /* Once the bounds are set, we can use bisection */ - if (Cₜ[0] > 1) { - μ₁ = μ; - } else { - μ₂ = μ; - } - μ = (μ₁ + μ₂) / 2; + it++; } + + if (it > 100) { + γ = std::max(0.5 * γ, 1e-3); + it = 0; + ΔC₀ = ΔC; } + */ + + std::cerr << "\x1b[2K" << "\r"; + std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0]; } /* Integrate the energy using Simpson's rule */ @@ -213,9 +178,9 @@ int main(int argc, char* argv[]) { } E *= β; + std::cerr << "\x1b[2K" << "\r"; std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl; β += Δβ; - μ *= 2; Cₜ₋₁ = Cₜ; Rₜ₋₁ = Rₜ; Ĉₜ₋₁ = Ĉₜ; -- cgit v1.2.3-70-g09d2