Added midMD checkpointing (for lattice only for now)

Grid/qcd/hmc/HMC.h

@@ -112,6 +112,8 @@ private:
   IntegratorType &TheIntegrator;
   ObsListType Observables;
 
+  int traj_num;
+
   /////////////////////////////////////////////////////////
   // Metropolis step
   /////////////////////////////////////////////////////////
@@ -202,7 +204,7 @@ private:
 
     std::cout << GridLogMessage << "--------------------------------------------------\n";
     std::cout << GridLogMessage << " Molecular Dynamics evolution ";
-    TheIntegrator.integrate(U);
+    TheIntegrator.integrate(U,traj_num);
     std::cout << GridLogMessage << "--------------------------------------------------\n";
 
     //////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -244,7 +246,7 @@ public:
   HybridMonteCarlo(HMCparameters _Pams, IntegratorType &_Int,
                    GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, 
                    ObsListType _Obs, Field &_U)
-    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
+    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U),traj_num(0) {}
   ~HybridMonteCarlo(){};
 
   void evolve(void) {
@@ -259,9 +261,10 @@ public:
     unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory;
 
     for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) {
+    
 
       std::cout << GridLogHMC << "-- # Trajectory = " << traj << "\n";
-
+      traj_num=traj;
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
       	std::cout << GridLogHMC << "-- Thermalization" << std::endl;
       }

Grid/qcd/hmc/integrators/Integrator.h

@@ -34,6 +34,7 @@ directory
 #define INTEGRATOR_INCLUDED
 
 #include <memory>
+#include <Grid/parallelIO/NerscIO.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -676,7 +677,7 @@ public:
   }
 
   
-  void integrate(Field& U) 
+  void integrate(Field& U, int traj=-1 ) 
   {
     // reset the clocks
     t_U = 0;
@@ -688,6 +689,12 @@ public:
       int first_step = (stp == 0);
       int last_step = (stp == Params.MDsteps - 1);
       this->step(U, 0, first_step, last_step);
+      if (traj>=0){
+        std::string file("./config."+std::to_string(traj)+"_"+std::to_string(stp+1) );
+        int precision32 = 0;
+        int tworow      = 0;
+        NerscIO::writeConfiguration(U,file,tworow,precision32);
+      }
     }
 
     // Check the clocks all match on all levels
@@ -697,7 +704,6 @@ public:
     }
 
     FieldImplementation::Project(U);
-
     // and that we indeed got to the end of the trajectory
     assert(fabs(t_U - Params.trajL) < 1.0e-6);
 

Grid/qcd/hmc/integrators/Integrator_algorithm.h

@@ -140,7 +140,7 @@ template <class FieldImplementation_, class SmearingPolicy, class Representation
 class MinimumNorm2 : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy> 
 {
 private:
-  const RealD lambda = 0.1931833275037836;
+//  const RealD lambda = 0.1931833275037836;
 
 public:
   typedef FieldImplementation_ FieldImplementation;
@@ -155,6 +155,11 @@ public:
     // level  : current level
     // fl     : final level
     // eps    : current step size
+    assert(level<3);
+    RealD lambda= this->Params.lambda0;
+    if (level>0) lambda= this->Params.lambda1;
+    if (level>1) lambda= this->Params.lambda2;
+    std::cout << GridLogMessage << "level: "<<level<< "lambda: "<<lambda<<std::endl;
 
     int fl = this->as.size() - 1;
 

benchmarks/Benchmark_ITT.cc

@@ -365,9 +365,15 @@ public:
     GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
     std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
 
+#if 1
     typedef DomainWallFermionF Action;
     typedef typename Action::FermionField Fermion;
     typedef LatticeGaugeFieldF Gauge;
+#else
+    typedef GparityDomainWallFermionF Action;
+    typedef typename Action::FermionField Fermion;
+    typedef LatticeGaugeFieldF Gauge;
+#endif
     
     ///////// Source preparation ////////////
     Gauge Umu(UGrid);  SU<Nc>::HotConfiguration(RNG4,Umu);