diff --git a/.gitignore b/.gitignore index 399f2f6b..dc59879f 100644 --- a/.gitignore +++ b/.gitignore @@ -123,3 +123,9 @@ make-bin-BUCK.sh ##################### lib/qcd/spin/gamma-gen/*.h lib/qcd/spin/gamma-gen/*.cc + +# vs code editor files # +######################## +.vscode/ +.vscode/settings.json +settings.json diff --git a/.travis.yml b/.travis.yml index 7d8203ce..ad4e5b73 100644 --- a/.travis.yml +++ b/.travis.yml @@ -44,3 +44,4 @@ script: - make -j4 - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals - make check + diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 6e123660..135934fb 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -47,6 +47,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include @@ -55,6 +56,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp new file mode 100644 index 00000000..c369f086 --- /dev/null +++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp @@ -0,0 +1,153 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MAction/Wilson.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Antonin Portelli + +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 */ + +#ifndef Hadrons_MAction_WilsonClover_hpp_ +#define Hadrons_MAction_WilsonClover_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TWilson quark action * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MAction) + +class WilsonCloverPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverPar, + std::string, gauge, + double , mass, + double , csw_r, + double , csw_t, + WilsonAnisotropyCoefficients ,clover_anisotropy, + std::string, boundary + ); +}; + +template +class TWilsonClover: public Module +{ +public: + FGS_TYPE_ALIASES(FImpl,); +public: + // constructor + TWilsonClover(const std::string name); + // destructor + virtual ~TWilsonClover(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(WilsonClover, TWilsonClover, MAction); + +/****************************************************************************** + * TWilsonClover template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWilsonClover::TWilsonClover(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWilsonClover::getInput(void) +{ + std::vector in = {par().gauge}; + + return in; +} + +template +std::vector TWilsonClover::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWilsonClover::setup(void) +{ + //unsigned int size; + + // size = 2*env().template lattice4dSize(); + // env().registerObject(getName(), size); + + + LOG(Message) << "Setting up TWilsonClover fermion matrix with m= " << par().mass + << " using gauge field '" << par().gauge << "'" << std::endl; + LOG(Message) << "Fermion boundary conditions: " << par().boundary + << std::endl; + LOG(Message) << "Clover term csw_r: " << par().csw_r + << " csw_t: " << par().csw_t + << std::endl; + auto &U = envGet(LatticeGaugeField, par().gauge); + auto &grid = *env().getGrid(); + auto &gridRb = *env().getRbGrid(); + std::vector boundary = strToVec(par().boundary); + typename WilsonCloverFermion::ImplParams implParams(boundary); + envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, + par().csw_r, + par().csw_t, + par().clover_anisotropy, + implParams); + + + //FMat *fMatPt = new WilsonCloverFermion(U, grid, gridRb, par().mass, + // par().csw_r, + // par().csw_t, + // par().clover_anisotropy, + // implParams); + //env().setObject(getName(), fMatPt); + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWilsonClover::execute() +{ + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WilsonClover_hpp_ diff --git a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp index 33787a0b..412e92d5 100644 --- a/extras/Hadrons/Modules/MFermion/GaugeProp.hpp +++ b/extras/Hadrons/Modules/MFermion/GaugeProp.hpp @@ -94,7 +94,6 @@ private: }; MODULE_REGISTER_NS(GaugeProp, TGaugeProp, MFermion); - /****************************************************************************** * TGaugeProp implementation * ******************************************************************************/ @@ -154,7 +153,7 @@ void TGaugeProp::execute(void) LOG(Message) << "Inverting using solver '" << par().solver << "' on source '" << par().source << "'" << std::endl; for (unsigned int s = 0; s < Ns; ++s) - for (unsigned int c = 0; c < Nc; ++c) + for (unsigned int c = 0; c < FImpl::Dimension; ++c) { LOG(Message) << "Inversion for spin= " << s << ", color= " << c << std::endl; @@ -163,11 +162,11 @@ void TGaugeProp::execute(void) { if (Ls_ == 1) { - PropToFerm(source, fullSrc, s, c); + PropToFerm(source, fullSrc, s, c); } else { - PropToFerm(tmp, fullSrc, s, c); + PropToFerm(tmp, fullSrc, s, c); make_5D(tmp, source, Ls_); } } @@ -180,18 +179,18 @@ void TGaugeProp::execute(void) } else { - PropToFerm(source, fullSrc, s, c); + PropToFerm(source, fullSrc, s, c); } } sol = zero; solver(sol, source); - FermToProp(prop, sol, s, c); + FermToProp(prop, sol, s, c); // create 4D propagators from 5D one if necessary if (Ls_ > 1) { PropagatorField &p4d = envGet(PropagatorField, getName()); make_4D(sol, tmp, Ls_); - FermToProp(p4d, tmp, s, c); + FermToProp(p4d, tmp, s, c); } } } diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.cc b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc new file mode 100644 index 00000000..f15a3b7c --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.cc @@ -0,0 +1,75 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc + +Copyright (C) 2015 +Copyright (C) 2016 + + +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 Grid; +using namespace Hadrons; +using namespace MGauge; + +// constructor ///////////////////////////////////////////////////////////////// +template +TFundtoHirep::TFundtoHirep(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFundtoHirep::getInput(void) +{ + std::vector in; + return in; +} + +template +std::vector TFundtoHirep::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TFundtoHirep::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFundtoHirep::execute(void) +{ + auto &U = *env().template getObject(par().gaugeconf); + LOG(Message) << "Transforming Representation" << std::endl; + + Rep TargetRepresentation(U._grid); + TargetRepresentation.update_representation(U); + + typename Rep::LatticeField &URep = *env().template createLattice(getName()); + URep = TargetRepresentation.U; +} diff --git a/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp new file mode 100644 index 00000000..6f072783 --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/FundtoHirep.hpp @@ -0,0 +1,77 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: David Preti + 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 */ + +#ifndef Hadrons_MGauge_FundtoHirep_hpp_ +#define Hadrons_MGauge_FundtoHirep_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Load a NERSC configuration * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class FundtoHirepPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FundtoHirepPar, + std::string, gaugeconf); +}; + +template +class TFundtoHirep: public Module +{ +public: + // constructor + TFundtoHirep(const std::string name); + // destructor + virtual ~TFundtoHirep(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + void setup(void); + // execution + void execute(void); +}; + +//MODULE_REGISTER_NS(FundtoAdjoint, TFundtoHirep, MGauge); +//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep, MGauge); +//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGauge_FundtoHirep_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 1c71301a..90602275 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,10 +1,13 @@ modules_cc =\ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc \ Modules/MContraction/WeakHamiltonianEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MGauge/Unit.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ + Modules/MGauge/FundtoHirep.cc \ Modules/MScalar/FreeProp.cc \ Modules/MScalar/ChargedProp.cc \ Modules/MIO/LoadNersc.cc @@ -31,6 +34,7 @@ modules_hpp =\ Modules/MGauge/Unit.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/StochEm.hpp \ + Modules/MGauge/FundtoHirep.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ @@ -39,6 +43,7 @@ modules_hpp =\ Modules/MScalar/ChargedProp.hpp \ Modules/MAction/DWF.hpp \ Modules/MAction/Wilson.hpp \ + Modules/MAction/WilsonClover.hpp \ Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/TwoPoint.hpp \ diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index 26746e6e..38444c04 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -183,11 +183,13 @@ namespace Grid { virtual RealD Mpc (const Field &in, Field &out) =0; virtual RealD MpcDag (const Field &in, Field &out) =0; virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { - Field tmp(in._grid); + Field tmp(in._grid); + tmp.checkerboard = in.checkerboard; ni=Mpc(in,tmp); no=MpcDag(tmp,out); } virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + out.checkerboard = in.checkerboard; MpcDagMpc(in,out,n1,n2); } virtual void HermOp(const Field &in, Field &out){ @@ -215,13 +217,15 @@ namespace Grid { public: SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){}; virtual RealD Mpc (const Field &in, Field &out) { - Field tmp(in._grid); -// std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; + Field tmp(in._grid); + tmp.checkerboard = !in.checkerboard; + //std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; _Mat.Meooe(in,tmp); _Mat.MooeeInv(tmp,out); _Mat.Meooe(out,tmp); + //std::cout << "cb in " << in.checkerboard << " cb out " << out.checkerboard << std::endl; _Mat.Mooee(in,out); return axpy_norm(out,-1.0,tmp,out); } diff --git a/lib/cartesian/Cartesian_base.h b/lib/cartesian/Cartesian_base.h index c49dac84..05a8a3da 100644 --- a/lib/cartesian/Cartesian_base.h +++ b/lib/cartesian/Cartesian_base.h @@ -79,6 +79,8 @@ public: std::vector _lstart; // local start of array in gcoors _processor_coor[d]*_ldimensions[d] std::vector _lend ; // local end of array in gcoors _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1 + bool _isCheckerBoarded; + public: //////////////////////////////////////////////////////////////// diff --git a/lib/cartesian/Cartesian_full.h b/lib/cartesian/Cartesian_full.h index b2372575..b6297d3d 100644 --- a/lib/cartesian/Cartesian_full.h +++ b/lib/cartesian/Cartesian_full.h @@ -97,6 +97,7 @@ public: /////////////////////// // Grid information /////////////////////// + _isCheckerBoarded = false; _ndimension = dimensions.size(); _fdimensions.resize(_ndimension); diff --git a/lib/cartesian/Cartesian_red_black.h b/lib/cartesian/Cartesian_red_black.h index ee424385..5a041f65 100644 --- a/lib/cartesian/Cartesian_red_black.h +++ b/lib/cartesian/Cartesian_red_black.h @@ -171,9 +171,8 @@ public: const std::vector &checker_dim_mask, int checker_dim) { - /////////////////////// - // Grid information - /////////////////////// + + _isCheckerBoarded = true; _checker_dim = checker_dim; assert(checker_dim_mask[checker_dim] == 1); _ndimension = dimensions.size(); diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 531b71bd..e1e38e82 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -39,6 +39,7 @@ namespace QCD { static const int Zdir = 2; static const int Tdir = 3; + static const int Xp = 0; static const int Yp = 1; static const int Zp = 2; @@ -420,15 +421,16 @@ namespace QCD { ////////////////////////////////////////////// // Fermion <-> propagator assignements ////////////////////////////////////////////// - template - void FermToProp(Prop &p, const Ferm &f, const int s, const int c) + //template + template + void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c) { - for(int j = 0; j < Ns; ++j) + for(int j = 0; j < Ns; ++j) { auto pjs = peekSpin(p, j, s); auto fj = peekSpin(f, j); - for(int i = 0; i < Nc; ++i) + for(int i = 0; i < Fimpl::Dimension; ++i) { pokeColour(pjs, peekColour(fj, i), i, c); } @@ -436,15 +438,16 @@ namespace QCD { } } - template - void PropToFerm(Ferm &f, const Prop &p, const int s, const int c) + //template + template + void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c) { for(int j = 0; j < Ns; ++j) { auto pjs = peekSpin(p, j, s); auto fj = peekSpin(f, j); - for(int i = 0; i < Nc; ++i) + for(int i = 0; i < Fimpl::Dimension; ++i) { pokeColour(fj, peekColour(pjs, i, c), i); } @@ -503,38 +506,6 @@ namespace QCD { } //namespace QCD } // Grid -/* -<<<<<<< HEAD -#include -#include -#include -#include -#include - -// Include representations -#include -#include -#include -#include - -// Scalar field -#include - -#include - -#include - -#include -#include -#include -#include - - -//#include -======= - ->>>>>>> develop -*/ #endif diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index ad2f383d..2a008cb7 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -50,11 +50,13 @@ Author: Peter Boyle //////////////////////////////////////////// #include // 4d wilson like -#include // 4d wilson like +#include // 4d wilson like +#include // 4d wilson clover fermions #include // 5d base used by all 5d overlap types -//#include + #include #include + #include // Cayley types #include #include @@ -104,10 +106,33 @@ typedef WilsonFermion WilsonTwoIndexSymmetricFermi typedef WilsonFermion WilsonTwoIndexSymmetricFermionF; typedef WilsonFermion WilsonTwoIndexSymmetricFermionD; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionR; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionF; +typedef WilsonFermion WilsonTwoIndexAntiSymmetricFermionD; + +// Twisted mass fermion typedef WilsonTMFermion WilsonTMFermionR; typedef WilsonTMFermion WilsonTMFermionF; typedef WilsonTMFermion WilsonTMFermionD; +// Clover fermions +typedef WilsonCloverFermion WilsonCloverFermionR; +typedef WilsonCloverFermion WilsonCloverFermionF; +typedef WilsonCloverFermion WilsonCloverFermionD; + +typedef WilsonCloverFermion WilsonCloverAdjFermionR; +typedef WilsonCloverFermion WilsonCloverAdjFermionF; +typedef WilsonCloverFermion WilsonCloverAdjFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionD; + +// Domain Wall fermions typedef DomainWallFermion DomainWallFermionR; typedef DomainWallFermion DomainWallFermionF; typedef DomainWallFermion DomainWallFermionD; diff --git a/lib/qcd/action/fermion/FermionCore.h b/lib/qcd/action/fermion/FermionCore.h index 17006961..60632c3a 100644 --- a/lib/qcd/action/fermion/FermionCore.h +++ b/lib/qcd/action/fermion/FermionCore.h @@ -70,7 +70,9 @@ Author: Peter Boyle #define TwoIndexFermOpTemplateInstantiate(A) \ template class A; \ - template class A; + template class A; \ + template class A; \ + template class A; #define FermOp5dVecTemplateInstantiate(A) \ template class A; \ diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index a9510222..c21a07ee 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -164,6 +164,7 @@ namespace QCD { public: static const int Dimension = Representation::Dimension; + static const bool isFundamental = Representation::isFundamental; static const bool LsVectorised=false; static const int Nhcs = Options::Nhcs; @@ -261,8 +262,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]; @@ -284,27 +299,28 @@ namespace QCD { //////////////////////////////////////////////////////////////////////////////////// // Single flavour four spinors with colour index, 5d redblack //////////////////////////////////////////////////////////////////////////////////// -template -class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > { +template +class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Representation::Dimension> > { public: - typedef PeriodicGaugeImpl > Gimpl; + typedef PeriodicGaugeImpl > Gimpl; INHERIT_GIMPL_TYPES(Gimpl); - static const int Dimension = Nrepresentation; + static const int Dimension = Representation::Dimension; + static const bool isFundamental = Representation::isFundamental; static const bool LsVectorised=true; static const int Nhcs = Options::Nhcs; typedef typename Options::_Coeff_t Coeff_t; typedef typename Options::template PrecisionMapper::LowerPrecVector SimdL; - template using iImplSpinor = iScalar, Ns> >; - template using iImplPropagator = iScalar, Ns> >; - template using iImplHalfSpinor = iScalar, Nhs> >; - template using iImplHalfCommSpinor = iScalar, Nhcs> >; - template using iImplDoubledGaugeField = iVector >, Nds>; - template using iImplGaugeField = iVector >, Nd>; - template using iImplGaugeLink = iScalar > >; + template using iImplSpinor = iScalar, Ns> >; + template using iImplPropagator = iScalar, Ns> >; + template using iImplHalfSpinor = iScalar, Nhs> >; + template using iImplHalfCommSpinor = iScalar, Nhcs> >; + template using iImplDoubledGaugeField = iVector >, Nds>; + template using iImplGaugeField = iVector >, Nd>; + template using iImplGaugeLink = iScalar > >; typedef iImplSpinor SiteSpinor; typedef iImplPropagator SitePropagator; @@ -340,8 +356,8 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres const SiteHalfSpinor &chi, int mu, StencilEntry *SE, StencilImpl &St) { SiteGaugeLink UU; - for (int i = 0; i < Nrepresentation; i++) { - for (int j = 0; j < Nrepresentation; j++) { + for (int i = 0; i < Dimension; i++) { + for (int j = 0; j < Dimension; j++) { vsplat(UU()()(i, j), U(mu)()(i, j)); } } @@ -353,8 +369,8 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres const SitePropagator &chi, int mu) { SiteGaugeLink UU; - for (int i = 0; i < Nrepresentation; i++) { - for (int j = 0; j < Nrepresentation; j++) { + for (int i = 0; i < Dimension; i++) { + for (int j = 0; j < Dimension; j++) { vsplat(UU()()(i, j), U(mu)()(i, j)); } } @@ -393,6 +409,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); @@ -445,25 +474,26 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres //////////////////////////////////////////////////////////////////////////////////////// // Flavour doubled spinors; is Gparity the only? what about C*? //////////////////////////////////////////////////////////////////////////////////////// -template -class GparityWilsonImpl : public ConjugateGaugeImpl > { +template +class GparityWilsonImpl : public ConjugateGaugeImpl > { public: - static const int Dimension = Nrepresentation; + static const int Dimension = Representation::Dimension; + static const bool isFundamental = Representation::isFundamental; static const int Nhcs = Options::Nhcs; static const bool LsVectorised=false; - typedef ConjugateGaugeImpl< GaugeImplTypes > Gimpl; + typedef ConjugateGaugeImpl< GaugeImplTypes > Gimpl; INHERIT_GIMPL_TYPES(Gimpl); typedef typename Options::_Coeff_t Coeff_t; typedef typename Options::template PrecisionMapper::LowerPrecVector SimdL; - template using iImplSpinor = iVector, Ns>, Ngp>; - template using iImplPropagator = iVector, Ns>, Ngp>; - template using iImplHalfSpinor = iVector, Nhs>, Ngp>; - template using iImplHalfCommSpinor = iVector, Nhcs>, Ngp>; - template using iImplDoubledGaugeField = iVector >, Nds>, Ngp>; + template using iImplSpinor = iVector, Ns>, Ngp>; + template using iImplPropagator = iVector, Ns>, Ngp>; + template using iImplHalfSpinor = iVector, Nhs>, Ngp>; + template using iImplHalfCommSpinor = iVector, Nhcs>, Ngp>; + template using iImplDoubledGaugeField = iVector >, Nds>, Ngp>; typedef iImplSpinor SiteSpinor; typedef iImplPropagator SitePropagator; @@ -636,6 +666,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]; @@ -665,6 +714,7 @@ class StaggeredImpl : public PeriodicGaugeImpl > Gimpl; @@ -776,8 +826,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 @@ -793,6 +843,7 @@ class StaggeredImpl : public PeriodicGaugeImpl > Gimpl; @@ -983,29 +1034,33 @@ typedef WilsonImpl Wilso typedef WilsonImpl WilsonTwoIndexSymmetricImplF; // Float typedef WilsonImpl WilsonTwoIndexSymmetricImplD; // Double -typedef DomainWallVec5dImpl DomainWallVec5dImplR; // Real.. whichever prec -typedef DomainWallVec5dImpl DomainWallVec5dImplF; // Float -typedef DomainWallVec5dImpl DomainWallVec5dImplD; // Double +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplR; // Real.. whichever prec +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplF; // Float +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplD; // Double + +typedef DomainWallVec5dImpl DomainWallVec5dImplR; // Real.. whichever prec +typedef DomainWallVec5dImpl DomainWallVec5dImplF; // Float +typedef DomainWallVec5dImpl DomainWallVec5dImplD; // Double -typedef DomainWallVec5dImpl DomainWallVec5dImplRL; // Real.. whichever prec -typedef DomainWallVec5dImpl DomainWallVec5dImplFH; // Float -typedef DomainWallVec5dImpl DomainWallVec5dImplDF; // Double +typedef DomainWallVec5dImpl DomainWallVec5dImplRL; // Real.. whichever prec +typedef DomainWallVec5dImpl DomainWallVec5dImplFH; // Float +typedef DomainWallVec5dImpl DomainWallVec5dImplDF; // Double -typedef DomainWallVec5dImpl ZDomainWallVec5dImplR; // Real.. whichever prec -typedef DomainWallVec5dImpl ZDomainWallVec5dImplF; // Float -typedef DomainWallVec5dImpl ZDomainWallVec5dImplD; // Double +typedef DomainWallVec5dImpl ZDomainWallVec5dImplR; // Real.. whichever prec +typedef DomainWallVec5dImpl ZDomainWallVec5dImplF; // Float +typedef DomainWallVec5dImpl ZDomainWallVec5dImplD; // Double -typedef DomainWallVec5dImpl ZDomainWallVec5dImplRL; // Real.. whichever prec -typedef DomainWallVec5dImpl ZDomainWallVec5dImplFH; // Float -typedef DomainWallVec5dImpl ZDomainWallVec5dImplDF; // Double +typedef DomainWallVec5dImpl ZDomainWallVec5dImplRL; // Real.. whichever prec +typedef DomainWallVec5dImpl ZDomainWallVec5dImplFH; // Float +typedef DomainWallVec5dImpl ZDomainWallVec5dImplDF; // Double -typedef GparityWilsonImpl GparityWilsonImplR; // Real.. whichever prec -typedef GparityWilsonImpl GparityWilsonImplF; // Float -typedef GparityWilsonImpl GparityWilsonImplD; // Double +typedef GparityWilsonImpl GparityWilsonImplR; // Real.. whichever prec +typedef GparityWilsonImpl GparityWilsonImplF; // Float +typedef GparityWilsonImpl GparityWilsonImplD; // Double -typedef GparityWilsonImpl GparityWilsonImplRL; // Real.. whichever prec -typedef GparityWilsonImpl GparityWilsonImplFH; // Float -typedef GparityWilsonImpl GparityWilsonImplDF; // Double +typedef GparityWilsonImpl GparityWilsonImplRL; // Real.. whichever prec +typedef GparityWilsonImpl GparityWilsonImplFH; // Float +typedef GparityWilsonImpl GparityWilsonImplDF; // Double typedef StaggeredImpl StaggeredImplR; // Real.. whichever prec typedef StaggeredImpl StaggeredImplF; // Float diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc new file mode 100644 index 00000000..3c082446 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -0,0 +1,243 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc + + Copyright (C) 2017 + + Author: paboyle + 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 +#include +#include + +namespace Grid +{ +namespace QCD +{ + +// *NOT* EO +template +RealD WilsonCloverFermion::M(const FermionField &in, FermionField &out) +{ + FermionField temp(out._grid); + + // Wilson term + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerNo); + + // Clover term + Mooee(in, temp); + + out += temp; + return norm2(out); +} + +template +RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) +{ + FermionField temp(out._grid); + + // Wilson term + out.checkerboard = in.checkerboard; + this->Dhop(in, out, DaggerYes); + + // Clover term + MooeeDag(in, temp); + + out += temp; + return norm2(out); +} + +template +void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) +{ + WilsonFermion::ImportGauge(_Umu); + GridBase *grid = _Umu._grid; + typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); + + // Compute the field strength terms mu>nu + WilsonLoops::FieldStrength(Bx, _Umu, Zdir, Ydir); + WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); + WilsonLoops::FieldStrength(Bz, _Umu, Ydir, Xdir); + WilsonLoops::FieldStrength(Ex, _Umu, Tdir, Xdir); + WilsonLoops::FieldStrength(Ey, _Umu, Tdir, Ydir); + WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); + + // Compute the Clover Operator acting on Colour and Spin + // multiply here by the clover coefficients for the anisotropy + CloverTerm = fillCloverYZ(Bx) * csw_r; + CloverTerm += fillCloverXZ(By) * csw_r; + CloverTerm += fillCloverXY(Bz) * csw_r; + CloverTerm += fillCloverXT(Ex) * csw_t; + CloverTerm += fillCloverYT(Ey) * csw_t; + CloverTerm += fillCloverZT(Ez) * csw_t; + CloverTerm += diag_mass; + + int lvol = _Umu._grid->lSites(); + int DimRep = Impl::Dimension; + + Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + + std::vector lcoor; + typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero; + + for (int site = 0; site < lvol; site++) + { + grid->LocalIndexToLocalCoor(site, lcoor); + EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + peekLocalSite(Qx, CloverTerm, lcoor); + Qxinv = zero; + //if (csw!=0){ + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b); + // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; + + EigenInvCloverOp = EigenCloverOp.inverse(); + //std::cout << EigenInvCloverOp << std::endl; + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); + // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; + // } + pokeLocalSite(Qxinv, CloverTermInv, lcoor); + } + + // Separate the even and odd parts + pickCheckerboard(Even, CloverTermEven, CloverTerm); + pickCheckerboard(Odd, CloverTermOdd, CloverTerm); + + pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); + pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm)); + + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); + pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv); + + pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); + pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); +} + +template +void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerNo, InverseNo); +} + +template +void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerYes, InverseNo); +} + +template +void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerNo, InverseYes); +} + +template +void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) +{ + this->MooeeInternal(in, out, DaggerYes, InverseYes); +} + +template +void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv) +{ + out.checkerboard = in.checkerboard; + CloverFieldType *Clover; + assert(in.checkerboard == Odd || in.checkerboard == Even); + + if (dag) + { + if (in._grid->_isCheckerBoarded) + { + if (in.checkerboard == Odd) + { + Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; + } + else + { + Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven; + } + out = *Clover * in; + } + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = adj(*Clover) * in; + } + } + else + { + if (in._grid->_isCheckerBoarded) + { + + if (in.checkerboard == Odd) + { + // std::cout << "Calling clover term Odd" << std::endl; + Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; + } + else + { + // std::cout << "Calling clover term Even" << std::endl; + Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; + } + out = *Clover * in; + // std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; + } + else + { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = *Clover * in; + } + } + +} // MooeeInternal + + +// Derivative parts +template +void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) +{ + assert(0); +} + +// Derivative parts +template +void WilsonCloverFermion::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) +{ + assert(0); // not implemented yet +} + +FermOpTemplateInstantiate(WilsonCloverFermion); +AdjointFermOpTemplateInstantiate(WilsonCloverFermion); +TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion); +//GparityFermOpTemplateInstantiate(WilsonCloverFermion); +} +} diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h new file mode 100644 index 00000000..268564c0 --- /dev/null +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -0,0 +1,366 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h + + Copyright (C) 2017 + + Author: Guido Cossu + Author: David Preti <> + + 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 */ + +#ifndef GRID_QCD_WILSON_CLOVER_FERMION_H +#define GRID_QCD_WILSON_CLOVER_FERMION_H + +#include + +namespace Grid +{ +namespace QCD +{ + +/////////////////////////////////////////////////////////////////// +// Wilson Clover +// +// Operator ( with anisotropy coefficients): +// +// Q = 1 + (Nd-1)/xi_0 + m +// + W_t + (nu/xi_0) * W_s +// - 1/2*[ csw_t * sum_s (sigma_ts F_ts) + (csw_s/xi_0) * sum_ss (sigma_ss F_ss) ] +// +// s spatial, t temporal directions. +// where W_t and W_s are the temporal and spatial components of the +// Wilson Dirac operator +// +// csw_r = csw_t to recover the isotropic version +////////////////////////////////////////////////////////////////// + +template +class WilsonCloverFermion : public WilsonFermion +{ +public: + // Types definitions + INHERIT_IMPL_TYPES(Impl); + template + using iImplClover = iScalar, Ns>>; + typedef iImplClover SiteCloverType; + typedef Lattice CloverFieldType; + +public: + typedef WilsonFermion WilsonBase; + + virtual void Instantiatable(void){}; + // Constructors + WilsonCloverFermion(GaugeField &_Umu, GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, + const RealD _mass, + const RealD _csw_r = 0.0, + const RealD _csw_t = 0.0, + const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(), + const ImplParams &impl_p = ImplParams()) : WilsonFermion(_Umu, + Fgrid, + Hgrid, + _mass, impl_p, clover_anisotropy), + CloverTerm(&Fgrid), + CloverTermInv(&Fgrid), + CloverTermEven(&Hgrid), + CloverTermOdd(&Hgrid), + CloverTermInvEven(&Hgrid), + CloverTermInvOdd(&Hgrid), + CloverTermDagEven(&Hgrid), + CloverTermDagOdd(&Hgrid), + CloverTermInvDagEven(&Hgrid), + CloverTermInvDagOdd(&Hgrid) + { + assert(Nd == 4); // require 4 dimensions + + if (clover_anisotropy.isAnisotropic) + { + csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0; + diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); + } + else + { + csw_r = _csw_r * 0.5; + diag_mass = 4.0 + _mass; + } + csw_t = _csw_t * 0.5; + + if (csw_r == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl; + if (csw_t == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl; + + ImportGauge(_Umu); + } + + virtual RealD M(const FermionField &in, FermionField &out); + virtual RealD Mdag(const FermionField &in, FermionField &out); + + virtual void Mooee(const FermionField &in, FermionField &out); + virtual void MooeeDag(const FermionField &in, FermionField &out); + virtual void MooeeInv(const FermionField &in, FermionField &out); + 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 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); + + // Guido: 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); + //std::cout << "Lambda:" << Lambda << std::endl; + + 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; + + RealD factor; + if (nu == 4 || mu == 4) + { + factor = 2.0 * csw_t; + } + else + { + factor = 2.0 * csw_r; + } + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= factor*Cmunu(U, lambda, mu, nu); // checked + count++; + } + + pokeLorentz(clover_force, U[mu] * force_mu, mu); + } + //clover_force *= csw; + 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::ShiftStaple(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::ShiftStaple(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_r; // Clover coefficient - spatial + RealD csw_t; // Clover coefficient - temporal + RealD diag_mass; // Mass term + CloverFieldType CloverTerm, CloverTermInv; // Clover term + CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO + CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO + CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO + CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO + + // eventually these can be compressed into 6x6 blocks instead of the 12x12 + // using the DeGrand-Rossi basis for the gamma matrices + CloverFieldType fillCloverYZ(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()()); + T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()()); + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverXZ(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = -F._odata[i]()(); + T._odata[i]()(1, 0) = F._odata[i]()(); + T._odata[i]()(2, 3) = -F._odata[i]()(); + T._odata[i]()(3, 2) = F._odata[i]()(); + } + + return T; + } + + CloverFieldType fillCloverXY(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + + T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()()); + T._odata[i]()(1, 1) = timesI(F._odata[i]()()); + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverXT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = timesI(F._odata[i]()()); + T._odata[i]()(1, 0) = timesI(F._odata[i]()()); + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverYT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 1) = -(F._odata[i]()()); + T._odata[i]()(1, 0) = (F._odata[i]()()); + T._odata[i]()(2, 3) = (F._odata[i]()()); + T._odata[i]()(3, 2) = -(F._odata[i]()()); + } + + return T; + } + + CloverFieldType fillCloverZT(const GaugeLinkField &F) + { + CloverFieldType T(F._grid); + T = zero; + PARALLEL_FOR_LOOP + for (int i = 0; i < CloverTerm._grid->oSites(); i++) + { + T._odata[i]()(0, 0) = timesI(F._odata[i]()()); + T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); + } + + return T; + } +}; +} +} + +#endif // GRID_QCD_WILSON_CLOVER_FERMION_H diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 1a020e8a..dfaa6758 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -47,7 +47,8 @@ int WilsonFermionStatic::HandOptDslash; template WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p) + const ImplParams &p, + const WilsonAnisotropyCoefficients &anis) : Kernels(p), _grid(&Fgrid), _cbgrid(&Hgrid), @@ -60,16 +61,41 @@ WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, Umu(&Fgrid), UmuEven(&Hgrid), UmuOdd(&Hgrid), - _tmp(&Hgrid) + _tmp(&Hgrid), + anisotropyCoeff(anis) { // Allocate the required comms buffer ImportGauge(_Umu); + if (anisotropyCoeff.isAnisotropic){ + diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0); + } else { + diag_mass = 4.0 + mass; + } + + } template void WilsonFermion::ImportGauge(const GaugeField &_Umu) { GaugeField HUmu(_Umu._grid); - HUmu = _Umu * (-0.5); + + //Here multiply the anisotropy coefficients + if (anisotropyCoeff.isAnisotropic) + { + + for (int mu = 0; mu < Nd; mu++) + { + GaugeLinkField U_dir = (-0.5)*PeekIndex(_Umu, mu); + if (mu != anisotropyCoeff.t_direction) + U_dir *= (anisotropyCoeff.nu / anisotropyCoeff.xi_0); + + PokeIndex(HUmu, U_dir, mu); + } + } + else + { + HUmu = _Umu * (-0.5); + } Impl::DoubleStore(GaugeGrid(), Umu, HUmu); pickCheckerboard(Even, UmuEven, Umu); pickCheckerboard(Odd, UmuOdd, Umu); @@ -83,14 +109,14 @@ template RealD WilsonFermion::M(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerNo); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template RealD WilsonFermion::Mdag(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; Dhop(in, out, DaggerYes); - return axpy_norm(out, 4 + mass, in, out); + return axpy_norm(out, diag_mass, in, out); } template @@ -114,7 +140,7 @@ void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - typename FermionField::scalar_type scal(4.0 + mass); + typename FermionField::scalar_type scal(diag_mass); out = scal * in; } @@ -127,7 +153,7 @@ void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { template void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { out.checkerboard = in.checkerboard; - out = (1.0/(4.0+mass))*in; + out = (1.0/(diag_mass))*in; } template @@ -204,7 +230,7 @@ void WilsonFermion::DerivInternal(StencilImpl &st, DoubledGaugeField &U, FermionField Btilde(B._grid); FermionField Atilde(B._grid); - Atilde = A; + Atilde = A;//redundant st.HaloExchange(B, compressor); @@ -393,7 +419,7 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); Lattice> ph(_grid), coor(_grid); - Complex i(0.0,1.0); + ComplexD i(0.0,1.0); PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLt = GridDefaultLatt()[Tp]; @@ -405,7 +431,7 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, LatticeCoordinate(coor, mu); ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu]))); } - ph = exp((Real)(2*M_PI)*i*ph); + ph = exp((RealD)(2*M_PI)*i*ph); q_out = zero; LatticeInteger coords(_grid); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index feba40ed..c0db827e 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -44,6 +44,21 @@ class WilsonFermionStatic { static const int npoint = 8; }; + struct WilsonAnisotropyCoefficients: Serializable + { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients, + bool, isAnisotropic, + int, t_direction, + double, xi_0, + double, nu); + + WilsonAnisotropyCoefficients(): + isAnisotropic(false), + t_direction(Nd-1), + xi_0(1.0), + nu(1.0){} +}; + template class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { public: @@ -65,8 +80,8 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // override multiply; cut number routines if pass dagger argument // and also make interface more uniformly consistent ////////////////////////////////////////////////////////////////// - RealD M(const FermionField &in, FermionField &out); - RealD Mdag(const FermionField &in, FermionField &out); + virtual RealD M(const FermionField &in, FermionField &out); + virtual RealD Mdag(const FermionField &in, FermionField &out); ///////////////////////////////////////////////////////// // half checkerboard operations @@ -117,8 +132,9 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // Constructor WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, - GridRedBlackCartesian &Hgrid, RealD _mass, - const ImplParams &p = ImplParams()); + GridRedBlackCartesian &Hgrid, RealD _mass, + const ImplParams &p = ImplParams(), + const WilsonAnisotropyCoefficients &anis = WilsonAnisotropyCoefficients() ); // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); @@ -130,6 +146,7 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { // protected: public: RealD mass; + RealD diag_mass; GridBase *_grid; GridBase *_cbgrid; @@ -146,6 +163,8 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; + + WilsonAnisotropyCoefficients anisotropyCoeff; /////////////////////////////////////////////////////////////// // Conserved current utilities diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 393ee7f3..3e58fed6 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -793,7 +793,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, Lattice> ph(FermionGrid()), coor(FermionGrid()); PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), tmp(FermionGrid()); - Complex i(0.0, 1.0); + ComplexD i(0.0, 1.0); unsigned int tshift = (mu == Tp) ? 1 : 0; unsigned int LLs = q_in._grid->_rdimensions[0]; unsigned int LLt = GridDefaultLatt()[Tp]; @@ -806,7 +806,7 @@ void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, LatticeCoordinate(coor, nu + 1); ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu]))); } - ph = exp((Real)(2*M_PI)*i*ph); + ph = exp((RealD)(2*M_PI)*i*ph); q_out = zero; LatticeInteger coords(_FourDimGrid); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index ed8d6be9..2369c98d 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -55,7 +55,7 @@ template class WilsonKernels : public FermionOperator , public public: template - typename std::enable_if::type + typename std::enable_if::type DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) { @@ -99,7 +99,7 @@ public: } template - typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type + typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool, void>::type DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1 ) { // no kernel choice @@ -116,7 +116,7 @@ public: } template - typename std::enable_if::type + typename std::enable_if::type DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) { @@ -161,7 +161,7 @@ public: } template - typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool,void>::type + typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool,void>::type DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,SiteHalfSpinor * buf, int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) { diff --git a/lib/qcd/action/fermion/WilsonKernelsHand.cc b/lib/qcd/action/fermion/WilsonKernelsHand.cc index 80b81714..aa6b5f6b 100644 --- a/lib/qcd/action/fermion/WilsonKernelsHand.cc +++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc @@ -946,5 +946,6 @@ INSTANTIATE_THEM(DomainWallVec5dImplFH); INSTANTIATE_THEM(DomainWallVec5dImplDF); INSTANTIATE_THEM(ZDomainWallVec5dImplFH); INSTANTIATE_THEM(ZDomainWallVec5dImplDF); - +INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplF); +INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplD); }} diff --git a/lib/qcd/modules/ObservableModules.h b/lib/qcd/modules/ObservableModules.h index 24511617..fbffc236 100644 --- a/lib/qcd/modules/ObservableModules.h +++ b/lib/qcd/modules/ObservableModules.h @@ -92,6 +92,19 @@ class PlaquetteMod: public ObservableModule, NoParameters> PlaquetteMod(): ObsBase(NoParameters()){} }; +template < class Impl > +class PolyakovMod: public ObservableModule, NoParameters>{ + typedef ObservableModule, NoParameters> ObsBase; + using ObsBase::ObsBase; // for constructors + + // acquire resource + virtual void initialize(){ + this->ObservablePtr.reset(new PolyakovLogger()); + } + public: + PolyakovMod(): ObsBase(NoParameters()){} +}; + template < class Impl > class TopologicalChargeMod: public ObservableModule, TopologyObsParameters>{ diff --git a/lib/qcd/observables/hmc_observable.h b/lib/qcd/observables/hmc_observable.h index db629ce7..fcf11774 100644 --- a/lib/qcd/observables/hmc_observable.h +++ b/lib/qcd/observables/hmc_observable.h @@ -45,5 +45,7 @@ class HmcObservable { #include "plaquette.h" #include "topological_charge.h" +#include "polyakov_loop.h" + #endif // HMC_OBSERVABLE_H diff --git a/lib/qcd/observables/polyakov_loop.h b/lib/qcd/observables/polyakov_loop.h new file mode 100644 index 00000000..d708b474 --- /dev/null +++ b/lib/qcd/observables/polyakov_loop.h @@ -0,0 +1,68 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/modules/polyakov_line.h + +Copyright (C) 2017 + +Author: David Preti + +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 */ + +#ifndef HMC_POLYAKOV_H +#define HMC_POLYAKOV_H + +namespace Grid { +namespace QCD { + +// this is only defined for a gauge theory +template +class PolyakovLogger : public HmcObservable { + public: + // here forces the Impl to be of gauge fields + // if not the compiler will complain + INHERIT_GIMPL_TYPES(Impl); + + // necessary for HmcObservable compatibility + typedef typename Impl::Field Field; + + void TrajectoryComplete(int traj, + Field &U, + GridSerialRNG &sRNG, + GridParallelRNG &pRNG) { + + ComplexD polyakov = WilsonLoops::avgPolyakovLoop(U); + + int def_prec = std::cout.precision(); + + std::cout << GridLogMessage + << std::setprecision(std::numeric_limits::digits10 + 1) + << "Polyakov Loop: [ " << traj << " ] "<< polyakov << std::endl; + + std::cout.precision(def_prec); + + } +}; + +} // namespace QCD +} // namespace Grid + +#endif // HMC_POLYAKOV_H diff --git a/lib/qcd/representations/adjoint.h b/lib/qcd/representations/adjoint.h index 078d12a1..052cd7a8 100644 --- a/lib/qcd/representations/adjoint.h +++ b/lib/qcd/representations/adjoint.h @@ -23,6 +23,7 @@ class AdjointRep { typedef typename SU_Adjoint::LatticeAdjMatrix LatticeMatrix; typedef typename SU_Adjoint::LatticeAdjField LatticeField; static const int Dimension = ncolour * ncolour - 1; + static const bool isFundamental = false; LatticeField U; diff --git a/lib/qcd/representations/fundamental.h b/lib/qcd/representations/fundamental.h index db52d893..9f039a07 100644 --- a/lib/qcd/representations/fundamental.h +++ b/lib/qcd/representations/fundamental.h @@ -19,6 +19,7 @@ template class FundamentalRep { public: static const int Dimension = ncolour; + static const bool isFundamental = true; // typdef to be used by the Representations class in HMC to get the // types for the higher representation fields diff --git a/lib/qcd/representations/two_index.h b/lib/qcd/representations/two_index.h index 082a52a5..2c7e8b3a 100644 --- a/lib/qcd/representations/two_index.h +++ b/lib/qcd/representations/two_index.h @@ -29,6 +29,7 @@ class TwoIndexRep { typedef typename SU_TwoIndex::LatticeTwoIndexMatrix LatticeMatrix; typedef typename SU_TwoIndex::LatticeTwoIndexField LatticeField; static const int Dimension = ncolour * (ncolour + S) / 2; + static const bool isFundamental = false; LatticeField U; diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 90d905d3..cdd76ecc 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -123,6 +123,28 @@ public: return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } + + ////////////////////////////////////////////////// + // average over all x,y,z the temporal loop + ////////////////////////////////////////////////// + static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4 + GaugeMat Ut(Umu._grid), P(Umu._grid); + ComplexD out; + int T = Umu._grid->GlobalDimensions()[3]; + int X = Umu._grid->GlobalDimensions()[0]; + int Y = Umu._grid->GlobalDimensions()[1]; + int Z = Umu._grid->GlobalDimensions()[2]; + + Ut = peekLorentz(Umu,3); //Select temporal direction + P = Ut; + for (int t=1;t--+ +-->--+(x) + // (x)+-->--+ +-->--+(x) - h.c. // | | | | // +--<--+ +--<--+ @@ -335,7 +359,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GaugeMat v = Vup - Vdn; GaugeMat u = PeekIndex(Umu, mu); // some redundant copies GaugeMat vu = v*u; - FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); + //FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); + FS = (u*v + Cshift(vu, mu, -1)); + FS = 0.125*(FS - adj(FS)); } static Real TopologicalCharge(GaugeLorentz &U){ @@ -360,6 +386,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { return TensorRemove(Tq).real(); } + ////////////////////////////////////////////////////// // Similar to above for rectangle is required ////////////////////////////////////////////////////// diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc new file mode 100644 index 00000000..9281e298 --- /dev/null +++ b/tests/core/Test_wilson_clover.cc @@ -0,0 +1,357 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_wilson.cc + + Copyright (C) 2015 + + 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 + +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::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl; + + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); + + typedef typename WilsonCloverFermionR::FermionField FermionField; + typename WilsonCloverFermionR::ImplParams params; + WilsonAnisotropyCoefficients anis; + + FermionField src(&Grid); + random(pRNG, src); + FermionField result(&Grid); + result = zero; + FermionField result2(&Grid); + result2 = zero; + FermionField ref(&Grid); + ref = zero; + FermionField tmp(&Grid); + tmp = zero; + FermionField err(&Grid); + err = zero; + FermionField err2(&Grid); + err2 = zero; + FermionField phi(&Grid); + random(pRNG, phi); + FermionField chi(&Grid); + random(pRNG, chi); + LatticeGaugeField Umu(&Grid); + SU3::HotConfiguration(pRNG, Umu); + std::vector U(4, &Grid); + + double volume = 1; + for (int mu = 0; mu < Nd; mu++) + { + volume = volume * latt_size[mu]; + } + + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + + WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + //Dwc.ImportGauge(Umu); // not necessary, included in the constructor + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + FermionField src_e(&RBGrid); + FermionField src_o(&RBGrid); + FermionField r_e(&RBGrid); + FermionField r_o(&RBGrid); + FermionField r_eo(&Grid); + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd, src_o, src); + + Dwc.Meooe(src_e, r_o); + std::cout << GridLogMessage << "Applied Meo" << std::endl; + Dwc.Meooe(src_o, r_e); + std::cout << GridLogMessage << "Applied Moe" << std::endl; + Dwc.Dhop(src, ref, DaggerNo); + + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; + std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + FermionField chi_e(&RBGrid); + FermionField chi_o(&RBGrid); + + FermionField dchi_e(&RBGrid); + FermionField dchi_o(&RBGrid); + + FermionField phi_e(&RBGrid); + FermionField phi_o(&RBGrid); + + FermionField dphi_e(&RBGrid); + FermionField dphi_o(&RBGrid); + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Meooe(chi_e, dchi_o); + Dwc.Meooe(chi_o, dchi_e); + Dwc.MeooeDag(phi_e, dphi_o); + Dwc.MeooeDag(phi_o, dphi_e); + + ComplexD pDce = innerProduct(phi_e, dchi_e); + ComplexD pDco = innerProduct(phi_o, dchi_o); + ComplexD cDpe = innerProduct(chi_e, dphi_e); + ComplexD cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv Mee = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.Mooee(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.Mooee(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInvDag(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = zero; + phi = zero; + tmp = zero; + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + ////////////////////// Gauge Transformation + std::vector seeds2({5, 6, 7, 8}); + GridParallelRNG pRNG2(&Grid); + pRNG2.SeedFixedIntegers(seeds2); + LatticeColourMatrix Omega(&Grid); + LatticeColourMatrix ShiftedOmega(&Grid); + LatticeGaugeField U_prime(&Grid); + U_prime = zero; + LatticeColourMatrix U_prime_mu(&Grid); + U_prime_mu = zero; + SU::LieRandomize(pRNG2, Omega, 1.0); + for (int mu = 0; mu < Nd; mu++) + { + U[mu] = peekLorentz(Umu, mu); + ShiftedOmega = Cshift(Omega, mu, 1); + U_prime_mu = Omega * U[mu] * adj(ShiftedOmega); + pokeLorentz(U_prime, U_prime_mu, mu); + } + ///////////////// + + WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + Dwc_prime.ImportGauge(U_prime); + + tmp = Omega * src; + pickCheckerboard(Even, src_e, tmp); + pickCheckerboard(Odd, src_o, tmp); + + Dwc_prime.Mooee(src_e, phi_e); + Dwc_prime.Mooee(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "=================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = zero; + phi = zero; + + WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params); + Dw.ImportGauge(Umu); + + Dw.M(src, result); + Dwc.M(src, chi); + + Dwc_prime.M(Omega * src, phi); + + WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params); + Dw_prime.ImportGauge(U_prime); + Dw_prime.M(Omega * src, result2); + + err = chi - adj(Omega) * phi; + err2 = result - adj(Omega) * result2; + std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff WilsonClover " << norm2(err2) << std::endl; + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = zero; + phi = zero; + err = zero; + WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0 + Dwc_csw0.ImportGauge(Umu); + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dw.Mooee(src_e, chi_e); + Dw.Mooee(src_o, chi_o); + Dwc_csw0.Mooee(src_e, phi_e); + Dwc_csw0.Mooee(src_o, phi_o); + + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + err = chi - phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = zero; + phi = zero; + err = zero; + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o) + + Dwc.M(src, ref); // Reference result from the unpreconditioned operator + + // EO matrix + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + Dwc.Meooe(src_o, phi_e); + Dwc.Meooe(src_e, phi_o); + + phi_o += chi_o; + phi_e += chi_e; + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = ref - phi; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff :" << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff :" << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff :" << norm2(err) << std::endl; + + Grid_finalize(); +} 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..f7090845 --- /dev/null +++ b/tests/forces/Test_wilsonclover_force.cc @@ -0,0 +1,194 @@ +/************************************************************************************* + + 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, 30, 50}); + + 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); + //pRNG.SeedFixedIntegers(seeds); + + LatticeFermion phi(&Grid); + gaussian(pRNG, phi); + LatticeFermion Mphi(&Grid); + LatticeFermion MphiPrime(&Grid); + + LatticeGaugeField U(&Grid); + + std::vector site = {0, 0, 0, 0}; + SU3::HotConfiguration(pRNG, U); + //SU3::ColdConfiguration(pRNG, U);// Clover term zero + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass = 0.1; + Real csw = 1.0; + WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw, csw); + Dw.ImportGauge(U); + 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; + ///////////////////////////////////////////// + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.00005; + 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; ss < mom._grid->oSites(); ss++) + { + 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 + ////////////////////////////////////////////// + + LatticeComplex dS(&Grid); + dS = zero; + LatticeComplex dSmom(&Grid); + dSmom = zero; + LatticeComplex dSmom2(&Grid); + dSmom2 = zero; + + 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" << 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 (S' - S) :" << Sprime - S << std::endl; + std::cout << GridLogMessage << "predict dS (force) :" << 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(); +} diff --git a/tests/hadrons/Test_hadrons_2AS_spectrum.cc b/tests/hadrons/Test_hadrons_2AS_spectrum.cc new file mode 100644 index 00000000..2f519834 --- /dev/null +++ b/tests/hadrons/Test_hadrons_2AS_spectrum.cc @@ -0,0 +1,168 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015 + + Author: Antonin Portelli + + 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. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + + + BEGIN_HADRONS_NAMESPACE + BEGIN_MODULE_NAMESPACE(MFermion) + MODULE_REGISTER_NS(GaugeProp2AS, TGaugeProp, MFermion); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSource) + MODULE_REGISTER_NS(Point2AS, TPoint, MSource); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MContraction) + MODULE_REGISTER_NS(Meson2AS, ARG(TMeson), MContraction); +// MODULE_REGISTER_NS(BaryonMultirep, ARG(TBaryon), MContraction); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSink) + MODULE_REGISTER_NS(ScalarPoint2AS, TPoint, MSink); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MSolver) + MODULE_REGISTER_NS(RBPrecCG2AS, TRBPrecCG, MSolver); + END_MODULE_NAMESPACE + BEGIN_MODULE_NAMESPACE(MAction) + MODULE_REGISTER_NS(WilsonClover2AS, TWilsonClover, MAction); + END_MODULE_NAMESPACE + END_HADRONS_NAMESPACE + + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"l", "s"}; + std::vector mass = {-0.01, -0.04}; + double csw = 1.0; + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + MSource::Point2AS::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::ScalarPoint2AS::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::WilsonClover2AS::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.mass = mass[i]; + actionPar.csw_r = csw; + actionPar.csw_t = csw; + actionPar.clover_anisotropy.isAnisotropic= false; + actionPar.clover_anisotropy.t_direction = Nd-1 ; + actionPar.clover_anisotropy.xi_0 = 1.0 ; + actionPar.clover_anisotropy.nu = 1.0 ; + actionPar.boundary = boundary; + application.createModule("WilsonClover2AS_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG2AS::Par solverPar; + solverPar.action = "WilsonClover2AS_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp2AS::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], quarkPar); + quarkPar.source = "z2"; + application.createModule("QZ2_" + flavour[i], quarkPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + MContraction::Meson2AS::Par mesPar; + + mesPar.output = "mesons2AS/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + + // mesPar.output = "mesons2AS/Z2_" + flavour[i] + flavour[j]; + // mesPar.q1 = "QZ2_" + flavour[i]; + // mesPar.q2 = "QZ2_" + flavour[j]; + // mesPar.gammas = "all"; + // mesPar.sink = "sink"; + // application.createModule("meson_Z2_" + // + flavour[i] + flavour[j], + // mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + for (unsigned int k = j; k < flavour.size(); ++k) + { + MContraction::Baryon::Par barPar; + + barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k]; + barPar.q1 = "Qpt_" + flavour[i]; + barPar.q2 = "Qpt_" + flavour[j]; + barPar.q3 = "Qpt_" + flavour[k]; + application.createModule( + "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); + } + + // execution + application.saveParameterFile("spectrum.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc index eac065e9..cf3d2590 100644 --- a/tests/hadrons/Test_hadrons_quark.cc +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -130,8 +130,8 @@ int main(int argc, char **argv) for (int c = 0; c < Nc; ++c) { ref = prop; - PropToFerm(ferm, prop, s, c); - FermToProp(prop, ferm, s, c); + PropToFerm(ferm, prop, s, c); + FermToProp(prop, ferm, s, c); std::cout << "Spin = " << s << ", Colour = " << c << std::endl; ref -= prop; diff --git a/tests/hadrons/Test_hadrons_wilsonFund.cc b/tests/hadrons/Test_hadrons_wilsonFund.cc new file mode 100644 index 00000000..083d3b8c --- /dev/null +++ b/tests/hadrons/Test_hadrons_wilsonFund.cc @@ -0,0 +1,157 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_spectrum.cc + + Copyright (C) 2015 + + Author: Antonin Portelli + + 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. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"l"}; + std::vector mass = {-0.1}; + double csw = 0.0; + + // global parameters + Application::GlobalPar globalPar; + + globalPar.trajCounter.start = 309; + globalPar.trajCounter.end = 310; + globalPar.trajCounter.step = 1; + + globalPar.seed = "1 2 3 4"; + + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + + // sources + //MSource::Z2::Par z2Par; + //z2Par.tA = 0; + //z2Par.tB = 0; + //application.createModule("z2", z2Par); + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::WilsonClover::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + actionPar.csw_r = csw; + actionPar.csw_t = csw; + + // !!!!! Check if Anisotropy works !!!!! + actionPar.clover_anisotropy.isAnisotropic= false; + actionPar.clover_anisotropy.t_direction = 3 ; // Explicit for D=4 + actionPar.clover_anisotropy.xi_0 = 1.0 ; + actionPar.clover_anisotropy.nu = 1.0 ; + + application.createModule("WilsonClover_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "WilsonClover_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], quarkPar); + // quarkPar.source = "z2"; + // application.createModule("QZ2_" + flavour[i], quarkPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "Fund_mesons/pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + // mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + // mesPar.q1 = "QZ2_" + flavour[i]; + // mesPar.q2 = "QZ2_" + flavour[j]; + // mesPar.gammas = "all"; + // mesPar.sink = "sink"; + // application.createModule("meson_Z2_" + // + flavour[i] + flavour[j], + // mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + for (unsigned int k = j; k < flavour.size(); ++k) + { + MContraction::Baryon::Par barPar; + + barPar.output = "Fund_baryons/pt_" + flavour[i] + flavour[j] + flavour[k]; + barPar.q1 = "Qpt_" + flavour[i]; + barPar.q2 = "Qpt_" + flavour[j]; + barPar.q3 = "Qpt_" + flavour[k]; + application.createModule( + "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); + } + + // execution + application.saveParameterFile("WilsonClover_spectrum.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc new file mode 100644 index 00000000..fbb8a4a3 --- /dev/null +++ b/tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc @@ -0,0 +1,139 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2016 + +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 + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + // Possibile to create the module by hand + // hardcoding parameters or using a Reader + + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6 ; + WilsonGaugeActionR Waction(beta); + + // temporarily need a gauge field + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + + ConjugateGradient CG(1.0e-8, 2000); + + TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + + + + + + + diff --git a/tests/hmc/Test_hmc_WC2ASFG_Production.cc b/tests/hmc/Test_hmc_WC2ASFG_Production.cc new file mode 100644 index 00000000..d255ab5d --- /dev/null +++ b/tests/hmc/Test_hmc_WC2ASFG_Production.cc @@ -0,0 +1,213 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonTwoIndexAntiSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverTwoIndexAntiSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + TwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + diff --git a/tests/hmc/Test_hmc_WC2SFG_Production.cc b/tests/hmc/Test_hmc_WC2SFG_Production.cc new file mode 100644 index 00000000..8d5fc458 --- /dev/null +++ b/tests/hmc/Test_hmc_WC2SFG_Production.cc @@ -0,0 +1,212 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, TwoIndexSymmetricRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonTwoIndexSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverTwoIndexSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + TwoIndexSymmetricRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCFG_Production.cc b/tests/hmc/Test_hmc_WCFG_Production.cc new file mode 100644 index 00000000..895e4f81 --- /dev/null +++ b/tests/hmc/Test_hmc_WCFG_Production.cc @@ -0,0 +1,210 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + //typedef TopologicalChargeMod QObs; + //TopologyObsParameters TopParams(Reader); + //TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc b/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc new file mode 100644 index 00000000..aa5cce85 --- /dev/null +++ b/tests/hmc/Test_hmc_WCMixedRepFG_Production.cc @@ -0,0 +1,224 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc + +Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +Author: neo +Author: paboyle + +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 "Grid/Grid.h" + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonCloverFund, + FermionParameters, WilsonCloverAS) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; +} + + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + // Here change the allowed (higher) representations + typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation> TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; + + typedef WilsonImplR FundImplPolicy; + typedef WilsonCloverFermionR FundFermionAction; + typedef typename FundFermionAction::FermionField FundFermionField; + + typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; + typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; + typedef typename ASymmFermionAction::FermionField ASymmFermionField; + + typedef Grid::JSONReader Serialiser; + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearingFund = MyParams.WilsonCloverFund.ApplySmearing; + bool ApplySmearingAS = MyParams.WilsonCloverAS.ApplySmearing; + + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + typedef TopologicalChargeMod QObs; + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + FundamentalRepresentation::LatticeField UF(GridPtr); + TwoIndexAntiSymmetricRepresentation::LatticeField UAS(GridPtr); + + + Real Fundmass = MyParams.WilsonCloverFund.mass; + Real Fundcsw = MyParams.WilsonCloverFund.csw; + Real ASmass = MyParams.WilsonCloverAS.mass; + Real AScsw = MyParams.WilsonCloverAS.csw; + + + + std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; + std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl; + + + FundFermionAction FundFermOp(UF, *GridPtr, *GridRBPtr, Fundmass, Fundcsw, Fundcsw); + ConjugateGradient CG_Fund(MyParams.WilsonCloverFund.StoppingCondition, MyParams.WilsonCloverFund.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2_Fund(FundFermOp, CG_Fund, CG_Fund); + + ASymmFermionAction ASFermOp(UAS, *GridPtr, *GridRBPtr, ASmass, AScsw, AScsw); + ConjugateGradient CG_AS(MyParams.WilsonCloverAS.StoppingCondition, MyParams.WilsonCloverAS.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2_AS(ASFermOp, CG_AS, CG_AS); + + Nf2_Fund.is_smeared = ApplySmearingFund; + Nf2_AS.is_smeared = ApplySmearingAS; + + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2_Fund); + Level1.push_back(&Nf2_AS); + + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; +/* + if (ApplySmearingFund || ApplySmearingAS){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } +*/ + TheHMC.Run(); + + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WCadjFG_Production.cc b/tests/hmc/Test_hmc_WCadjFG_Production.cc new file mode 100644 index 00000000..48cea756 --- /dev/null +++ b/tests/hmc/Test_hmc_WCadjFG_Production.cc @@ -0,0 +1,213 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.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 + + +namespace Grid{ + struct FermionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters, + double, mass, + double, csw, + double, StoppingCondition, + int, MaxCGIterations, + bool, ApplySmearing); + }; + + + struct WilsonCloverHMCParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters, + double, gauge_beta, + FermionParameters, WilsonClover) + + template + WilsonCloverHMCParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + }; + + struct SmearingParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters, + double, rho, + Integer, Nsmear) + + template + SmearingParameters(Reader& Reader){ + read(Reader, "StoutSmearing", *this); + } + + }; + + +} + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunnerHirep HMCWrapper; // Uses the default minimum norm + typedef WilsonAdjImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions + typedef WilsonCloverAdjFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...) + typedef typename FermionAction::FermionField FermionField; + typedef Grid::JSONReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.ReadCommandLine(argc, argv); + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + WilsonCloverHMCParameters MyParams(Reader); + + // Apply smearing to the fermionic action + bool ApplySmearing = MyParams.WilsonClover.ApplySmearing; + + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + + /* + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + */ + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + /* + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + */ + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + + typedef TopologicalChargeMod QObs; + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + + //RealD beta = 5.6; + WilsonGaugeActionR Waction(MyParams.gauge_beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + AdjointRepresentation::LatticeField U(GridPtr); + + //Real mass = 0.01; + //Real csw = 1.0; + + Real mass = MyParams.WilsonClover.mass; + Real csw = MyParams.WilsonClover.csw; + + std::cout << "mass and csw" << mass << " and " << csw << std::endl; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw); + ConjugateGradient CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations); + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = ApplySmearing; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + + TheHMC.Parameters.initialize(Reader); + //TheHMC.Parameters.MD.MDsteps = 20; + //TheHMC.Parameters.MD.trajL = 1.0; + + if (ApplySmearing){ + SmearingParameters SmPar(Reader); + //double rho = 0.1; // smearing parameter + //int Nsmear = 3; // number of smearing levels + Smear_Stout Stout(SmPar.rho); + SmearedConfiguration SmearingPolicy(GridPtr, SmPar.Nsmear, Stout); + TheHMC.Run(SmearingPolicy); // for smearing + } else { + TheHMC.Run(); // no smearing + } + + //TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + //TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + diff --git a/tests/hmc/Test_hmc_WG_Production.cc b/tests/hmc/Test_hmc_WG_Production.cc new file mode 100644 index 00000000..b99446d5 --- /dev/null +++ b/tests/hmc/Test_hmc_WG_Production.cc @@ -0,0 +1,114 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2015 + +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 + +namespace Grid{ + struct ActionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(ActionParameters, + double, beta) + + + template + ActionParameters(Reader& Reader){ + read(Reader, "Action", *this); + } + + }; + +} + + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + GridLogLayout(); + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + HMCWrapper TheHMC; + typedef Grid::JSONReader Serialiser; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + // Reader, file should come from command line + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + + + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + typedef TopologicalChargeMod QObs; + TheHMC.Resources.AddObservable(); + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + ActionParameters WilsonPar(Reader); + //RealD beta = 6.4 ; + WilsonGaugeActionR Waction(WilsonPar.beta); + + ActionLevel Level1(1); + Level1.push_back(&Waction); + //Level1.push_back(WGMod.getPtr()); + TheHMC.TheAction.push_back(Level1); + ///////////////////////////////////////////////////////////// + + // HMC parameters are serialisable + TheHMC.Parameters.initialize(Reader); + + //TheHMC.Parameters.MD.MDsteps = 17; + //TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.Run(); // no smearing + + Grid_finalize(); + +} // main diff --git a/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc new file mode 100644 index 00000000..8d5701b8 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc @@ -0,0 +1,129 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.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 + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + + typedef PolyakovMod PolyakovObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6; + WilsonGaugeActionR Waction(beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + ConjugateGradient CG(1.0e-8, 5000); + + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/lanczos/Test_WCMultiRep_lanczos.cc b/tests/lanczos/Test_WCMultiRep_lanczos.cc new file mode 100644 index 00000000..b6d69aee --- /dev/null +++ b/tests/lanczos/Test_WCMultiRep_lanczos.cc @@ -0,0 +1,178 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_dwf_lanczos.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; + + +//typedef WilsonCloverFermionR FermionOp; +//typedef typename WilsonFermionR::FermionField FermionField; + +typedef WilsonImplR FundImplPolicy; +typedef WilsonCloverFermionR FundFermionAction; +typedef typename FundFermionAction::FermionField FundFermionField; + +typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy; +typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction; +typedef typename ASymmFermionAction::FermionField ASymmFermionField; + + +RealD AllZero(RealD x) { return 0.; } + +int main(int argc, char** argv) { + Grid_init(&argc, &argv); + + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = + SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = UGrid; + GridRedBlackCartesian* FrbGrid = UrbGrid; + printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, + FrbGrid); + + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); + RNG5.SeedFixedIntegers(seeds5); + + GridParallelRNG pRNG(UGrid); + GridSerialRNG sRNG; + + FundamentalRepresentation::LatticeField Umu(UGrid); + + TwoIndexAntiSymmetricRepresentation HiRep(UGrid); + TwoIndexAntiSymmetricRepresentation::LatticeField UmuAS(UGrid); + + + CheckpointerParameters CPparams; + + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.format = "IEEE64BIG"; + +//NerscHmcCheckpointer Checkpoint(std::string("ckpoint_lat"), + // std::string("ckpoint_rng"), 1); + +NerscHmcCheckpointer Checkpoint(CPparams); + + int CNFGSTART=1; + int CNFGEND=2; + int CNFGSTEP=1; + + Real Fundmass = -0.1; + Real Fundcsw = 1.0; + Real ASmass = -0.1; + Real AScsw = 1.0; + + std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl; + std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl; + + const int Nstop = 30; + const int Nk = 40; + const int Np = 40; + const int Nm = Nk + Np; + const int MaxIt = 10000; + RealD resid = 1.0e-8; + + for (int cnfg=CNFGSTART;cnfg<=CNFGEND;cnfg+=CNFGSTEP){ + Checkpoint.CheckpointRestore(cnfg,Umu, sRNG, pRNG); + + //SU4::HotConfiguration(RNG4, Umu); // temporary, then read. + + HiRep.update_representation(Umu); + UmuAS = HiRep.U; + + FundFermionAction FundFermOp(Umu,*FGrid,*FrbGrid, Fundmass, Fundcsw, Fundcsw); + MdagMLinearOperator HermOpFund(FundFermOp); /// <----- + + ASymmFermionAction ASFermOp(UmuAS,*FGrid,*FrbGrid, ASmass, AScsw, AScsw); + MdagMLinearOperator HermOpAS(ASFermOp); /// <----- + + std::vector Coeffs{0, -1.}; + Polynomial FundPolyX(Coeffs); + //Chebyshev FundCheb(0.0, 10., 12); + + FunctionHermOp FundPolyXOp(FundPolyX,HermOpFund); + PlainHermOp FundOp (HermOpFund); + + ImplicitlyRestartedLanczos IRL_Fund(FundOp, FundPolyXOp, Nstop, Nk, Nm, + resid, MaxIt); + + Polynomial ASPolyX(Coeffs); + //Chebyshev ASCheb(0.0, 10., 12); + + FunctionHermOp ASPolyXOp(ASPolyX,HermOpAS); + PlainHermOp ASOp (HermOpAS); + + ImplicitlyRestartedLanczos IRL_AS(ASOp, ASPolyXOp, Nstop, Nk, Nm, + resid, MaxIt); + + std::vector Fundeval(Nm); + std::vector ASeval(Nm); + + FundFermionField Fundsrc(FGrid); + ASymmFermionField ASsrc(FGrid); + + gaussian(RNG5, Fundsrc); + gaussian(RNG5, ASsrc); + + std::vector Fundevec(Nm, FGrid); + std::vector ASevec(Nm, FGrid); + + for (int i = 0; i < 1; i++) { + std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i]._grid + << std::endl; + }; + for (int i = 0; i < 1; i++) { + std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i]._grid + << std::endl; + }; + + int FundNconv, ASNconv; + IRL_Fund.calc(Fundeval, Fundevec, Fundsrc, FundNconv); + IRL_AS.calc(ASeval, ASevec, ASsrc, ASNconv); + + for (int i=0;i +#include +#include +#include + +// Mass +double mq = 0.1; + +// Define Wilson Types +typedef Grid::QCD::WilsonImplR::FermionField FermionField; +typedef Grid::QCD::LatticeGaugeField GaugeField; + +enum ChromaAction +{ + Wilson, // Wilson + WilsonClover // CloverFermions +}; + +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; + AnisoParam_t _apar; + _apar.anisoP = true; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; + + 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(2.0); + p.u0 = QDP::Real(1.0); + AnisoParam_t _apar; + _apar.anisoP = true; + _apar.t_dir = 3; // in 4d + _apar.xi_0 = 2.0; + _apar.nu = 1.0; + p.anisoParam = _apar; + + 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 + +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); + + // Fermion field + Grid::gaussian(RNG4, src); + /* + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + + + std::vector x(4); // 4d fermions + std::vector 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 < 4; sp++) + { + for (int j = 0; j < 3; j++) // colours + { + F()(sp)(j) = Grid::Complex(0.0,0.0); + if (((sp == 0)|| (sp==3)) && (j==2)) + { + 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) +{ + 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) + { + WilsonAnisotropyCoefficients anis; + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams iParam; + Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis); + + 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_r = 1.0; + Grid::RealD _csw_t = 2.0; + WilsonAnisotropyCoefficients anis; + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams CloverImplParam; + Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam); + 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); +} + +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); + FermionField only_wilson(UGrid); + FermionField difference(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; + + difference = res_chroma - res_grid; + std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl; + } + } + + std::cout << "Finished test " << std::endl; + + Chroma::finalize(); + } + } +}