From e19171523be4e66e59baa3408e6098a6695fb0f1 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Fri, 13 Oct 2023 17:57:56 +0300
Subject: [PATCH] FTHMC Status at lattice conference commit

---
 HMC/FTHMC2p1f.cc                    | 18 ++++++-----
 HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc |  2 +-
 HMC/Mobius2p1f_EOFA_96I_hmc.cc      | 46 +++++++++++++++--------------
 3 files changed, 35 insertions(+), 31 deletions(-)

diff --git a/HMC/FTHMC2p1f.cc b/HMC/FTHMC2p1f.cc
index dd824138..7d93d168 100644
--- a/HMC/FTHMC2p1f.cc
+++ b/HMC/FTHMC2p1f.cc
@@ -54,15 +54,16 @@ int main(int argc, char **argv)
   //  MD.name    = std::string("Force Gradient");
   typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
   MD.name    = std::string("MinimumNorm2");
-  MD.MDsteps = 12;
+  MD.MDsteps = 24;
   MD.trajL   = 1.0;
 
   HMCparameters HMCparams;
-  HMCparams.StartTrajectory  = 0;
+  HMCparams.StartTrajectory  = 104;
   HMCparams.Trajectories     = 200;
   HMCparams.NoMetropolisUntil=  20;
   // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
-  HMCparams.StartingType     =std::string("HotStart");
+  //  HMCparams.StartingType     =std::string("HotStart");
+  HMCparams.StartingType     =std::string("CheckpointStart");
   HMCparams.MD = MD;
   HMCWrapper TheHMC(HMCparams);
 
@@ -87,6 +88,7 @@ int main(int argc, char **argv)
   // here there is too much indirection
   typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
   TheHMC.Resources.AddObservable<PlaqObs>();
+
   //////////////////////////////////////////////
 
   const int Ls      = 16;
@@ -134,7 +136,6 @@ int main(int argc, char **argv)
   ////////////////////////////////////
   ActionLevel<HMCWrapper::Field> Level1(1);
   ActionLevel<HMCWrapper::Field> Level2(2);
-  ActionLevel<HMCWrapper::Field> Level3(4);
 
   ////////////////////////////////////
   // Strange action
@@ -191,7 +192,7 @@ int main(int argc, char **argv)
   Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
   SmearedConfigurationMasked<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, Nstep, Stout);
   JacobianAction<HMCWrapper::ImplPolicy> Jacobian(&SmearingPolicy);
-  if( ApplySmearing ) Level2.push_back(&Jacobian);
+  if( ApplySmearing ) Level1.push_back(&Jacobian);
   std::cout << GridLogMessage << " Built the Jacobian "<< std::endl;
 
 
@@ -200,7 +201,7 @@ int main(int argc, char **argv)
   /////////////////////////////////////////////////////////////
   //  GaugeAction.is_smeared = ApplySmearing;
   GaugeAction.is_smeared = true;
-  Level3.push_back(&GaugeAction);
+  Level2.push_back(&GaugeAction);
 
   std::cout << GridLogMessage << " ************************************************"<< std::endl;
   std::cout << GridLogMessage << " Action complete -- NO FERMIONS FOR NOW -- FIXME"<< std::endl;
@@ -210,10 +211,11 @@ int main(int argc, char **argv)
 
 
   std::cout << GridLogMessage << " Running the FT HMC "<< std::endl;
-
   TheHMC.TheAction.push_back(Level1);
   TheHMC.TheAction.push_back(Level2);
-  TheHMC.TheAction.push_back(Level3);
+
+  TheHMC.ReadCommandLine(argc,argv);  // params on CML or from param file
+  TheHMC.initializeGaugeFieldAndRNGs(U);
 
   TheHMC.Run(SmearingPolicy); // for smearing
 
diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc b/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc
index 24d8951e..83f20b92 100644
--- a/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc
+++ b/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc
@@ -343,7 +343,7 @@ int main(int argc, char **argv) {
   // Probably dominates the force - back to EOFA.
   OneFlavourRationalParams SFRp;
   SFRp.lo       = 0.1;
-  SFRp.hi       = 25.0;
+  SFRp.hi       = 30.0;
   SFRp.MaxIter  = 10000;
   SFRp.tolerance= 1.0e-5;
   SFRp.mdtolerance= 2.0e-4;
diff --git a/HMC/Mobius2p1f_EOFA_96I_hmc.cc b/HMC/Mobius2p1f_EOFA_96I_hmc.cc
index 3d674db4..91f0bd95 100644
--- a/HMC/Mobius2p1f_EOFA_96I_hmc.cc
+++ b/HMC/Mobius2p1f_EOFA_96I_hmc.cc
@@ -128,7 +128,7 @@ template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, c
       ////////////////////////////////////////////////////////////////////////////////////
       // Make a mixed precision conjugate gradient
       ////////////////////////////////////////////////////////////////////////////////////
-#if 1
+#if 0
       RealD delta=1.e-4;
       std::cout << GridLogMessage << "Calling reliable update Conjugate Gradient" <<std::endl;
       ConjugateGradientReliableUpdate<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations*MaxOuterIterations,delta,SinglePrecGrid5,LinOpF,LinOpD);
@@ -180,7 +180,7 @@ int main(int argc, char **argv) {
   // 4/2 => 0.6 dH
   // 3/3 => 0.8 dH .. depth 3, slower
   //MD.MDsteps =  4;
-  MD.MDsteps =  14;
+  MD.MDsteps =  12;
   MD.trajL   = 0.5;
 
   HMCparameters HMCparams;
@@ -204,7 +204,7 @@ int main(int argc, char **argv) {
   TheHMC.Resources.LoadNerscCheckpointer(CPparams);
   std::cout << "loaded NERSC checpointer"<<std::endl;
   RNGModuleParameters RNGpar;
-  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.serial_seeds = "1 2 3 4 5 6 7 8 9 10";
   RNGpar.parallel_seeds = "6 7 8 9 10";
   TheHMC.Resources.SetRNGSeeds(RNGpar);
 
@@ -218,15 +218,14 @@ int main(int argc, char **argv) {
   RealD M5  = 1.8;
   RealD b   = 1.5;
   RealD c   = 0.5;
-  Real beta         = 2.13;
+  RealD beta         = 2.13;
   //  Real light_mass   = 5.4e-4;
   Real light_mass     = 7.8e-4;
+  //  Real light_mass     = 7.8e-3;
   Real strange_mass = 0.0362;
   Real pv_mass      = 1.0;
-  //  std::vector<Real> hasenbusch({ 0.01, 0.045, 0.108, 0.25, 0.51 , pv_mass });
-  //  std::vector<Real> hasenbusch({ light_mass, 0.01, 0.045, 0.108, 0.25, 0.51 , pv_mass });
-  std::vector<Real> hasenbusch({ 0.005, 0.0145, 0.045, 0.108, 0.25, 0.51 }); // Updated
-  //  std::vector<Real> hasenbusch({ light_mass, 0.0145, 0.045, 0.108, 0.25, 0.51 , 0.75 , pv_mass });
+  std::vector<Real> hasenbusch({ 0.005, 0.0145, 0.045, 0.108, 0.25, 0.35 , 0.51, 0.6, 0.8 }); // Updated
+  //std::vector<Real> hasenbusch({ 0.0145, 0.045, 0.108, 0.25, 0.35 , 0.51, 0.6, 0.8 }); // Updated
 
   auto GridPtr   = TheHMC.Resources.GetCartesian();
   auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
@@ -277,20 +276,20 @@ int main(int argc, char **argv) {
 
   //  double StoppingCondition = 1e-14;
   //  double MDStoppingCondition = 1e-9;
-  double StoppingCondition = 1e-9;
-  double MDStoppingCondition = 1e-8;
-  double MDStoppingConditionLoose = 1e-8;
-  double MDStoppingConditionStrange = 1e-8;
-  double MaxCGIterations = 300000;
+  double StoppingCondition = 1e-14;
+  double MDStoppingCondition = 1e-9;
+  double MDStoppingConditionLoose = 1e-9;
+  double MDStoppingConditionStrange = 1e-9;
+  double MaxCGIterations = 50000;
   ConjugateGradient<FermionField>  CG(StoppingCondition,MaxCGIterations);
   ConjugateGradient<FermionField>  MDCG(MDStoppingCondition,MaxCGIterations);
 
   ////////////////////////////////////
   // Collect actions
   ////////////////////////////////////
-  //  ActionLevel<HMCWrapper::Field> Level1(1);
-  ActionLevel<HMCWrapper::Field> Level2(1);
-  ActionLevel<HMCWrapper::Field> Level3(15);
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  ActionLevel<HMCWrapper::Field> Level2(2);
+  ActionLevel<HMCWrapper::Field> Level3(4);
 
   ////////////////////////////////////
   // Strange action
@@ -300,11 +299,11 @@ int main(int argc, char **argv) {
 
   // Probably dominates the force - back to EOFA.
   OneFlavourRationalParams SFRp;
-  SFRp.lo       = 0.1;
+  SFRp.lo       = 0.8;
   SFRp.hi       = 30.0;
   SFRp.MaxIter  = 10000;
-  SFRp.tolerance= 1.0e-8;
-  SFRp.mdtolerance= 2.0e-6;
+  SFRp.tolerance= 1.0e-12;
+  SFRp.mdtolerance= 1.0e-9;
   SFRp.degree   = 10;
   SFRp.precision= 50;
   
@@ -355,8 +354,10 @@ int main(int argc, char **argv) {
   ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 
     EOFA(Strange_Op_L, Strange_Op_R, 
 	 ActionCG, 
-	 ActionCGL, ActionCGR,
-	 DerivativeCGL, DerivativeCGR,
+	 //	 ActionCGL, ActionCGR,
+	 //	 DerivativeCGL, DerivativeCGR,
+	 ActionCG, ActionCG,
+	 DerivativeCG, DerivativeCG,
 	 SFRp, true);
   Level2.push_back(&EOFA);
 
@@ -443,13 +444,14 @@ int main(int argc, char **argv) {
   }
   int nquo=Quotients.size();
   for(int h=0;h<nquo;h++){
-    Level2.push_back(Quotients[h]);
+    Level1.push_back(Quotients[h]);
   }
 
   /////////////////////////////////////////////////////////////
   // Gauge action
   /////////////////////////////////////////////////////////////
   Level3.push_back(&GaugeAction);
+  TheHMC.TheAction.push_back(Level1);
   TheHMC.TheAction.push_back(Level2);
   TheHMC.TheAction.push_back(Level3);
   std::cout << GridLogMessage << " Action complete "<< std::endl;