fix test of the integrator

claudiopica · Jul 11, 2024 · 8adb063 · 8adb063
1 parent 5b90509
commit 8adb063
Showing 1 changed file with 49 additions and 51 deletions.
diff --git a/TestProgram/Integrators/check_integrator_1.c b/TestProgram/Integrators/check_integrator_1.c
@@ -26,7 +26,7 @@
 #include "spectrum.h"
 #include "setup.h"
 #include "linear_algebra.h"
-#define SCALING_RANGE 20
+#define SCALING_RANGE 5
 
 /* flow control variable */
 hmc_flow flow = init_hmc_flow(flow);
@@ -86,59 +86,57 @@ int main(int argc, char *argv[])
     lprintf("MAIN", 0, "Checking that the action doesn't change on reusing the integrator twice from the same starting fields %.2e\n", res);
     lprintf("MAIN", 0, "(should be around 1*10^(-15) or so)\n\n\n");
 
-    if (res > 1.0e-12)
+    if (res > 1.0e-11)
       return_value++;
   }
-  if (1 == 1)
+
+  double rt0[SCALING_RANGE], rt1[SCALING_RANGE], rt2[SCALING_RANGE];
+
+  for (int i = 0; i < SCALING_RANGE; i++)
+  {
+    set_integrator_type(&(flow.hmc_v->hmc_p), LEAPFROG);
+    set_first_integrator_nsteps(&(flow.hmc_v->hmc_p), i * 10 + 20);
+    gettimeofday(&start, 0);
+    rt0[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
+    gettimeofday(&end, 0);
+    timeval_subtract(&etime, &end, &start);
+    elapsed_lf = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
+
+    set_integrator_type(&(flow.hmc_v->hmc_p), O2MN);
+    gettimeofday(&start, 0);
+    rt1[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
+    gettimeofday(&end, 0);
+    timeval_subtract(&etime, &end, &start);
+    elapsed_o2 = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
+
+    set_integrator_type(&(flow.hmc_v->hmc_p), O4MN);
+    gettimeofday(&start, 0);
+    rt2[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
+    gettimeofday(&end, 0);
+    timeval_subtract(&etime, &end, &start);
+    elapsed_o4 = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
+
+    lprintf("MAIN", 0, "Delta H for first level monomial nsteps: %d dt: %lf LF: %1.16e %1.16e msec  O2: %1.16e %1.16e msec  O4: %1.16e %1.16e msec\n\n",
+            i * 10 + 20, 1.0 / ((double)(i * 10 + 20)),
+            fabs(rt0[i]), elapsed_lf,
+            fabs(rt1[i]), elapsed_o2,
+            fabs(rt2[i]), elapsed_o4);
+  }
+
+  double scaling_deviation1;
+  double scaling_deviation2;
+  double dt, ldt = 1.0 / ((SCALING_RANGE - 1) * 10 + 20);
+  for (int i = 0; i < SCALING_RANGE - 1; i++)
   {
-    double rt0[SCALING_RANGE], rt1[SCALING_RANGE], rt2[SCALING_RANGE];
-
-    for (int i = 0; i < SCALING_RANGE; i++)
-    {
-      set_integrator_type(&(flow.hmc_v->hmc_p), LEAPFROG);
-      set_first_integrator_nsteps(&(flow.hmc_v->hmc_p), i * 10 + 20);
-      gettimeofday(&start, 0);
-      rt0[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
-      gettimeofday(&end, 0);
-      timeval_subtract(&etime, &end, &start);
-      elapsed_lf = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
-
-      set_integrator_type(&(flow.hmc_v->hmc_p), O2MN);
-      gettimeofday(&start, 0);
-      rt1[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
-      gettimeofday(&end, 0);
-      timeval_subtract(&etime, &end, &start);
-      elapsed_o2 = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
-
-      set_integrator_type(&(flow.hmc_v->hmc_p), O4MN);
-      gettimeofday(&start, 0);
-      rt2[i] = integrate_ghmc(1, &(flow.hmc_v->hmc_p));
-      gettimeofday(&end, 0);
-      timeval_subtract(&etime, &end, &start);
-      elapsed_o4 = etime.tv_sec * 1000. + etime.tv_usec * 0.001;
-
-      lprintf("MAIN", 0, "Delta H for first level monomial nsteps: %d dt: %lf LF: %1.16e %1.16e msec  O2: %1.16e %1.16e msec  O4: %1.16e %1.16e msec\n\n",
-              i * 10 + 20, 1.0 / ((double)(i * 10 + 20)),
-              fabs(rt0[i]), elapsed_lf,
-              fabs(rt1[i]), elapsed_o2,
-              fabs(rt2[i]), elapsed_o4);
-    }
-
-    double scaling_deviation1;
-    double scaling_deviation2;
-    double dt, ldt = 1.0 / ((SCALING_RANGE - 1) * 10 + 20);
-    for (int i = 0; i < SCALING_RANGE - 1; i++)
-    {
-
-      dt = 1.0 / (i * 10 + 20);
-      scaling_deviation1 = fabs(1.0 - rt0[i] * rt1[SCALING_RANGE - 1] / rt1[i] / rt0[SCALING_RANGE - 1]);
-      scaling_deviation2 = fabs(1.0 - dt * dt * rt1[i] * rt2[SCALING_RANGE - 1] / rt2[i] / rt1[SCALING_RANGE - 1] / ldt / ldt);
-      lprintf("MAIN", 0, "Scaling deviation for LF/O2 (normalized at the smallest dt) at dt=%lf %1.16e \n(should be around 1*10^(-2) or so)\n\n", dt, scaling_deviation1);
-      lprintf("MAIN", 0, "Scaling deviation for dt^2*O2/O4 (normalized at the smallest dt) at dt=%lf %1.16e \n(should be around 1*10^(-2) or so)\n\n", dt, scaling_deviation2);
-
-      if (scaling_deviation1 > 0.5 || scaling_deviation2 > 0.5)
-        return_value++;
-    }
+
+    dt = 1.0 / (i * 10 + 20);
+    scaling_deviation1 = fabs(1.0 - rt0[i] * rt1[SCALING_RANGE - 1] / rt1[i] / rt0[SCALING_RANGE - 1]);
+    scaling_deviation2 = fabs(1.0 - dt * dt * rt1[i] * rt2[SCALING_RANGE - 1] / rt2[i] / rt1[SCALING_RANGE - 1] / ldt / ldt);
+    lprintf("MAIN", 0, "Scaling deviation for LF/O2 (normalized at the smallest dt) at dt=%lf %1.16e \n(should be around 1*10^(-2) or so)\n\n", dt, scaling_deviation1);
+    lprintf("MAIN", 0, "Scaling deviation for dt^2*O2/O4 (normalized at the smallest dt) at dt=%lf %1.16e \n(should be around 1*10^(-2) or so)\n\n", dt, scaling_deviation2);
+
+    if (scaling_deviation1 > 0.5 || scaling_deviation2 > 0.5)
+      return_value++;
   }
 
   /* close communications */