diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 1e6da0b8..75026e2a 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -138,9 +138,11 @@ protected: std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl; auto name = as[level].actions.at(a)->action_name(); if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); + DumpSliceNorm("force before Ta",force,Nd-1); force = FieldImplementation::projectForce(force); // Ta for gauge fields double end_force = usecond(); + DumpSliceNorm("force before filter",force,Nd-1); MomFilter->applyFilter(force); Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm. nb. norm2(latt) = \sum_x norm2(latt[x]) @@ -162,20 +164,12 @@ protected: double time_force = (end_force - start_force) / 1e3; std::cout << GridLogMessage << "["<(force,Nd-1); - DumpSliceNorm("force_t",pol); - pol=Zero(); - PokeIndex(force,pol,Nd-1); - DumpSliceNorm("force_xyz",force); - */ + DumpSliceNorm("force after filter",force,Nd-1); } // Force from the other representations as[level].apply(update_P_hireps, Representations, Mom, U, ep); - MomFilter->applyFilter(Mom); } void update_U(Field& U, double ep) @@ -189,8 +183,12 @@ protected: void update_U(MomentaField& Mom, Field& U, double ep) { + MomentaField MomFiltered(Mom.Grid()); + MomFiltered = Mom; + MomFilter->applyFilter(MomFiltered); + // exponential of Mom*U in the gauge fields case - FieldImplementation::update_field(Mom, U, ep); + FieldImplementation::update_field(MomFiltered, U, ep); // Update the smeared fields, can be implemented as observer Smearer.set_Field(U); @@ -366,7 +364,6 @@ public: as[level].apply(refresh_hireps, Representations, sRNG, pRNG); } - MomFilter->applyFilter(P); } // to be used by the actionlevel class to iterate diff --git a/HMC/Mobius2f.cc b/HMC/Mobius2f.cc index 78b60775..3383e97c 100644 --- a/HMC/Mobius2f.cc +++ b/HMC/Mobius2f.cc @@ -56,12 +56,12 @@ int main(int argc, char **argv) { MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 0; + HMCparams.StartTrajectory = 17; HMCparams.Trajectories = 200; - HMCparams.NoMetropolisUntil= 20; + HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - HMCparams.StartingType =std::string("ColdStart"); - // HMCparams.StartingType =std::string("CheckpointStart"); + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); @@ -94,7 +94,7 @@ int main(int argc, char **argv) { RealD b = 1.0; RealD c = 0.0; - std::vector hasenbusch({ 0.1 }); + std::vector hasenbusch({ 0.1, 0.4, 0.7 }); auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); diff --git a/HMC/Mobius2f_DDHMC_mixed.cc b/HMC/Mobius2f_DDHMC_mixed.cc index 3a762bf0..5d48d1bb 100644 --- a/HMC/Mobius2f_DDHMC_mixed.cc +++ b/HMC/Mobius2f_DDHMC_mixed.cc @@ -51,7 +51,7 @@ public: {}; virtual void refreshRestrict(FermionField &eta) { - Domains.ProjectDomain(eta,1); + Domains.ProjectDomain(eta,0); DumpSliceNorm("refresh Restrict eta",eta); }; }; @@ -97,20 +97,29 @@ int main(int argc, char **argv) // typedef GenericHMCRunner HMCWrapper; // MD.name = std::string("Force Gradient"); typedef GenericHMCRunner HMCWrapper; + /* MD.name = std::string("MinimumNorm2"); MD.MDsteps = 4; // dH = 0.08 // MD.MDsteps = 3; // dH = 0.8 MD.trajL = 1.0; + */ HMCparameters HMCparams; - HMCparams.StartTrajectory = 48; - HMCparams.Trajectories = 20; + { + XmlReader HMCrd("HMCparameters.xml"); + read(HMCrd,"HMCparameters",HMCparams); + std::cout << GridLogMessage<< HMCparams < LinOpD; std::vector LinOpF; - const int MX_inner = 1000; - const RealD MX_tol = 1.0e-8; + int MX_inner = 1000; + RealD MX_tol = 1.0e-5; for(int h=0;h (BoundaryDenominator, @@ -344,13 +357,13 @@ int main(int argc, char **argv) DomainDecomposedBoundaryTwoFlavourBosonPseudoFermion (BoundaryNumerator, BoundaryDerivativeStoppingCondition,ActionStoppingCondition,MX_tol)); - /* +#else Level1.push_back(new DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion (BoundaryNumerator, BoundaryDenominator, BoundaryDerivativeStoppingCondition,ActionStoppingCondition)); - */ +#endif ///////////////////////////////////////////////////////////// // Gauge action diff --git a/HMC/Mobius2p1fEOFA.cc b/HMC/Mobius2p1fEOFA.cc index 68ecc355..5239ab19 100644 --- a/HMC/Mobius2p1fEOFA.cc +++ b/HMC/Mobius2p1fEOFA.cc @@ -59,20 +59,20 @@ int main(int argc, char **argv) { IntegratorParameters MD; // typedef GenericHMCRunner HMCWrapper; // MD.name = std::string("Leap Frog"); - // typedef GenericHMCRunner HMCWrapper; - // MD.name = std::string("Force Gradient"); - typedef GenericHMCRunner HMCWrapper; - MD.name = std::string("MinimumNorm2"); - MD.MDsteps = 12; + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("Force Gradient"); + //typedef GenericHMCRunner HMCWrapper; + //MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 15; MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 211; + HMCparams.StartTrajectory = 0; HMCparams.Trajectories = 1000; - HMCparams.NoMetropolisUntil= 0; + HMCparams.NoMetropolisUntil= 10; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - // HMCparams.StartingType =std::string("ColdStart"); - HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.StartingType =std::string("ColdStart"); + //HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); @@ -97,16 +97,16 @@ int main(int argc, char **argv) { TheHMC.Resources.AddObservable(); ////////////////////////////////////////////// - const int Ls = 16; + const int Ls = 24; Real beta = 2.13; - Real light_mass = 0.01; - Real strange_mass = 0.04; + Real light_mass = 0.005; + Real strange_mass = 0.0362; Real pv_mass = 1.0; RealD M5 = 1.8; - RealD b = 1.0; - RealD c = 0.0; + RealD b = 1.5; + RealD c = 0.5; - std::vector hasenbusch({ 0.1, 0.3, 0.6 }); + std::vector hasenbusch({ 0.02, 0.2, 0.6 }); auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian();