From 147e2025b991bde0324d685d55ae45a95da0fc2d Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Mon, 8 Aug 2016 16:54:22 +0100 Subject: [PATCH] Added unit tests on the representation transformations Status: Passing all tests --- lib/qcd/QCD.h | 2 +- lib/qcd/action/pseudofermion/TwoFlavour.h | 4 +- lib/qcd/hmc/HmcRunner.h | 4 +- lib/qcd/hmc/integrators/Integrator.h | 5 +- lib/qcd/utils/SUn.h | 14 +++ lib/qcd/utils/SUnAdjoint.h | 23 ++++- tests/Test_hmc_WilsonAdjointFermionGauge.cc | 5 +- tests/Test_lie_generators.cc | 96 ++++++++++++++++++++- tests/Test_wilson_force.cc | 4 +- 9 files changed, 143 insertions(+), 14 deletions(-) diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index ddb0e217..1503bdf4 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -45,7 +45,7 @@ namespace QCD { static const int Zm = 6; static const int Tm = 7; - static const int Nc=3; + static const int Nc=2; static const int Ns=4; static const int Nd=4; static const int Nhs=2; // half spinor diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index d7acc8e7..6eaec166 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -130,8 +130,8 @@ class TwoFlavourPseudoFermionAction : public Action { MdagMLinearOperator, FermionField> MdagMOp(FermOp); X = zero; - DerivativeSolver(MdagMOp, Phi, X); - MdagMOp.Op(X, Y); + DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi + MdagMOp.Op(X, Y); // Y = M X = (Mdag)^-1 phi // Check hermiticity diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h index 1245b6ff..dafd2a01 100644 --- a/lib/qcd/hmc/HmcRunner.h +++ b/lib/qcd/hmc/HmcRunner.h @@ -103,7 +103,7 @@ class NerscHmcRunnerTemplate { // Create integrator, including the 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 @@ -111,7 +111,7 @@ 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 > diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index 8ae31574..4a9b2c7a 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -120,6 +120,7 @@ class Integrator { FieldType forceR(U._grid); // Implement smearing only for the fundamental representation now repr_set.at(a)->deriv(Rep.U, forceR); + forceR -= adj(forceR); GF force = Rep.RtoFundamentalProject(forceR); // Ta for the fundamental rep std::cout << GridLogIntegrator << "Hirep Force average: " @@ -200,10 +201,12 @@ class Integrator { template void operator()(std::vector*> repr_set, Repr& Rep, GridParallelRNG& pRNG) { - for (int a = 0; a < repr_set.size(); ++a) + 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{}; // Initialization of momenta and actions diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 26c7f32d..53bbb5f6 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -672,6 +672,20 @@ class SU { } } + // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) + // inverse operation: FundamentalLieAlgebraMatrix + static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { + conformable(h_out, in); + h_out = zero; + Matrix Ta; + + for (int a = 0; a < AdjointDimension; a++) { + generator(a, Ta); + auto tmp = - 2.0 * (trace(timesI(Ta) * in)) * scale;// 2.0 for the normalization of the trace in the fundamental rep + pokeColour(h_out, tmp, a); + } + std::cout << "h_out " << h_out << std::endl; + } template diff --git a/lib/qcd/utils/SUnAdjoint.h b/lib/qcd/utils/SUnAdjoint.h index 70f22f4c..6b206898 100644 --- a/lib/qcd/utils/SUnAdjoint.h +++ b/lib/qcd/utils/SUnAdjoint.h @@ -110,6 +110,22 @@ class SU_Adjoint : public SU { std::cout << GridLogMessage << std::endl; } + static void AdjointLieAlgebraMatrix( + const typename SU::LatticeAlgebraVector &h, + LatticeAdjMatrix &out, Real scale = 1.0) { + conformable(h, out); + GridBase *grid = out._grid; + LatticeAdjMatrix la(grid); + AMatrix iTa; + + out = zero; + for (int a = 0; a < Dimension; a++) { + generator(a, iTa); + la = peekColour(h, a) * iTa * scale; + out += la; + } + } + // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) static void projectOnAlgebra(typename SU::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) { conformable(h_out, in); @@ -118,9 +134,10 @@ class SU_Adjoint : public SU { for (int a = 0; a < Dimension; a++) { generator(a, iTa); - auto tmp = - 2.0 * real(trace(iTa * in)) * scale; + auto tmp = - 1.0/(ncolour) * (trace(iTa * in)) * scale;// 1/Nc for the normalization of the trace in the adj rep pokeColour(h_out, tmp, a); } + //std::cout << "h_out " << h_out << std::endl; } // a projector that keeps the generators stored to avoid the overhead of recomputing. @@ -130,12 +147,12 @@ class SU_Adjoint : public SU { h_out = zero; static bool precalculated = false; if (!precalculated){ - precalculated = true; + precalculated = true; for (int a = 0; a < Dimension; a++) generator(a, iTa[a]); } for (int a = 0; a < Dimension; a++) { - auto tmp = - 2.0*real(trace(iTa[a] * in)) * scale; + auto tmp = - 1.0/(ncolour) * (trace(iTa[a] * in)) * scale; pokeColour(h_out, tmp, a); } } diff --git a/tests/Test_hmc_WilsonAdjointFermionGauge.cc b/tests/Test_hmc_WilsonAdjointFermionGauge.cc index 935667a8..9c43587d 100644 --- a/tests/Test_hmc_WilsonAdjointFermionGauge.cc +++ b/tests/Test_hmc_WilsonAdjointFermionGauge.cc @@ -65,12 +65,13 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > { AdjointRepresentation::LatticeField U(UGrid); // Gauge action - WilsonGaugeActionR Waction(5.6); + WilsonGaugeActionR Waction(2.0); - Real mass = -0.77; + Real mass = -1.0; FermionAction FermOp(U, *FGrid, *FrbGrid, mass); ConjugateGradient CG(1.0e-8, 10000); + ConjugateResidual CR(1.0e-8, 10000); TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); diff --git a/tests/Test_lie_generators.cc b/tests/Test_lie_generators.cc index b30be869..533ccf9a 100644 --- a/tests/Test_lie_generators.cc +++ b/tests/Test_lie_generators.cc @@ -109,10 +109,104 @@ int main(int argc, char** argv) { // Testing HMC representation classes - AdjointRep<3> AdjRep(grid); + AdjointRep AdjRep(grid); // AdjointRepresentation has the predefined number of colours Nc Representations RepresentationTypes(grid); + LatticeGaugeField U(grid), V(grid); + SU::HotConfiguration(gridRNG, U); + SU::HotConfiguration(gridRNG, V); + + // Test group structure + // (U_f * V_f)_r = U_r * V_r + LatticeGaugeField UV(grid); + for (int mu = 0; mu < Nd; mu++) { + SU::LatticeMatrix Umu = peekLorentz(U,mu); + SU::LatticeMatrix Vmu = peekLorentz(V,mu); + pokeLorentz(UV,Umu*Vmu, mu); + } + + AdjRep.update_representation(UV); + typename AdjointRep::LatticeField UVr = AdjRep.U; // (U_f * V_f)_r + + + AdjRep.update_representation(U); + typename AdjointRep::LatticeField Ur = AdjRep.U; // U_r + + AdjRep.update_representation(V); + typename AdjointRep::LatticeField Vr = AdjRep.U; // V_r + + typename AdjointRep::LatticeField UrVr(grid); + for (int mu = 0; mu < Nd; mu++) { + typename AdjointRep::LatticeMatrix Urmu = peekLorentz(Ur,mu); + typename AdjointRep::LatticeMatrix Vrmu = peekLorentz(Vr,mu); + pokeLorentz(UrVr,Urmu*Vrmu, mu); + } + + typename AdjointRep::LatticeField Diff_check = UVr - UrVr; + std::cout << GridLogMessage << "Group structure check difference : " << norm2(Diff_check) << std::endl; + + // Check correspondence of algebra and group transformations + // Create a random vector + SU::LatticeAlgebraVector h_adj(grid); + typename AdjointRep::LatticeMatrix Ar(grid); + random(gridRNG,h_adj); + h_adj = real(h_adj); + SU_Adjoint::AdjointLieAlgebraMatrix(h_adj,Ar); + + // Re-extract h_adj + SU::LatticeAlgebraVector h_adj2(grid); + SU_Adjoint::projectOnAlgebra(h_adj2, Ar); + SU::LatticeAlgebraVector h_diff = h_adj - h_adj2; + std::cout << GridLogMessage << "Projections structure check vector difference : " << norm2(h_diff) << std::endl; + + // Exponentiate + typename AdjointRep::LatticeMatrix Uadj(grid); + Uadj = expMat(Ar, 1.0, 16); + + typename AdjointRep::LatticeMatrix uno(grid); + uno = 1.0; + // Check matrix Uadj, must be real orthogonal + typename AdjointRep::LatticeMatrix Ucheck = Uadj - conjugate(Uadj); + std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) + << std::endl; + + Ucheck = Uadj * adj(Uadj) - uno; + std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck) + << std::endl; + Ucheck = adj(Uadj) * Uadj - uno; + std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck) + << std::endl; + + + // Construct the fundamental matrix in the group + SU::LatticeMatrix Af(grid); + SU::FundamentalLieAlgebraMatrix(h_adj,Af); + SU::LatticeMatrix Ufund(grid); + Ufund = expMat(Af, 1.0, 16); + // Check unitarity + SU::LatticeMatrix uno_f(grid); + uno_f = 1.0; + SU::LatticeMatrix UnitCheck(grid); + UnitCheck = Ufund * adj(Ufund) - uno_f; + std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) + << std::endl; + UnitCheck = adj(Ufund) * Ufund - uno_f; + std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck) + << std::endl; + + + // Tranform to the adjoint representation + U = zero; // fill this with only one direction + pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields + + AdjRep.update_representation(U); + Ur = AdjRep.U; // U_r + typename AdjointRep::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj + + typename AdjointRep::LatticeMatrix Diff_check_mat = Ur0 - Uadj; + std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl; + Grid_finalize(); } diff --git a/tests/Test_wilson_force.cc b/tests/Test_wilson_force.cc index 0432fbfa..a6117897 100644 --- a/tests/Test_wilson_force.cc +++ b/tests/Test_wilson_force.cc @@ -58,7 +58,7 @@ int main (int argc, char ** argv) LatticeGaugeField U(&Grid); - SU3::HotConfiguration(pRNG,U); + SU2::HotConfiguration(pRNG,U); // SU3::ColdConfiguration(pRNG,U); //////////////////////////////////// @@ -95,7 +95,7 @@ int main (int argc, char ** argv) for(int mu=0;mu