From 089f0ab582071b76769a8f2be7697922b016fe70 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 28 Jul 2016 16:44:41 +0100 Subject: [PATCH] Debugged HMC for Creutz relation --- lib/qcd/action/ActionBase.h | 4 ++-- lib/qcd/action/fermion/WilsonKernels.h | 18 +++++++-------- lib/qcd/action/pseudofermion/TwoFlavour.h | 22 ++++++++++++++++++ lib/qcd/hmc/HmcRunner.h | 11 +++++---- lib/qcd/hmc/integrators/Integrator.h | 15 ++++++++----- lib/qcd/representations/adjoint.h | 15 ++++++++++++- lib/qcd/smearing/GaugeConfiguration.h | 25 +++++++++++++++++++++ lib/qcd/smearing/StoutSmearing.h | 6 ++--- lib/qcd/utils/SUn.h | 6 ++--- lib/qcd/utils/SUnAdjoint.h | 4 ++-- lib/tensors/Tensor_exp.h | 2 +- tests/Test_hmc_WilsonAdjointFermionGauge.cc | 7 +++--- 12 files changed, 100 insertions(+), 35 deletions(-) diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h index 5c9ca6b5..ad8d20ae 100644 --- a/lib/qcd/action/ActionBase.h +++ b/lib/qcd/action/ActionBase.h @@ -150,11 +150,11 @@ struct ActionLevelHirep { // Loop on tuple for a callable function template inline typename std::enable_if::value, void>::type apply( - Callable, Repr& R,Args...) const {} + Callable, Repr& R,Args&...) const {} template inline typename std::enable_if::value, void>::type apply( - Callable fn, Repr& R, Args... arguments) const { + Callable fn, Repr& R, Args&... arguments) const { fn(std::get(actions_hirep), std::get(R.rep), arguments...); apply(fn, R, arguments...); } diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 59c71bfd..23e590b9 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -54,7 +54,7 @@ class WilsonKernels : public FermionOperator, public WilsonKernelsStatic { public: template - typename std::enable_if::type + typename std::enable_if::type DiracOptDhopSite( StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, std::vector > &buf, @@ -85,7 +85,7 @@ class WilsonKernels : public FermionOperator, public WilsonKernelsStatic { } template - typename std::enable_if::type + typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type DiracOptDhopSite( StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, std::vector > &buf, @@ -102,7 +102,7 @@ class WilsonKernels : public FermionOperator, public WilsonKernelsStatic { } template - typename std::enable_if::type + typename std::enable_if::type DiracOptDhopSiteDag( StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, std::vector > &buf, @@ -126,11 +126,11 @@ class WilsonKernels : public FermionOperator, public WilsonKernelsStatic { } template - typename std::enable_if::type - DiracOptDhopSiteDag( - StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, - std::vector > &buf, - int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) { + typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type + DiracOptDhopSiteDag( + StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, + std::vector > &buf, + int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) { for (int site = 0; site < Ns; site++) { for (int s = 0; s < Ls; s++) { WilsonKernels::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF, sU, @@ -140,7 +140,7 @@ class WilsonKernels : public FermionOperator, public WilsonKernelsStatic { sU++; } } - + void DiracOptDhopDir( StencilImpl &st, DoubledGaugeField &U, std::vector > &buf, diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index 3595dc67..01e56daf 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -133,6 +133,28 @@ class TwoFlavourPseudoFermionAction : public Action { DerivativeSolver(MdagMOp, Phi, X); MdagMOp.Op(X, Y); + // Check hermiticity + /* + std::vector seeds({1,2,3,4}); + GridParallelRNG RNG(U._grid); RNG.SeedFixedIntegers(seeds); + FermionField RNGphi(FermOp.FermionGrid()); + FermionField RNGchi(FermOp.FermionGrid()); + FermionField Achi(FermOp.FermionGrid()); + FermionField Aphi(FermOp.FermionGrid()); + + random(RNG, RNGphi); + random(RNG, RNGchi); + MdagMOp.Op(RNGchi, Achi); + MdagMOp.Op(RNGphi, Aphi); + ComplexD pAc = innerProduct(RNGphi, Achi); + ComplexD cAp = innerProduct(RNGchi, Aphi); + //these should be real + ComplexD pAp = innerProduct(RNGphi, Aphi); + ComplexD cAc = innerProduct(RNGchi, Achi); + + std::cout< ParSeed({6, 7, 8, 9, 10}); // Create integrator, including the smearing policy - // Smearing policy + // Smearing policy, only defined for Nc=3 + std::cout << GridLogDebug << " Creating the Stout class\n"; double rho = 0.1; // smearing parameter, now hardcoded int Nsmear = 1; // number of smearing levels @@ -110,7 +111,9 @@ class NerscHmcRunnerTemplate { std::cout << GridLogDebug << " Creating the SmearedConfiguration class\n"; SmearedConfiguration SmearingPolicy(UGrid, Nsmear, Stout); std::cout << GridLogDebug << " done\n"; + ////////////// + //NoSmearing SmearingPolicy; typedef MinimumNorm2, RepresentationsPolicy > IntegratorType; // change here to change the algorithm IntegratorParameters MDpar(20, 1.0); @@ -131,19 +134,19 @@ class NerscHmcRunnerTemplate { HMCpar.MetropolisTest = true; sRNG.SeedFixedIntegers(SerSeed); pRNG.SeedFixedIntegers(ParSeed); - SU3::HotConfiguration(pRNG, U); + SU::HotConfiguration(pRNG, U); } else if (StartType == ColdStart) { // Cold start HMCpar.MetropolisTest = true; sRNG.SeedFixedIntegers(SerSeed); pRNG.SeedFixedIntegers(ParSeed); - SU3::ColdConfiguration(pRNG, U); + SU::ColdConfiguration(pRNG, U); } else if (StartType == TepidStart) { // Tepid start HMCpar.MetropolisTest = true; sRNG.SeedFixedIntegers(SerSeed); pRNG.SeedFixedIntegers(ParSeed); - SU3::TepidConfiguration(pRNG, U); + SU::TepidConfiguration(pRNG, U); } else if (StartType == CheckpointStart) { HMCpar.MetropolisTest = true; // CheckpointRestart diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index 67725621..80ce9e53 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -149,7 +149,7 @@ class Integrator { } // Force from the other representations - as[level].apply(update_P_hireps, Representations, Mom, U, ep); + //as[level].apply(update_P_hireps, Representations, Mom, U, ep); } void update_U(GaugeField& U, double ep) { @@ -173,7 +173,7 @@ class Integrator { // Update the smeared fields, can be implemented as observer Smearer.set_GaugeField(U); // Update the higher representations fields - Representations.update(U); // void functions if fundamental representation + //Representations.update(U); // void functions if fundamental representation } virtual void step(GaugeField& U, int level, int first, int last) = 0; @@ -203,6 +203,7 @@ class Integrator { GridParallelRNG& pRNG) { for (int a = 0; a < repr_set.size(); ++a) repr_set.at(a)->refresh(Rep.U, pRNG); + std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl; } } refresh_hireps{}; @@ -216,7 +217,7 @@ class Integrator { // of the Metropolis Smearer.set_GaugeField(U); // Set the (eventual) representations gauge fields - // Representations.update(U); + //Representations.update(U); // The Smearer is attached to a pointer of the gauge field // automatically gets the correct field @@ -230,7 +231,8 @@ class Integrator { as[level].actions.at(actionID)->refresh(Us, pRNG); } - as[level].apply(refresh_hireps, Representations, pRNG); + // Refresh the higher representation actions + //as[level].apply(refresh_hireps, Representations, pRNG); } } @@ -240,12 +242,13 @@ class Integrator { template void operator()(std::vector*> repr_set, Repr& Rep, int level, RealD& H) { - RealD H_hirep = 0.0; + for (int a = 0; a < repr_set.size(); ++a) { RealD Hterm = repr_set.at(a)->S(Rep.U); std::cout << GridLogMessage << "S Level " << level << " term " << a << " H Hirep = " << Hterm << std::endl; H += Hterm; + } } } S_hireps{}; @@ -278,7 +281,7 @@ class Integrator { << actionID << " H = " << Hterm << std::endl; H += Hterm; } - as[level].apply(S_hireps, Representations, level, H); + //as[level].apply(S_hireps, Representations, level, H); } return H; diff --git a/lib/qcd/representations/adjoint.h b/lib/qcd/representations/adjoint.h index d873353f..34cf2a42 100644 --- a/lib/qcd/representations/adjoint.h +++ b/lib/qcd/representations/adjoint.h @@ -28,6 +28,7 @@ class AdjointRep { explicit AdjointRep(GridBase *grid) : U(grid) {} void update_representation(const LatticeGaugeField &Uin) { + std::cout << GridLogDebug << "Updating adjoint representation\n" ; // Uin is in the fundamental representation // get the U in AdjointRep // (U_adj)_B = tr[e^a U e^b U^dag] @@ -49,9 +50,18 @@ class AdjointRep { auto U_mu = peekLorentz(U, mu); for (int a = 0; a < Dimension; a++) { tmp = 2.0 * adj(Uin_mu) * ta[a] * Uin_mu; - for (int b = 0; b < (ncolour * ncolour - 1); b++) + for (int b = 0; b < Dimension; b++) pokeColour(U_mu, trace(tmp * ta[b]), a, b); } + // Check matrix U_mu, must be real orthogonal + //reality + LatticeMatrix Ucheck = U_mu - conjugate(U_mu); + std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) << std::endl; + LatticeMatrix uno(Uin._grid); + uno = 1.0; + Ucheck = U_mu * adj(U_mu) - uno; + std::cout << GridLogMessage << "orthogonality check: " << norm2(Ucheck) << std::endl; + pokeLorentz(U, U_mu, mu); } } @@ -59,6 +69,7 @@ class AdjointRep { LatticeGaugeField RtoFundamentalProject(const LatticeField &in, Real scale = 1.0) const { LatticeGaugeField out(in._grid); + out = zero; for (int mu = 0; mu < Nd; mu++) { LatticeColourMatrix out_mu(in._grid); // fundamental representation @@ -70,6 +81,8 @@ class AdjointRep { projectOnAlgebra(h, in_mu, scale); FundamentalLieAlgebraMatrix(h, out_mu, 1.0); // apply scale only once pokeLorentz(out, out_mu, mu); + // Returns traceless antihermitian matrix Nc * Nc. + // Confirmed } return out; } diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h index 6219702b..57cccb0a 100644 --- a/lib/qcd/smearing/GaugeConfiguration.h +++ b/lib/qcd/smearing/GaugeConfiguration.h @@ -10,6 +10,29 @@ namespace Grid { namespace QCD { + //trivial class for no smearing + template< class Gimpl > +class NoSmearing { +public: + INHERIT_GIMPL_TYPES(Gimpl); + + GaugeField* + ThinLinks; + + NoSmearing(): ThinLinks(NULL) {} + + void set_GaugeField(GaugeField& U) { ThinLinks = &U; } + + void smeared_force(GaugeField& SigmaTilde) const {} + + GaugeField& get_SmearedU() { return *ThinLinks; } + + GaugeField& get_U(bool smeared = false) { + return *ThinLinks; + } + +}; + /*! @brief Smeared configuration container @@ -201,6 +224,8 @@ class SmearedConfiguration { SmearedConfiguration() : smearingLevels(0), StoutSmearing(), SmearedSet(), ThinLinks(NULL) {} + + // attach the smeared routines to the thin links U and fill the smeared set void set_GaugeField(GaugeField& U) { fill_smearedSet(U); } diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index 50a09972..b50ffe21 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -18,14 +18,12 @@ class Smear_Stout : public Smear { INHERIT_GIMPL_TYPES(Gimpl) Smear_Stout(Smear* base) : SmearBase(base) { - static_assert(Nc == 3, - "Stout smearing currently implemented only for Nc==3"); + assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3"); } /*! Default constructor */ Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE(rho)) { - static_assert(Nc == 3, - "Stout smearing currently implemented only for Nc==3"); + assert(Nc == 3);// "Stout smearing currently implemented only for Nc==3"); } ~Smear_Stout() {} // delete SmearBase... diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index c3e8295f..26c7f32d 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -48,7 +48,7 @@ class SU { using iSU2Matrix = iScalar > >; template using iSUnAlgebraVector = - iScalar > >; + iScalar > >; ////////////////////////////////////////////////////////////////////////////////////////////////// // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, @@ -656,7 +656,7 @@ class SU { } } - static void FundamentalLieAlgebraMatrix(LatticeAlgebraVector &h, + static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h, LatticeMatrix &out, Real scale = 1.0) { conformable(h, out); @@ -667,7 +667,7 @@ class SU { out = zero; for (int a = 0; a < AdjointDimension; a++) { generator(a, ta); - la = peekColour(h, a) * scale * ta; + la = peekColour(h, a) * timesI(ta) * scale; out += la; } } diff --git a/lib/qcd/utils/SUnAdjoint.h b/lib/qcd/utils/SUnAdjoint.h index 93a7601e..70f22f4c 100644 --- a/lib/qcd/utils/SUnAdjoint.h +++ b/lib/qcd/utils/SUnAdjoint.h @@ -118,7 +118,7 @@ class SU_Adjoint : public SU { for (int a = 0; a < Dimension; a++) { generator(a, iTa); - auto tmp = real(trace(iTa * in)) * scale; + auto tmp = - 2.0 * real(trace(iTa * in)) * scale; pokeColour(h_out, tmp, a); } } @@ -135,7 +135,7 @@ class SU_Adjoint : public SU { } for (int a = 0; a < Dimension; a++) { - auto tmp = real(trace(iTa[a] * in)) * scale; + auto tmp = - 2.0*real(trace(iTa[a] * in)) * scale; pokeColour(h_out, tmp, a); } } diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index 058e524d..ad0d2071 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -56,7 +56,7 @@ namespace Grid { temp = unit + temp*arg; } - return ProjectOnGroup(temp);//maybe not strictly necessary + return temp; } diff --git a/tests/Test_hmc_WilsonAdjointFermionGauge.cc b/tests/Test_hmc_WilsonAdjointFermionGauge.cc index 0de43427..e5125756 100644 --- a/tests/Test_hmc_WilsonAdjointFermionGauge.cc +++ b/tests/Test_hmc_WilsonAdjointFermionGauge.cc @@ -39,7 +39,8 @@ namespace Grid { namespace QCD { // Here change the allowed (higher) representations -typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations; +//typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations; +typedef Representations< FundamentalRepresentation > TheRepresentations; class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > { @@ -61,7 +62,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > { // temporarily need a gauge field LatticeGaugeField U(UGrid); - //AdjointRepresentation::LatticeField Ua(UGrid); + //AdjointRepresentation::LatticeField U(UGrid); // Gauge action WilsonGaugeActionR Waction(5.6); @@ -69,7 +70,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > { Real mass = -0.77; FermionAction FermOp(U, *FGrid, *FrbGrid, mass); - ConjugateGradient CG(1.0e-6, 10000); + ConjugateGradient CG(1.0e-8, 10000); TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG);