From 6391b2a1d0ae8128dafa25349be754a198b978d3 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Mon, 23 Oct 2017 14:42:35 +0100 Subject: [PATCH] Added test for Wilson and Clover fermions --- .vscode/settings.json | 3 +- tests/qdpxx/Test_qdpxx_wilson.cc | 463 +++++++++++++++++++++++++++++++ 2 files changed, 465 insertions(+), 1 deletion(-) create mode 100644 tests/qdpxx/Test_qdpxx_wilson.cc diff --git a/.vscode/settings.json b/.vscode/settings.json index 3e49029b..f13d503b 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -45,6 +45,7 @@ "istream": "cpp", "ostream": "cpp", "sstream": "cpp", - "streambuf": "cpp" + "streambuf": "cpp", + "algorithm": "cpp" } } \ No newline at end of file diff --git a/tests/qdpxx/Test_qdpxx_wilson.cc b/tests/qdpxx/Test_qdpxx_wilson.cc new file mode 100644 index 00000000..70a8b1bf --- /dev/null +++ b/tests/qdpxx/Test_qdpxx_wilson.cc @@ -0,0 +1,463 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc + + Copyright (C) 2017 + + Author: Guido Cossu + + 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 + +// Mass +double mq = 0.1; + +// Define Wilson Types +typedef Grid::QCD::WilsonImplR::FermionField FermionField; +typedef Grid::QCD::LatticeGaugeField GaugeField; + +#include +#include +#include + +enum ChromaAction +{ + Wilson, // Wilson + WilsonClover // CloverFermions +}; + +void make_gauge(GaugeField &lat, FermionField &src); +void calc_grid(ChromaAction CA, GaugeField &lat, FermionField &src, FermionField &res, int dag); +void calc_chroma(ChromaAction CA, GaugeField &lat, FermionField &src, FermionField &res, int dag); + +namespace Chroma +{ + +class ChromaWrapper +{ +public: + typedef multi1d U; + typedef LatticeFermion T4; + + static void ImportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd = gr._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]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + Grid::peekSite(LCM, gr, x); + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cc = LCM(mu)()(i, j); + c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc))); + QDP::pokeColor(cm, c, i, j); + } + } + QDP::pokeSite(ch[mu], cm, cx); + } + } + } + } + } + } + + static void ExportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::QCD::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + std::vector gd = gr._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]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cm = QDP::peekSite(ch[mu], cx); + c = QDP::peekColor(cm, i, j); + cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c))); + LCM(mu) + ()(i, j) = cc; + } + } + } + Grid::pokeSite(LCM, gr, x); + } + } + } + } + } + + // Specific for Wilson Fermions + static void ImportFermion(Grid::QCD::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // explicit 4d fermions in Grid + QDP::multi1d cx(4); + std::vector gd = gr._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]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + Grid::peekSite(F, gr, x); + + for (int j = 0; j < 3; j++) + { + for (int sp = 0; sp < 4; sp++) + { + + c = F()(sp)(j); + + cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c))); + + QDP::pokeSpin(cS, cc, sp); + } + QDP::pokeColor(cF, cS, j); + } + QDP::pokeSite(ch, cF, cx); + } + } + } + } + } + + // Specific for 4d Wilson fermions + static void ExportFermion(Grid::QCD::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // 4d fermions + QDP::multi1d cx(4); + std::vector gd = gr._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]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + cF = QDP::peekSite(ch, cx); + for (int sp = 0; sp < 4; sp++) + { + for (int j = 0; j < 3; j++) + { + cS = QDP::peekColor(cF, j); + cc = QDP::peekSpin(cS, sp); + c = Grid::Complex(QDP::toDouble(QDP::real(cc)), + QDP::toDouble(QDP::imag(cc))); + F() + (sp)(j) = c; + } + } + Grid::pokeSite(F, gr, x); + } + } + } + } + } + + static Handle> GetLinOp(U &u, ChromaAction params) + { + QDP::Real _mq(mq); + QDP::multi1d bcs(QDP::Nd); + + // Boundary conditions + bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1; + + if (params == Wilson) + { + + Chroma::WilsonFermActParams p; + p.Mass = _mq; + + Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); + Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); + Chroma::UnprecWilsonFermAct S_f(cfs, p); + Chroma::Handle> ffs(S_f.createState(u)); + return S_f.linOp(ffs); + } + + if (params == WilsonClover) + { + Chroma::CloverFermActParams p; + p.Mass = _mq; + p.clovCoeffR = QDP::Real(1.0); + p.clovCoeffT = QDP::Real(1.0); + Real u0 = QDP::Real(0.0); + + + Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); + Chroma::Handle> cfs(new Chroma::CreateSimpleFermState(fbc)); + Chroma::UnprecCloverFermAct S_f(cfs, p); + Chroma::Handle> ffs(S_f.createState(u)); + return S_f.linOp(ffs); + } + } +}; +} // namespace Chroma + +int main(int argc, char **argv) +{ + + /******************************************************** + * Setup QDP + *********************************************************/ + Chroma::initialize(&argc, &argv); + Chroma::WilsonTypeFermActs4DEnv::registerAll(); + + /******************************************************** + * Setup Grid + *********************************************************/ + Grid::Grid_init(&argc, &argv); + Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), + Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()), + Grid::GridDefaultMpi()); + + std::vector gd = UGrid->GlobalDimensions(); + QDP::multi1d nrow(QDP::Nd); + for (int mu = 0; mu < 4; mu++) + nrow[mu] = gd[mu]; + + QDP::Layout::setLattSize(nrow); + QDP::Layout::create(); + + GaugeField Ug(UGrid); + FermionField src(UGrid); + FermionField res_chroma(UGrid); + FermionField res_grid(UGrid); + + std::vector ActionList({Wilson, WilsonClover}); + std::vector ActionName({"Wilson", "WilsonClover"}); + + { + + for (int i = 0; i < ActionList.size(); i++) + { + std::cout << "*****************************" << std::endl; + std::cout << "Action " << ActionName[i] << std::endl; + std::cout << "*****************************" << std::endl; + make_gauge(Ug, src); // fills the gauge field and the fermion field with random numbers + + for (int dag = 0; dag < 2; dag++) + { + + { + + std::cout << "Dag = " << dag << std::endl; + + calc_chroma(ActionList[i], Ug, src, res_chroma, dag); + + // Remove the normalisation of Chroma Gauge links ???????? + std::cout << "Norm of Chroma " << ActionName[i] << " multiply " << Grid::norm2(res_chroma) << std::endl; + calc_grid(ActionList[i], Ug, src, res_grid, dag); + + std::cout << "Norm of gauge " << Grid::norm2(Ug) << std::endl; + + std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl; + + res_chroma = res_chroma - res_grid; + std::cout << "Norm of difference " << Grid::norm2(res_chroma) << std::endl; + } + } + + std::cout << "Finished test " << std::endl; + + Chroma::finalize(); + } + } +} + + +void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag) +{ + QDP::multi1d u(4); + Chroma::ChromaWrapper::ImportGauge(lat, u); + + QDP::LatticeFermion check; + QDP::LatticeFermion result; + QDP::LatticeFermion psi; + + Chroma::ChromaWrapper::ImportFermion(src, psi); + + for (int mu = 0; mu < 4; mu++) + { + std::cout << "Imported Gauge norm [" << mu << "] " << QDP::norm2(u[mu]) << std::endl; + } + std::cout << "Imported Fermion norm " << QDP::norm2(psi) << std::endl; + + typedef QDP::LatticeFermion T; + typedef QDP::multi1d U; + + auto linop = Chroma::ChromaWrapper::GetLinOp(u, action); + + printf("Calling Chroma Linop\n"); + fflush(stdout); + + if (dag) + (*linop)(check, psi, Chroma::MINUS); + else + (*linop)(check, psi, Chroma::PLUS); + + printf("Called Chroma Linop\n"); + fflush(stdout); + + // std::cout << "Calling Chroma Linop " << std::endl; + // linop->evenEvenLinOp(tmp, psi, isign); + // check[rb[0]] = tmp; + // linop->oddOddLinOp(tmp, psi, isign); + // check[rb[1]] = tmp; + // linop->evenOddLinOp(tmp, psi, isign); + // check[rb[0]] += tmp; + // linop->oddEvenLinOp(tmp, psi, isign); + // check[rb[1]] += tmp; + + Chroma::ChromaWrapper::ExportFermion(res, check); +} + +void make_gauge(GaugeField &Umu, FermionField &src) +{ + using namespace Grid; + using namespace Grid::QCD; + + std::vector seeds4({1, 2, 3, 4}); + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; + Grid::GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + Grid::QCD::SU3::HotConfiguration(RNG4, Umu); + Grid::gaussian(RNG4, src); +} + +void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; + Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + Grid::RealD _mass = mq; + + if (action == Wilson) + { + + Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass); + + std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl; + + if (dag) + Wf.Mdag(src, res); + else + Wf.M(src, res); + return; + } + + if (action == WilsonClover) + { + Grid::RealD _csw = 1.0; + + Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw); + Wf.ImportGauge(Umu); + + std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl; + + if (dag) + Wf.Mdag(src, res); + else + Wf.M(src, res); + return; + } + + assert(0); +}