Automaticlly adjust γ

1 files changed, 13 insertions, 8 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 2fe7e49..9e41b5f 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -107,9 +107,10 @@ int main(int argc, char* argv[]) {
   }
 
   while (y += Δy, y <= yₘₐₓ) {
-    Real ΔC = 1;;
+    Real ΔC = 1;
+    Real ΔCprev = 1000;
     unsigned it = 0;
-    while (sqrt(ΔC / C.size()) > ε) {
+    while (sqrt(2 * ΔC / C.size()) > ε) {
       it++;
       auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
 
@@ -130,7 +131,7 @@ int main(int argc, char* argv[]) {
       }
 
       ΔC = 0;
-      for (unsigned i = 0; i < Cnew.size(); i++) {
+      for (unsigned i = 0; i < Cnew.size() / 2; i++) {
         ΔC += pow(Cnew[i] - C[i], 2);
         ΔC += pow(Rnew[i] - R[i], 2);
       }
@@ -139,18 +140,22 @@ int main(int argc, char* argv[]) {
         C[i] += γ * (Cnew[i] - C[i]);
       }
 
-      for (unsigned i = 0; i < Rnew.size(); i++) {
+      for (unsigned i = 0; i < Rnew.size() / 2; i++) {
         R[i] += γ * (Rnew[i] - R[i]);
       }
 
       z *= Cnew[0];
 
-      std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << sqrt(ΔC / C.size()) << std::endl;
-
-      if (it > maxIterations) {
-        it = 0;
+      if (ΔC < ΔCprev) {
+        γ = std::min(1.1 * γ, 1.0);
+      } else {
         γ /= 2;
       }
+
+      ΔCprev = ΔC;
+
+      std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
+
     }
 
     Real e = energy(C, R, p, s, λ, y, Δτ);