portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Debugging process for the clover term

Browse Source
This commit is contained in:
Guido Cossu 2017-10-23 18:27:34 +01:00
parent 6391b2a1d0
commit 031c94e02e
3 changed files with 65 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
lib/qcd/action/fermion/WilsonCloverFermion.cc
View File

@ -71,10 +71,10 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
GridBase *grid = _Umu._grid; GridBase *grid = _Umu._grid;
typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
// Compute the field strength terms // Compute the field strength terms mu>nu
WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Ydir, Zdir); WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir); WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Xdir, Ydir); WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir); WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir); WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir); WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
@ -86,7 +86,7 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
CloverTerm += fillCloverXT(Ex); CloverTerm += fillCloverXT(Ex);
CloverTerm += fillCloverYT(Ey); CloverTerm += fillCloverYT(Ey);
CloverTerm += fillCloverZT(Ez); CloverTerm += fillCloverZT(Ez);
CloverTerm *= 0.5 * csw; // FieldStrength normalization? should be ( -i/8 ). Is it the anti-symmetric combination? CloverTerm *= (0.5) * csw;
int lvol = _Umu._grid->lSites(); int lvol = _Umu._grid->lSites();
@ -232,7 +232,8 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
out = *Clover * in; out = *Clover * in;
} else { } else {
Clover = (inv) ? &CloverTermInv : &CloverTerm; Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = adj(*Clover) * in; //out = adj(*Clover) * in;
out = adj(CloverTerm) * in;
} }

24
lib/qcd/action/fermion/WilsonCloverFermion.h
View File

@ -156,10 +156,10 @@ private:
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++) for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{ {
T._odata[i]()(0, 1) = timesI(F._odata[i]()()); //fixed T._odata[i]()(0, 1) = timesI(F._odata[i]()());
T._odata[i]()(1, 0) = timesI(F._odata[i]()()); //fixed T._odata[i]()(1, 0) = timesI(F._odata[i]()());
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); //fixed T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); //fixed T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
} }
return T; return T;
@ -172,10 +172,10 @@ private:
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++) for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{ {
T._odata[i]()(0, 1) = -(F._odata[i]()()); //fixed T._odata[i]()(0, 1) = -(F._odata[i]()());
T._odata[i]()(1, 0) = (F._odata[i]()()); //fixed T._odata[i]()(1, 0) = (F._odata[i]()());
T._odata[i]()(2, 3) = (F._odata[i]()()); //fixed T._odata[i]()(2, 3) = (F._odata[i]()());
T._odata[i]()(3, 2) = -(F._odata[i]()()); //fixed T._odata[i]()(3, 2) = -(F._odata[i]()());
} }
return T; return T;
@ -188,10 +188,10 @@ private:
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++) for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{ {
T._odata[i]()(0, 0) = timesI(F._odata[i]()()); //fixed T._odata[i]()(0, 0) = timesI(F._odata[i]()());
T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); //fixed T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); //fixed T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 3) = timesI(F._odata[i]()()); //fixed T._odata[i]()(3, 3) = timesI(F._odata[i]()());
} }
return T; return T;

52
tests/qdpxx/Test_qdpxx_wilson.cc
View File

@ -29,7 +29,7 @@
#include <Grid/Grid.h> #include <Grid/Grid.h>
// Mass // Mass
double mq = 0.1; double mq = 0.0;
// Define Wilson Types // Define Wilson Types
typedef Grid::QCD::WilsonImplR::FermionField FermionField; typedef Grid::QCD::WilsonImplR::FermionField FermionField;
@ -274,7 +274,7 @@ public:
p.Mass = _mq; p.Mass = _mq;
p.clovCoeffR = QDP::Real(1.0); p.clovCoeffR = QDP::Real(1.0);
p.clovCoeffT = QDP::Real(1.0); p.clovCoeffT = QDP::Real(1.0);
Real u0 = QDP::Real(0.0); Real u0 = QDP::Real(1.0);
Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs)); Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
@ -316,6 +316,8 @@ int main(int argc, char **argv)
FermionField src(UGrid); FermionField src(UGrid);
FermionField res_chroma(UGrid); FermionField res_chroma(UGrid);
FermionField res_grid(UGrid); FermionField res_grid(UGrid);
FermionField only_wilson(UGrid);
FermionField difference(UGrid);
std::vector<ChromaAction> ActionList({Wilson, WilsonClover}); std::vector<ChromaAction> ActionList({Wilson, WilsonClover});
std::vector<std::string> ActionName({"Wilson", "WilsonClover"}); std::vector<std::string> ActionName({"Wilson", "WilsonClover"});
@ -346,8 +348,19 @@ int main(int argc, char **argv)
std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl; std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl;
res_chroma = res_chroma - res_grid; difference = res_chroma - res_grid;
std::cout << "Norm of difference " << Grid::norm2(res_chroma) << std::endl; std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl;
// Isolate Clover term
calc_grid(Wilson, Ug, src, only_wilson, dag);// Wilson term
res_grid -= only_wilson;
res_chroma -= only_wilson;
std::cout << "Chroma:" << res_chroma << std::endl;
std::cout << "Grid :" << res_grid << std::endl;
} }
} }
@ -416,7 +429,36 @@ void make_gauge(GaugeField &Umu, FermionField &src)
Grid::GridParallelRNG RNG4(UGrid); Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4, Umu); Grid::QCD::SU3::HotConfiguration(RNG4, Umu);
Grid::gaussian(RNG4, src);
// Fermion field
//Grid::gaussian(RNG4, src);
Grid::QCD::SpinColourVector F;
Grid::Complex c;
std::vector<int> x(4); // 4d fermions
std::vector<int> gd = src._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
for (int sp = 0; sp < 1; sp++)
{
for (int j = 1; j < 2; j++)// colours
{
c = Grid::Complex(1.0, 0.0);
F()(sp)(j) = c;
}
}
Grid::pokeSite(F, src, x);
}
}
}
}
} }
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag) void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)