diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 33f6b5ea..e75584bc 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -247,8 +247,9 @@ template void WilsonFermion::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { conformable(U._grid, _cbgrid); conformable(U._grid, V._grid); - conformable(U._grid, mat._grid); - + //conformable(U._grid, mat._grid); not general, leaving as a comment (Guido) + // Motivation: look at the SchurDiff operator + assert(V.checkerboard == Even); assert(U.checkerboard == Odd); mat.checkerboard = Odd; @@ -260,7 +261,7 @@ template void WilsonFermion::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { conformable(U._grid, _cbgrid); conformable(U._grid, V._grid); - conformable(U._grid, mat._grid); + //conformable(U._grid, mat._grid); assert(V.checkerboard == Odd); assert(U.checkerboard == Even); diff --git a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h index 367a735d..90b6dbaa 100644 --- a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h +++ b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -67,6 +67,7 @@ namespace Grid{ // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // it is not conformable with the HMC force field + // Case: Ls vectorised fields // INHERIT FROM THE Force field instead GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); GaugeField ForceO(forcecb); @@ -76,9 +77,8 @@ namespace Grid{ // X^dag Der_oe MeeInv Meo Y // Use Mooee as nontrivial but gauge field indept this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied - this->_Mat.MooeeInv(tmp1,tmp2); // even->even + this->_Mat.MooeeInv(tmp1,tmp2); // even->even this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); - // Accumulate X^dag M_oe MeeInv Der_eo Y this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 009899ea..21faa0f7 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -31,82 +31,89 @@ directory #define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H namespace Grid { -namespace QCD { + namespace QCD { -//////////////////////////////////////////////////////////////////////// -// Two flavour pseudofermion action for any EO prec dop -//////////////////////////////////////////////////////////////////////// -template -class TwoFlavourEvenOddPseudoFermionAction - : public Action { - public: - INHERIT_IMPL_TYPES(Impl); + //////////////////////////////////////////////////////////////////////// + // Two flavour pseudofermion action for any EO prec dop + //////////////////////////////////////////////////////////////////////// + template + class TwoFlavourEvenOddPseudoFermionAction + : public Action { + public: + INHERIT_IMPL_TYPES(Impl); - private: - FermionOperator &FermOp; // the basic operator + private: + FermionOperator &FermOp; // the basic operator - OperatorFunction &DerivativeSolver; - OperatorFunction &ActionSolver; + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; - FermionField PhiOdd; // the pseudo fermion field for this trajectory - FermionField PhiEven; // the pseudo fermion field for this trajectory + FermionField PhiOdd; // the pseudo fermion field for this trajectory + FermionField PhiEven; // the pseudo fermion field for this trajectory - public: - ///////////////////////////////////////////////// - // Pass in required objects. - ///////////////////////////////////////////////// - TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, - OperatorFunction &DS, - OperatorFunction &AS) - : FermOp(Op), - DerivativeSolver(DS), - ActionSolver(AS), - PhiEven(Op.FermionRedBlackGrid()), - PhiOdd(Op.FermionRedBlackGrid()) - {}; + public: + ///////////////////////////////////////////////// + // Pass in required objects. + ///////////////////////////////////////////////// + TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, + OperatorFunction &DS, + OperatorFunction &AS) + : FermOp(Op), + DerivativeSolver(DS), + ActionSolver(AS), + PhiEven(Op.FermionRedBlackGrid()), + PhiOdd(Op.FermionRedBlackGrid()) + {}; - virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";} + virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";} + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5. - + RealD scale = std::sqrt(0.5); - + FermionField eta (FermOp.FermionGrid()); FermionField etaOdd (FermOp.FermionRedBlackGrid()); FermionField etaEven(FermOp.FermionRedBlackGrid()); - + gaussian(pRNG,eta); pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Odd,etaOdd,eta); - + FermOp.ImportGauge(U); SchurDifferentiableOperator PCop(FermOp); - - + + PCop.MpcDag(etaOdd,PhiOdd); - + FermOp.MooeeDag(etaEven,PhiEven); - + PhiOdd =PhiOdd*scale; PhiEven=PhiEven*scale; - + }; - + ////////////////////////////////////////////////////// // S = phi^dag (Mdag M)^-1 phi (odd) // + phi^dag (Mdag M)^-1 phi (even) ////////////////////////////////////////////////////// virtual RealD S(const GaugeField &U) { - + FermOp.ImportGauge(U); FermionField X(FermOp.FermionRedBlackGrid()); @@ -137,7 +144,7 @@ class TwoFlavourEvenOddPseudoFermionAction // ////////////////////////////////////////////////////// virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - + std::cout << GridLogDebug << "Calling deriv" << std::endl; FermOp.ImportGauge(U); FermionField X(FermOp.FermionRedBlackGrid()); @@ -151,10 +158,17 @@ class TwoFlavourEvenOddPseudoFermionAction X=zero; DerivativeSolver(Mpc,PhiOdd,X); + std::cout << GridLogDebug << "Calling deriv 2 " << std::endl; Mpc.Mpc(X,Y); - Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; + std::cout << GridLogDebug << "Calling deriv 3 " << std::endl; + Mpc.MpcDeriv(tmp , Y, X ); + std::cout << GridLogDebug << "Calling deriv 4 " << std::endl; + dSdU=tmp; + std::cout << GridLogDebug << "Calling deriv 5 " << std::endl; Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; + std::cout << GridLogDebug << "Calling deriv 6" << std::endl; + // Treat the EE case. (MdagM)^-1 = Minv Minvdag // Deriv defaults to zero. // FermOp.MooeeInvDag(PhiOdd,Y); @@ -165,10 +179,10 @@ class TwoFlavourEvenOddPseudoFermionAction assert(FermOp.ConstEE() == 1); /* - FermOp.MooeeInvDag(PhiOdd,Y); - FermOp.MooeeInv(Y,X); - FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; - FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; + FermOp.MooeeInvDag(PhiOdd,Y); + FermOp.MooeeInv(Y,X); + FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; + FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; */ //dSdU = Ta(dSdU); diff --git a/lib/qcd/hmc/HMCModules.h b/lib/qcd/hmc/HMCModules.h index 6fc48dc2..01549a3c 100644 --- a/lib/qcd/hmc/HMCModules.h +++ b/lib/qcd/hmc/HMCModules.h @@ -72,8 +72,12 @@ public: class GridModule { public: - GridCartesian* get_full() { return grid_.get(); } - GridRedBlackCartesian* get_rb() { return rbgrid_.get(); } + GridCartesian* get_full() { + std::cout << GridLogDebug << "Getting cartesian in module"<< std::endl; + return grid_.get(); } + GridRedBlackCartesian* get_rb() { + std::cout << GridLogDebug << "Getting rb-cartesian in module"<< std::endl; + return rbgrid_.get(); } void set_full(GridCartesian* grid) { grid_.reset(grid); } void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); } diff --git a/lib/qcd/modules/ActionModules.h b/lib/qcd/modules/ActionModules.h index 7b0ca963..cc1382bb 100644 --- a/lib/qcd/modules/ActionModules.h +++ b/lib/qcd/modules/ActionModules.h @@ -40,11 +40,11 @@ namespace Grid { // Actions ////////////////////////////////////////////// -template +template class ActionModuleBase: public HMCModuleBase{ public: - typedef Resource Res; - virtual void acquireResource(Resource& ){}; + typedef R Resource; + virtual void acquireResource(R& ){}; }; @@ -193,57 +193,95 @@ class DBW2GModule: public ActionModule, BetaGaugeActionPar ///////////////////////////////////////// -template -class TwoFlavourFModule: public ActionModule, NoParameters> { - typedef ActionModule, NoParameters> ActionBase; +template class FermionA > +class PseudoFermionModuleBase: public ActionModule, NoParameters> { +protected: + typedef ActionModule, NoParameters> ActionBase; using ActionBase::ActionBase; // for constructors - std::unique_ptr > >solver_mod; - std::unique_ptr> > fop_mod; -public: + typedef std::unique_ptr> > operator_type; + typedef std::unique_ptr > > solver_type; - virtual void acquireResource(typename ActionBase::Res& GridMod){ - fop_mod->AddGridPair(GridMod); + template + void getFermionOperator(Reader& Reader, operator_type &fo, std::string section_name){ + auto &FOFactory = HMC_FermionOperatorModuleFactory::getInstance(); + Reader.push(section_name); + std::string op_name; + read(Reader,"name", op_name); + fo = FOFactory.create(op_name, Reader); + Reader.pop(); } template - TwoFlavourFModule(Reader& Reader) : ActionBase(Reader){ - std::cout << "Constructing TwoFlavourFModule" << std::endl; + void getSolverOperator(Reader& Reader, solver_type &so, std::string section_name){ auto& SolverFactory = HMC_SolverModuleFactory::getInstance(); - Reader.push("Solver"); + Reader.push(section_name); std::string solv_name; read(Reader,"name", solv_name); - solver_mod = SolverFactory.create(solv_name, Reader); - std::cout << "Registered types " << std::endl; - std::cout << SolverFactory.getBuilderList() << std::endl; - solver_mod->print_parameters(); - Reader.pop(); + so = SolverFactory.create(solv_name, Reader); + Reader.pop(); + } +}; +template +class TwoFlavourFModule: public PseudoFermionModuleBase{ + typedef PseudoFermionModuleBase Base; + using Base::Base; - auto &FOFactory = HMC_FermionOperatorModuleFactory::getInstance(); - Reader.push("Operator"); - std::string op_name; - read(Reader,"name", op_name); - fop_mod = FOFactory.create(op_name, Reader); - std::cout << "Registered types " << std::endl; - std::cout << FOFactory.getBuilderList() << std::endl; + typename Base::operator_type fop_mod; + typename Base::solver_type solver_mod; - fop_mod->print_parameters(); - Reader.pop(); + public: + virtual void acquireResource(typename Base::Resource& GridMod){ + fop_mod->AddGridPair(GridMod); + } + // constructor + template + TwoFlavourFModule(Reader& R): PseudoFermionModuleBase(R) { + this->getSolverOperator(R, solver_mod, "Solver"); + this->getFermionOperator(R, fop_mod, "Operator"); + R.pop(); + } - - }; - -private: // acquire resource virtual void initialize() { - this->ActionPtr.reset(new TwoFlavourPseudoFermionAction(*fop_mod->getPtr(), *solver_mod->getPtr(), *solver_mod->getPtr())); + // here temporarily assuming that the force and action solver are the same + this->ActionPtr.reset(new TwoFlavourPseudoFermionAction(*(this->fop_mod->getPtr()), *(this->solver_mod->getPtr()), *(this->solver_mod->getPtr()))); } }; +// very similar, I could have templated this but it is overkilling +template +class TwoFlavourEOFModule: public PseudoFermionModuleBase{ + typedef PseudoFermionModuleBase Base; + using Base::Base; + + typename Base::operator_type fop_mod; + typename Base::solver_type solver_mod; + + public: + virtual void acquireResource(typename Base::Resource& GridMod){ + fop_mod->AddGridPair(GridMod); + } + + // constructor + template + TwoFlavourEOFModule(Reader& R): PseudoFermionModuleBase(R) { + this->getSolverOperator(R, solver_mod, "Solver"); + this->getFermionOperator(R, fop_mod, "Operator"); + R.pop(); + } + + // acquire resource + virtual void initialize() { + // here temporarily assuming that the force and action solver are the same + this->ActionPtr.reset(new TwoFlavourEvenOddPseudoFermionAction(*(this->fop_mod->getPtr()), *(this->solver_mod->getPtr()), *(this->solver_mod->getPtr()))); + } + +}; @@ -308,6 +346,9 @@ static Registrar , HMC_LGTActionModuleFactory > __TwoFlavourFmodXMLInit("TwoFlavours"); +static Registrar , HMC_LGTActionModuleFactory > __TwoFlavourEOFmodXMLInit("TwoFlavoursEvenOdd"); + + // add here the registration for other implementations and readers diff --git a/lib/qcd/modules/FermionOperatorModules.h b/lib/qcd/modules/FermionOperatorModules.h index d48a5757..d0b821cd 100644 --- a/lib/qcd/modules/FermionOperatorModules.h +++ b/lib/qcd/modules/FermionOperatorModules.h @@ -124,8 +124,9 @@ class WilsonFermionModule: public FermionOperatorModuleGridRefs[0].get().get_full()); - this->FOPtr.reset(new WilsonFermion(U, *(this->GridRefs[0].get().get_full()), *(this->GridRefs[0].get().get_rb()), this->Par_.mass)); + auto &GridMod = this->GridRefs[0].get(); + typename FermionImpl::GaugeField U(GridMod.get_full()); + this->FOPtr.reset(new WilsonFermion(U, *(GridMod.get_full()), *(GridMod.get_rb()), this->Par_.mass)); } }; diff --git a/lib/qcd/modules/Modules.h b/lib/qcd/modules/Modules.h index adb88ca1..3d564117 100644 --- a/lib/qcd/modules/Modules.h +++ b/lib/qcd/modules/Modules.h @@ -43,16 +43,6 @@ class NoParameters{}; /* Base class for modules with parameters */ - - - -class ObjectInfo: Serializable { -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectInfo, - std::string, name); -}; - - template < class P > class Parametrized{ public: @@ -101,9 +91,6 @@ class Parametrized{ /* Lowest level abstract module class */ - - - template class HMCModuleBase { public: diff --git a/tests/hmc/Test_hmc_EOWilsonFermionGauge.cc b/tests/hmc/Test_hmc_EOWilsonFermionGauge.cc index e79d7435..26affc69 100644 --- a/tests/hmc/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/hmc/Test_hmc_EOWilsonFermionGauge.cc @@ -1,103 +1,141 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/Test_hmc_EOWilsonFermionGauge.cc +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc - Copyright (C) 2015 +Copyright (C) 2016 -Author: Peter Boyle -Author: Peter Boyle -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 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. +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. +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 */ +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; - -namespace Grid { - namespace QCD { - - -class HmcRunner : public NerscHmcRunner { -public: - - void BuildTheAction (int argc, char **argv) - - { - typedef WilsonImplR ImplPolicy; - typedef WilsonFermionR FermionAction; - typedef typename FermionAction::FermionField FermionField; - - UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - - FGrid = UGrid; - FrbGrid = UrbGrid; - - // temporarily need a gauge field - LatticeGaugeField U(UGrid); - - // Gauge action - WilsonGaugeActionR Waction(5.6); - - Real mass=-0.77; - FermionAction FermOp(U,*FGrid,*FrbGrid,mass); - - ConjugateGradient CG(1.0e-8,10000); - - TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp,CG,CG); - - //Set smearing (true/false), default: false - Nf2.is_smeared=true; - - //Collect actions - ActionLevel Level1(1); - Level1.push_back(&Nf2); - - ActionLevel Level2(4); - Level2.push_back(&Waction); - - TheAction.push_back(Level1); - TheAction.push_back(Level2); - - Run(argc,argv); - }; - -}; - -}} - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); - std::cout< HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonFermionR 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.BuildTheAction(argc,argv); + TheHMC.Resources.LoadBinaryCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.SerialSeed = {1,2,3,4,5}; + RNGpar.ParallelSeed = {6,7,8,9,10}; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + PlaquetteObsParameters PlPar; + PlPar.output_prefix = "Plaquette"; + PlaquetteMod PlaqModule(PlPar); + TheHMC.Resources.AddObservable(&PlaqModule); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // 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.77; + + // Can we define an overloaded operator that does not need U and initialises + // it with zeroes? + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass); + + 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_Factories.cc b/tests/hmc/Test_hmc_Factories.cc index a10c0ba7..da1ab75d 100644 --- a/tests/hmc/Test_hmc_Factories.cc +++ b/tests/hmc/Test_hmc_Factories.cc @@ -68,7 +68,10 @@ int main(int argc, char **argv) { // Construct the module auto myHMCmodule = HMCfactory.create(HMCpar.MD.name, Reader); + myHMCmodule->getPtr()->initialize(Reader); + myHMCmodule->getPtr()->Run(); + Grid_finalize(); /* @@ -100,10 +103,6 @@ int main(int argc, char **argv) { */ - myHMCmodule->getPtr()->initialize(Reader); - myHMCmodule->getPtr()->Run(); - Grid_finalize(); - } // main