diff --git a/lib/qcd/hmc/HMCModules.h b/lib/qcd/hmc/HMCModules.h index 768e8a47..0535dd54 100644 --- a/lib/qcd/hmc/HMCModules.h +++ b/lib/qcd/hmc/HMCModules.h @@ -37,25 +37,19 @@ namespace Grid { namespace QCD { //////////////////////////////////////////////////////////////////// -class RNGModuleParameters: Serializable { - +struct RNGModuleParameters: Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(RNGModuleParameters, std::string, serial_seeds, std::string, parallel_seeds,); -public: - std::vector SerialSeed; - std::vector ParallelSeed; - RNGModuleParameters(const std::vector S = std::vector(), - const std::vector P = std::vector()) - : SerialSeed(S), ParallelSeed(P) {} + std::vector getSerialSeeds(){return strToVec(serial_seeds);} + std::vector getParallelSeeds(){return strToVec(parallel_seeds);} + RNGModuleParameters(): serial_seeds("1"), parallel_seeds("1"){} template RNGModuleParameters(Reader& Reader){ read(Reader, "RandomNumberGenerator", *this); - SerialSeed = strToVec(serial_seeds); - ParallelSeed = strToVec(parallel_seeds); } }; @@ -82,12 +76,14 @@ public: GridParallelRNG& get_pRNG() { return *pRNG_.get(); } void seed() { - if (Params_.SerialSeed.size() == 0 && Params_.ParallelSeed.size() == 0) { - std::cout << "Seeds not initialized" << std::endl; + auto SerialSeeds = Params_.getSerialSeeds(); + auto ParallelSeeds = Params_.getParallelSeeds(); + if (SerialSeeds.size() == 0 && ParallelSeeds.size() == 0) { + std::cout << GridLogError << "Seeds not initialized" << std::endl; exit(1); } - sRNG_.SeedFixedIntegers(Params_.SerialSeed); - pRNG_->SeedFixedIntegers(Params_.ParallelSeed); + sRNG_.SeedFixedIntegers(SerialSeeds); + pRNG_->SeedFixedIntegers(ParallelSeeds); } }; diff --git a/lib/qcd/hmc/HMCResourceManager.h b/lib/qcd/hmc/HMCResourceManager.h index aaa49c8c..b98e083a 100644 --- a/lib/qcd/hmc/HMCResourceManager.h +++ b/lib/qcd/hmc/HMCResourceManager.h @@ -75,6 +75,8 @@ class HMCResourceManager { bool have_RNG; bool have_CheckPointer; + // NOTE: operator << is not overloaded for std::vector + // so thsi function is necessary void output_vector_string(const std::vector &vs){ for (auto &i: vs) std::cout << i << " "; @@ -85,13 +87,13 @@ class HMCResourceManager { public: HMCResourceManager() : have_RNG(false), have_CheckPointer(false) {} - template + template void initialize(ReaderClass &Read){ // assumes we are starting from the main node // Geometry GridModuleParameters GridPar(Read); - GridFourDimModule GridMod( GridPar) ; + GridFourDimModule GridMod( GridPar) ; AddGrid("gauge", GridMod); // Checkpointer @@ -100,9 +102,6 @@ class HMCResourceManager { std::string cp_type; read(Read,"name", cp_type); std::cout << "Registered types " << std::endl; - // NOTE: operator << is not overloaded for std::vector - // so it complains here - //std::cout << CPfactory.getBuilderList() << std::endl; output_vector_string(CPfactory.getBuilderList()); @@ -178,7 +177,7 @@ class HMCResourceManager { // Add a named grid set, 4d shortcut void AddFourDimGrid(std::string s) { - GridFourDimModule Mod; + GridFourDimModule Mod; AddGrid(s, Mod); } diff --git a/lib/qcd/hmc/HMC_GridModules.h b/lib/qcd/hmc/HMC_GridModules.h index 6e8d6195..fc4bfba6 100644 --- a/lib/qcd/hmc/HMC_GridModules.h +++ b/lib/qcd/hmc/HMC_GridModules.h @@ -31,41 +31,48 @@ directory #define HMC_GRID_MODULES namespace Grid { -namespace QCD { // Resources // Modules for grids -class GridModuleParameters: Serializable{ - +// Introduce another namespace HMCModules? + +class GridModuleParameters: Serializable{ +public: GRID_SERIALIZABLE_CLASS_MEMBERS(GridModuleParameters, std::string, lattice, - std::string, mpi); + std::string, mpi); -public: - // these namings are ugly - // also ugly the distinction between the serializable members - // and this - std::vector lattice_v; - std::vector mpi_v; + std::vector getLattice(){return strToVec(lattice);} + std::vector getMpi() {return strToVec(mpi);} - GridModuleParameters(const std::vector l_ = std::vector(), - const std::vector mpi_ = std::vector()):lattice_v(l_), mpi_v(mpi_){} - - template - GridModuleParameters(Reader& Reader) { - read(Reader, "LatticeGrid", *this); - lattice_v = strToVec(lattice); - mpi_v = strToVec(mpi); - if (mpi_v.size() != lattice_v.size()) { - std::cout << "Error in GridModuleParameters: lattice and mpi dimensions " + void check(){ + if (getLattice().size() != getMpi().size()) { + std::cout << GridLogError + << "Error in GridModuleParameters: lattice and mpi dimensions " "do not match" << std::endl; exit(1); } + } + + template + GridModuleParameters(Reader& Reader) { + read(Reader, name, *this); + check(); } + + // Save on file + template< class WriterClass> + void save(Writer& Writer){ + check(); + write(Writer, name, *this); + } +private: + std::string name = "LatticeGrid"; }; +// Lower level class class GridModule { public: GridCartesian* get_full() { @@ -84,27 +91,33 @@ class GridModule { }; -// helpers -// FIXME define a class accepting also real vtypes +//////////////////////////////////// +// Classes for the user +//////////////////////////////////// +// Note: the space time grid must be out of the QCD namespace +template< class vector_type> class GridFourDimModule : public GridModule { public: - // add a function to create the module from a Reader GridFourDimModule() { + using namespace QCD; set_full(SpaceTimeGrid::makeFourDimGrid( - GridDefaultLatt(), GridDefaultSimd(4, vComplex::Nsimd()), + GridDefaultLatt(), GridDefaultSimd(4, vector_type::Nsimd()), GridDefaultMpi())); set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get())); } - template GridFourDimModule(GridModuleParameters Params) { - if (Params.lattice_v.size() == 4) { + using namespace QCD; + Params.check(); + std::vector lattice_v = Params.getLattice(); + std::vector mpi_v = Params.getMpi(); + if (lattice_v.size() == 4) { set_full(SpaceTimeGrid::makeFourDimGrid( - Params.lattice_v, GridDefaultSimd(4, vector_type::Nsimd()), - Params.mpi_v)); + lattice_v, GridDefaultSimd(4, vector_type::Nsimd()), + mpi_v)); set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get())); } else { - std::cout + std::cout << GridLogError << "Error in GridFourDimModule: lattice dimension different from 4" << std::endl; exit(1); @@ -112,10 +125,9 @@ class GridFourDimModule : public GridModule { } }; +typedef GridFourDimModule GridDefaultFourDimModule; - -} // namespace QCD } // namespace Grid #endif // HMC_GRID_MODULES \ No newline at end of file diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index 8033bd41..6acb872c 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -131,6 +131,49 @@ class Integrator { as[level].apply(update_P_hireps, Representations, Mom, U, ep); } + void implicit_update_P(MomentaField& Mom, Field& U, int level, double ep) { + // Fundamental updates, include smearing + MomentaField Msum(Mom._grid); + Msum = zero; + for (int a = 0; a < as[level].actions.size(); ++a) { + // Compute the force + // We need to compute the derivative of the actions + // only once + Field force(U._grid); + conformable(U._grid, Mom._grid); + Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); + as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta + + std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl; + if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); + force = FieldImplementation::projectForce(force); // Ta for gauge fields + Real force_abs = std::sqrt(norm2(force)/U._grid->gSites()); + std::cout << GridLogIntegrator << "Force average: " << force_abs << std::endl; + Msum += force; + } + + MomentaField NewMom = Mom; + MomentaField OldMom = Mom; + double threshold = 1e-6; + // Here run recursively + do{ + MomentaField MomDer(Mom._grid); + OldMom = NewMom; + // Compute the derivative of the kinetic term + // with respect to the gauge field + + // Laplacian.Mder(NewMom, MomDer); + // NewMom = Mom - ep*(MomDer + Msum); + + } while (norm2(NewMom - OldMom) > threshold); + + Mom = NewMom; + + + // update the auxiliary fields momenta + } + + void update_U(Field& U, double ep) { update_U(P, U, ep); diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h index f9709bec..d3477077 100644 --- a/lib/qcd/hmc/integrators/Integrator_algorithm.h +++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h @@ -285,6 +285,74 @@ class ForceGradient : public Integrator > +class ImplicitLeapFrog : public Integrator { + public: + typedef ImplicitLeapFrog + Algorithm; + INHERIT_FIELD_TYPES(FieldImplementation); + + // Riemannian manifold metric operator + // Hermitian operator Fisher + + std::string integrator_name(){return "ImplicitLeapFrog";} + + ImplicitLeapFrog(GridBase* grid, IntegratorParameters Par, + ActionSet& Aset, SmearingPolicy& Sm) + : Integrator( + grid, Par, Aset, Sm){}; + + void step(Field& U, int level, int _first, int _last) { + int fl = this->as.size() - 1; + // level : current level + // fl : final level + // eps : current step size + + // Get current level step size + RealD eps = this->Params.trajL/this->Params.MDsteps; + for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier; + + int multiplier = this->as[level].multiplier; + for (int e = 0; e < multiplier; ++e) { + int first_step = _first && (e == 0); + int last_step = _last && (e == multiplier - 1); + + if (first_step) { // initial half step + this->implicit_update_P(U, level, eps / 2.0); + } + + if (level == fl) { // lowest level + this->update_U(U, eps); + } else { // recursive function call + this->step(U, level + 1, first_step, last_step); + } + + int mm = last_step ? 1 : 2; + this->update_P(U, level, mm * eps / 2.0); + } + } +}; + + + + + + + + } } diff --git a/lib/qcd/utils/CovariantLaplacian.h b/lib/qcd/utils/CovariantLaplacian.h index 1b40d32c..b8c264a7 100644 --- a/lib/qcd/utils/CovariantLaplacian.h +++ b/lib/qcd/utils/CovariantLaplacian.h @@ -33,13 +33,27 @@ directory namespace Grid { namespace QCD { +//////////////////////////////////////////////////////////// +// Laplacian operator L on adjoint fields +// +// phi: adjoint field +// L: D_mu^dag D_mu +// +// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag + +// U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu) +// -2phi(x)] +// +// Operator designed to be encapsulated by +// an HermitianLinearOperator<.. , ..> +//////////////////////////////////////////////////////////// + template class LaplacianAdjointField { public: INHERIT_GIMPL_TYPES(Impl); - typedef SU::LatticeAlgebraVector AVector; - - LaplacianAdjointField(GridBase* grid) : U(Nd, grid){}; + + LaplacianAdjointField(GridBase* grid, const RealD k = 1.0) : + U(Nd, grid), kappa(k){}; void ImportGauge(const GaugeField& _U) { for (int mu = 0; mu < Nd; mu++) { @@ -49,61 +63,48 @@ class LaplacianAdjointField { void Mdiag(const GaugeLinkField& in, GaugeLinkField& out) { assert(0); } - void Mdir(const GaugeLinkField& in, GaugeLinkField& out, int dir, int disp) { assert(0); } - -/* - // Operator with algebra vector inputs and outputs - void M2(const AVector& in, AVector& out) { - double kappa = 0.9; - //Reconstruct matrix - - GaugeLinkField tmp(in._grid); - GaugeLinkField tmp2(in._grid); - GaugeLinkField sum(in._grid); - GaugeLinkField out_mat(in._grid); - GaugeLinkField in_mat(in._grid); - SU::FundamentalLieAlgebraMatrix(in, in_mat); - - sum = zero; - for (int mu = 0; mu < Nd; mu++) { - tmp = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]); - tmp2 = adj(U[mu]) * in_mat * U[mu]; - sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_mat; - } - out_mat = (1.0 - kappa) * in_mat - kappa/(double(4*Nd)) * sum; - // Project - SU::projectOnAlgebra(out, out_mat); + void Mdir(const GaugeLinkField& in, GaugeLinkField& out, int dir, int disp) { + assert(0); } -*/ - - void M(const GaugeLinkField& in, GaugeLinkField& out) { - double kappa = 0.999; - //Reconstruct matrix + void M(const GaugeLinkField& in, GaugeLinkField& out) { GaugeLinkField tmp(in._grid); GaugeLinkField tmp2(in._grid); GaugeLinkField sum(in._grid); sum = zero; for (int mu = 0; mu < Nd; mu++) { - tmp = U[mu] * Cshift(in, mu, +1) * adj(U[mu]); + tmp = U[mu] * Cshift(in, mu, +1) * adj(U[mu]); tmp2 = adj(U[mu]) * in * U[mu]; sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in; } - out = (1.0 - kappa) * in - kappa/(double(4*Nd)) * sum; + out = (1.0 - kappa) * in - kappa / (double(4 * Nd)) * sum; + } + + void MDeriv(const GaugeLinkField& in, GaugeLinkField& out, bool dag){ + RealD factor = - kappa / (double(4 * Nd)) + if (!dag) + out = factor * Cshift(in, mu, +1) * adj(U[mu]) + adj(U[mu]) * in; + else + out = factor * U[mu] * Cshift(in, mu, +1) + in * U[mu]; } private: + RealD kappa; std::vector U; }; +// This is just a debug tests +// not meant to be used + template class LaplacianAlgebraField { public: INHERIT_GIMPL_TYPES(Impl); typedef SU::LatticeAlgebraVector AVector; - LaplacianAlgebraField(GridBase* grid) : U(Nd, grid){}; + LaplacianAlgebraField(GridBase* grid, const RealD k) : + U(Nd, grid), kappa(k){}; void ImportGauge(const GaugeField& _U) { for (int mu = 0; mu < Nd; mu++) { @@ -117,28 +118,28 @@ class LaplacianAlgebraField { // Operator with algebra vector inputs and outputs void M(const AVector& in, AVector& out) { - double kappa = 0.999; - //Reconstruct matrix - GaugeLinkField tmp(in._grid); GaugeLinkField tmp2(in._grid); GaugeLinkField sum(in._grid); GaugeLinkField out_mat(in._grid); GaugeLinkField in_mat(in._grid); + + // Reconstruct matrix SU::FundamentalLieAlgebraMatrix(in, in_mat); sum = zero; for (int mu = 0; mu < Nd; mu++) { - tmp = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]); + tmp = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]); tmp2 = adj(U[mu]) * in_mat * U[mu]; sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_mat; } - out_mat = (1.0 - kappa) * in_mat - kappa/(double(4*Nd)) * sum; + out_mat = (1.0 - kappa) * in_mat - kappa / (double(4 * Nd)) * sum; // Project SU::projectOnAlgebra(out, out_mat); } private: + RealD kappa; std::vector U; }; diff --git a/tests/hmc/Test_hmc_WilsonGauge.cc b/tests/hmc/Test_hmc_WilsonGauge.cc index 56c72bd5..84240ff1 100644 --- a/tests/hmc/Test_hmc_WilsonGauge.cc +++ b/tests/hmc/Test_hmc_WilsonGauge.cc @@ -59,8 +59,8 @@ int main(int argc, char **argv) { TheHMC.Resources.LoadBinaryCheckpointer(CPparams); RNGModuleParameters RNGpar; - RNGpar.SerialSeed = {1,2,3,4,5}; - RNGpar.ParallelSeed = {6,7,8,9,10}; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; TheHMC.Resources.SetRNGSeeds(RNGpar); // Construct observables diff --git a/tests/hmc/Test_hmc_WilsonTMFermionGauge.cc b/tests/hmc/Test_hmc_WilsonTMFermionGauge.cc new file mode 100644 index 00000000..abb77867 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonTMFermionGauge.cc @@ -0,0 +1,147 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2015 + +Author: Peter Boyle +Author: neo +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 WilsonTMFermionR 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.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 = 3.9 ; + SymanzikGaugeActionR Waction(beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + Real mass = -0.89163; + Real mu = 0.01; + + // Can we define an overloaded operator that does not need U and initialises + // it with zeroes? + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, mu); + + ConjugateGradient CG(1.0e-8, 2000); + + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // With modules + /* + + TwoFlavourFmodule TwoFMod(Reader); + + */ + + // 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/solver/Test_laplacian.cc b/tests/solver/Test_laplacian.cc index 6660a4ec..2ba14fc2 100644 --- a/tests/solver/Test_laplacian.cc +++ b/tests/solver/Test_laplacian.cc @@ -48,28 +48,33 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(&Grid); SU::HotConfiguration(pRNG,Umu); + double Kappa = 0.9999; + typedef SU::LatticeAlgebraVector AVector; // Source and result in the algebra AVector src_vec(&Grid); random(pRNG, src_vec); AVector result_vec(&Grid); result_vec = zero; + LatticeColourMatrix src(&Grid); SU::FundamentalLieAlgebraMatrix(src_vec, src); LatticeColourMatrix result(&Grid); result=zero; - LaplacianAdjointField Laplacian(&Grid); + LaplacianAdjointField Laplacian(&Grid, Kappa); Laplacian.ImportGauge(Umu); HermitianLinearOperator,LatticeColourMatrix> HermOp(Laplacian); ConjugateGradient CG(1.0e-8,10000); + std::cout << GridLogMessage << "Testing the Laplacian using the full matrix" < LaplacianAlgebra(&Grid); + LaplacianAlgebraField LaplacianAlgebra(&Grid, Kappa); LaplacianAlgebra.ImportGauge(Umu); HermitianLinearOperator,AVector> HermOpAlg(LaplacianAlgebra); ConjugateGradient CG_Algebra(1.0e-8,10000); + std::cout << GridLogMessage << "Testing the Laplacian using the algebra vectors" <