Adjust damping with failures

1 files changed, 17 insertions, 11 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 4ee9150..e2446fd 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -71,12 +71,13 @@ int main(int argc, char* argv[]) {
   unsigned log2n = 8;
   Real τₘₐₓ = 20 / M_PI;
 
-  unsigned maxIterations = 100;
+  unsigned maxIterations = 1000;
   Real ε = 1e-14;
+  Real γ = 0;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:2:T:t:y:d:I:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:2:T:t:y:d:I:g:")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
@@ -99,6 +100,9 @@ int main(int argc, char* argv[]) {
     case 'I':
       maxIterations = (unsigned)atof(optarg);
       break;
+    case 'g':
+      γ = atof(optarg);
+      break;
     default:
       exit(1);
     }
@@ -132,8 +136,9 @@ int main(int argc, char* argv[]) {
   Real y = 0;
 
   while (y += Δy, y <= yₘₐₓ) {
+    Real ΔC = 1;;
     unsigned it = 0;
-    while (true) {
+    while (sqrt(ΔC / C.size()) > ε) {
       it++;
       std::vector<Real> RddfC(C.size());
       for (unsigned i = 0; i < C.size(); i++) {
@@ -155,24 +160,28 @@ int main(int argc, char* argv[]) {
         Rtnew[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
       }
 
-      Rt = Rtnew;
+      for (unsigned i = 0; i < Rt.size(); i++) {
+        Rt[i] += γ * (Rtnew[i] - Rt[i]);
+      }
 
       for (unsigned i = 0; i < Rt.size(); i++) {
         Real ω = i * Δω;
         Ctnew[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
       }
 
-      Ct = Ctnew;
+      for (unsigned i = 0; i < Ct.size(); i++) {
+        Ct[i] += γ * (Ctnew[i] - Ct[i]);
+      }
 
       std::vector<Real> Cnew = fft.inverse(Ctnew);
       std::vector<Real> Rnew = fft.inverse(Rtnew);
 
-      Real ΔC = 0;
+      ΔC = 0;
       for (unsigned i = 0; i < Cnew.size(); i++) {
         ΔC += pow(Cnew[i] - C[i], 2);
       }
 
-      z *= sqrt(Cnew[0]);
+      z *= Cnew[0];
 
       Real energy = 0;
 
@@ -185,12 +194,9 @@ int main(int argc, char* argv[]) {
 
       std::cerr << "y " << y << " iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << pow(y, 2) / z << " " << energy << std::endl;
 
-      if (sqrt(ΔC / C.size()) < ε) break;
       if (it > maxIterations) {
-        y -= Δy;
-        Δy /= 2;
-        y += Δy;
         it = 0;
+        γ /= 2;
       }
     }
   }