From 91b8bf061323e5aeb73f84c41345a07447a90e47 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 26 Oct 2017 18:23:55 +0100 Subject: [PATCH] Debugging force term --- lib/qcd/action/fermion/FermionOperatorImpl.h | 52 ++++- lib/qcd/action/fermion/WilsonCloverFermion.cc | 30 +-- lib/qcd/action/fermion/WilsonCloverFermion.h | 128 ++++++++++- tests/forces/Test_wilson_force.cc | 7 +- tests/forces/Test_wilsonclover_force.cc | 209 ++++++++++++++++++ 5 files changed, 395 insertions(+), 31 deletions(-) create mode 100644 tests/forces/Test_wilsonclover_force.cc diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 9d24deb2..89bd9a15 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -254,8 +254,22 @@ namespace QCD { GaugeLinkField link(mat._grid); link = TraceIndex(outerProduct(Btilde,A)); PokeIndex(mat,link,mu); - } + } + + inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){ + mat = outerProduct(B,A); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + mat = TraceIndex(P); + } + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + for (int mu = 0; mu < Nd; mu++) + mat[mu] = PeekIndex(Uds, mu); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ int Ls=Btilde._grid->_fdimensions[0]; @@ -373,6 +387,19 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres assert(0); } + inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){ + assert(0); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + assert(0); + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + assert(0); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã, int mu) { assert(0); @@ -611,6 +638,25 @@ class GparityWilsonImpl : public ConjugateGaugeImpl(P); + parallel_for(auto ss = tmp.begin(); ss < tmp.end(); ss++) { + mat[ss]() = tmp[ss](0, 0) + conjugate(tmp[ss](1, 1)); + } + */ + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + assert(0); + } + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã, int mu) { int Ls = Btilde._grid->_fdimensions[0]; @@ -751,8 +797,8 @@ class StaggeredImpl : public PeriodicGaugeImpl(outerProduct(Btilde,A)); PokeIndex(mat,link,mu); - } - + } + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ assert (0); // Must never hit diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index 1dd12f52..2159fffc 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -43,11 +43,15 @@ RealD WilsonCloverFermion::M(const FermionField &in, FermionField &out) // Wilson term out.checkerboard = in.checkerboard; - this->Dhop(in, out, DaggerNo); + //this->Dhop(in, out, DaggerNo); // Clover term Mooee(in, temp); + //hack + out = zero; + + out += temp; return norm2(out); } @@ -59,11 +63,14 @@ RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) // Wilson term out.checkerboard = in.checkerboard; - this->Dhop(in, out, DaggerYes); + //this->Dhop(in, out, DaggerYes); // Clover term MooeeDag(in, temp); + //hack + out = zero; + out += temp; return norm2(out); } @@ -84,7 +91,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); // Compute the Clover Operator acting on Colour and Spin - CloverTerm = fillCloverYZ(Bx); + CloverTerm = fillCloverYZ(Bx); CloverTerm += fillCloverXZ(By); CloverTerm += fillCloverXY(Bz); CloverTerm += fillCloverXT(Ex); @@ -223,23 +230,6 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie } // MooeeInternal -// Derivative parts -template -void WilsonCloverFermion::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) -{ - - GaugeField tmp(mat._grid); - - conformable(U._grid, V._grid); - conformable(U._grid, mat._grid); - - mat.checkerboard = U.checkerboard; - tmp.checkerboard = U.checkerboard; - - this->DhopDeriv(mat, U, V, dag); - MooDeriv(tmp, U, V, dag); - mat += tmp; -} // Derivative parts template diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index 34482941..d8a42129 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -67,15 +67,18 @@ public: CloverTermOdd(&Hgrid), CloverTermInvEven(&Hgrid), CloverTermInvOdd(&Hgrid), - CloverTermDagEven(&Hgrid), - CloverTermDagOdd(&Hgrid), - CloverTermInvDagEven(&Hgrid), - CloverTermInvDagOdd(&Hgrid) + CloverTermDagEven(&Hgrid), + CloverTermDagOdd(&Hgrid), + CloverTermInvDagEven(&Hgrid), + CloverTermInvDagOdd(&Hgrid) { csw = _csw; assert(Nd == 4); // require 4 dimensions - if (csw == 0) std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl; + if (csw == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl; + + ImportGauge(_Umu); } virtual RealD M(const FermionField &in, FermionField &out); @@ -87,16 +90,127 @@ public: virtual void MooeeInvDag(const FermionField &in, FermionField &out); virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv); - virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); + //virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); virtual void MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); virtual void MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag); void ImportGauge(const GaugeField &_Umu); + // Derivative parts unpreconditioned pseudofermions + void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) + { + conformable(X._grid, Y._grid); + conformable(X._grid, force._grid); + GaugeLinkField force_mu(force._grid), lambda(force._grid); + GaugeField clover_force(force._grid); + PropagatorField Lambda(force._grid); + + // Here we are hitting some performance issues: + // need to extract the components of the DoubledGaugeField + // for each call + // Possible solution + // Create a vector object to store them? (cons: wasting space) + std::vector U(Nd, this->Umu._grid); + + Impl::extractLinkField(U, this->Umu); + + force = zero; + // Derivative of the Wilson hopping term + //this->DhopDeriv(force, X, Y, dag); + + /////////////////////////////////////////////////////////// + // Clover term derivative + /////////////////////////////////////////////////////////// + Impl::outerProductImpl(Lambda, X, Y); + + Gamma::Algebra sigma[] = { + Gamma::Algebra::SigmaXY, + Gamma::Algebra::SigmaXZ, + Gamma::Algebra::SigmaXT, + Gamma::Algebra::MinusSigmaXY, + Gamma::Algebra::SigmaYZ, + Gamma::Algebra::SigmaYT, + Gamma::Algebra::MinusSigmaXZ, + Gamma::Algebra::MinusSigmaYZ, + Gamma::Algebra::SigmaZT, + Gamma::Algebra::MinusSigmaXT, + Gamma::Algebra::MinusSigmaYT, + Gamma::Algebra::MinusSigmaZT}; + + /* + sigma_{\mu \nu}= + | 0 sigma[0] sigma[1] sigma[2] | + | sigma[3] 0 sigma[4] sigma[5] | + | sigma[6] sigma[7] 0 sigma[8] | + | sigma[9] sigma[10] sigma[11] 0 | + */ + + int count = 0; + clover_force = zero; + for (int mu = 0; mu < 4; mu++) + { + force_mu = zero; + for (int nu = 0; nu < 4; nu++) + { + if (mu == nu) continue; + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; + Impl::TraceSpinImpl(lambda, Slambda); //traceSpin + force_mu += Cmunu(U, lambda, mu, nu); + count++; + } + + pokeLorentz(clover_force, U[mu] * force_mu, mu); + } + clover_force *= csw / 8.; + force += clover_force; + } + + // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis + GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) + { + conformable(lambda._grid, U[0]._grid); + GaugeLinkField out(lambda._grid), tmp(lambda._grid); + + // insertion in upper staple + // please check redundancy of shift operations + + // C1+ + tmp = lambda * U[nu]; + out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C2+ + tmp = U[mu] * Impl::CovShiftIdentityForward(adj(lambda), mu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C3+ + tmp = U[nu] * Impl::CovShiftIdentityForward(adj(lambda), nu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); + + // C4+ + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda; + + // insertion in lower staple + // C1- + out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C2- + tmp = adj(lambda) * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C3- + tmp = lambda * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); + + // C4- + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda; + + return out; + } + private: // here fixing the 4 dimensions, make it more general? - RealD csw; // Clover coefficient + RealD csw; // Clover coefficient CloverFieldType CloverTerm, CloverTermInv; // Clover term CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO diff --git a/tests/forces/Test_wilson_force.cc b/tests/forces/Test_wilson_force.cc index 1f34a48a..f834726b 100644 --- a/tests/forces/Test_wilson_force.cc +++ b/tests/forces/Test_wilson_force.cc @@ -50,7 +50,12 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); - pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); + std::vector vrand(4); + std::srand(std::time(0)); + std::generate(vrand.begin(), vrand.end(), std::rand); + std::cout << GridLogMessage << vrand << std::endl; + pRNG.SeedFixedIntegers(vrand); + //pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); LatticeFermion phi (&Grid); gaussian(pRNG,phi); LatticeFermion Mphi (&Grid); diff --git a/tests/forces/Test_wilsonclover_force.cc b/tests/forces/Test_wilsonclover_force.cc new file mode 100644 index 00000000..c99cfa98 --- /dev/null +++ b/tests/forces/Test_wilsonclover_force.cc @@ -0,0 +1,209 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_force.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size, simd_layout, mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + std::vector seeds({1, 2, 3, 4}); + + GridParallelRNG pRNG(&Grid); + std::vector vrand(4); + std::srand(std::time(0)); + std::generate(vrand.begin(), vrand.end(), std::rand); + std::cout << GridLogMessage << vrand << std::endl; + pRNG.SeedFixedIntegers(vrand); + + LatticeFermion phi(&Grid); + gaussian(pRNG, phi); + LatticeFermion Mphi(&Grid); + LatticeFermion MphiPrime(&Grid); + + LatticeGaugeField U(&Grid); + + SU3::HotConfiguration(pRNG,U); + //SU3::ColdConfiguration(pRNG, U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass = -4.0; //kills the diagonal term + Real csw = 1.0; + WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw); + + Dw.M(phi, Mphi); + ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(&Grid); + LatticeGaugeField tmp(&Grid); + + Dw.MDeriv(tmp, Mphi, phi, DaggerNo); UdSdU = tmp; + Dw.MDeriv(tmp, phi, Mphi, DaggerYes); UdSdU += tmp; + // Take the traceless antihermitian component + UdSdU = Ta(UdSdU); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.0001; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; + RealD Hmompp = 0.0; + LatticeColourMatrix mommu(&Grid); + LatticeColourMatrix forcemu(&Grid); + LatticeGaugeField mom(&Grid); + LatticeGaugeField Uprime(&Grid); + + for (int mu = 0; mu < Nd; mu++) { + // Traceless antihermitian momentum; gaussian in lie alg + SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); + Hmom -= real(sum(trace(mommu * mommu))); + PokeIndex(mom, mommu, mu); + } + + parallel_for(int ss=0;ssoSites();ss++){ + for (int mu = 0; mu < Nd; mu++) + Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]); + } + + std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl; + + // New action + Dw.ImportGauge(Uprime); + Dw.M(phi, MphiPrime); + + ComplexD Sprime = innerProduct(MphiPrime, MphiPrime); + + ////////////////////////////////////////////// + // Use derivative to estimate dS + ////////////////////////////////////////////// + + /////////////////////////////////////////////////////// + std::cout << GridLogMessage << "Antihermiticity tests - 1 " << std::endl; + for (int mu = 0; mu < Nd; mu++) + { + mommu = PeekIndex(mom, mu); + std::cout << GridLogMessage << " Mommu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " Test: Mommu + Mommudag " << norm2(mommu) << std::endl; + mommu = PeekIndex(UdSdU, mu); + std::cout << GridLogMessage << " dsdumu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " Test: dsdumu + dag " << norm2(mommu) << std::endl; + std::cout << "" << std::endl; + } + //////////////////////////////////////////////////////// + + LatticeComplex dS(&Grid); + dS = zero; + LatticeComplex dSmom(&Grid); + dSmom = zero; + LatticeComplex dSmom2(&Grid); + dSmom2 = zero; + + + // need for this??? + // ultimately it is just a 2.0 factor in UdSdU + for (int mu = 0; mu < Nd; mu++) + { + mommu = PeekIndex(UdSdU, mu); // P_mu = + mommu = Ta(mommu) * 2.0; // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag) + PokeIndex(UdSdU, mommu, mu); // UdSdU_mu = Mom + } + + std::cout << GridLogMessage<< "Antihermiticity tests - 2 " << std::endl; + for (int mu = 0; mu < Nd; mu++) + { + mommu = PeekIndex(mom, mu); + std::cout << GridLogMessage << " Mommu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " Mommu + Mommudag " << norm2(mommu) << std::endl; + mommu = PeekIndex(UdSdU, mu); + std::cout << GridLogMessage << " dsdumu " << norm2(mommu) << std::endl; + mommu = mommu + adj(mommu); + std::cout << GridLogMessage << " dsdumu + dag " << norm2(mommu) << std::endl; + std::cout << "" << std::endl; + } + ///////////////////////////////////////////////////// + + + for (int mu = 0; mu < Nd; mu++) + { + forcemu = PeekIndex(UdSdU, mu); + mommu = PeekIndex(mom, mu); + + // Update PF action density + dS = dS + trace(mommu * forcemu) * dt; + + dSmom = dSmom - trace(mommu * forcemu) * dt; + dSmom2 = dSmom2 - trace(forcemu * forcemu) * (0.25 * dt * dt); + + // Update mom action density + mommu = mommu + forcemu * (dt * 0.5); + + Hmomprime -= real(sum(trace(mommu * mommu))); + } + + ComplexD dSpred = sum(dS); + ComplexD dSm = sum(dSmom); + ComplexD dSm2 = sum(dSmom2); + + std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl; + std::cout << GridLogMessage << "Final mom hamiltonian is " << Hmomprime << std::endl; + std::cout << GridLogMessage << "Delta mom hamiltonian is " << Hmomprime - Hmom << std::endl; + + std::cout << GridLogMessage << " S " << S << std::endl; + std::cout << GridLogMessage << " Sprime " << Sprime << std::endl; + std::cout << GridLogMessage << "dS " << Sprime - S << std::endl; + std::cout << GridLogMessage << "predict dS " << dSpred << std::endl; + std::cout << GridLogMessage << "dSm " << dSm << std::endl; + std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl; + + std::cout << GridLogMessage << "Total dS " << Hmomprime - Hmom + Sprime - S << std::endl; + + assert(fabs(real(Sprime - S - dSpred)) < 1.0); + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +}