portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Merge branch 'master' of https://github.com/paboyle/Grid

Browse Source 
Conflicts:
	lib/Make.inc
	lib/qcd/hmc/HMC.h
	tests/Make.inc
	tests/Test_hmc_WilsonFermionGauge.cc
This commit is contained in:
Peter Boyle
2015-08-01 22:24:54 +09:00
parent f37552a385 3dd846c93c
commit 35feb93f56
16 changed files with 8261 additions and 155 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
lib/qcd/hmc/HMC.h
View File

@ -1,15 +1,16 @@
//--------------------------------------------------------------------
/*! @file HMC.h
* @brief Declaration of classes for Hybrid Monte Carlo update
* @brief Classes for Hybrid Monte Carlo update
*
* @author Guido Cossu
* Time-stamp: <2015-07-30 16:58:26 neo>
*/
//--------------------------------------------------------------------
#ifndef HMC_INCLUDED
#define HMC_INCLUDED
#include <string>
#include <memory>
namespace Grid{
namespace QCD{
@ -20,7 +21,7 @@ namespace Grid{
Integer ThermalizationSteps;
Integer StartingConfig;
Integer SaveInterval; //Setting to 0 does not save configurations
std::string Filename_prefix; // To save configurations
std::string Filename_prefix; // To save configurations and rng seed
HMCparameters();
};
@ -46,8 +47,8 @@ namespace Grid{
random(sRNG,rn_test);
std::cout<<GridLogMessage<< "--------------------------------------------\n";
std::cout<<GridLogMessage<< "dH = "<<DeltaH << " Random = "<< rn_test
<< "\nAcc. Probability = " << ((prob<1.0)? prob: 1.0)<< " ";
std::cout<<GridLogMessage<< "dH = "<<DeltaH << " Random = "<< rn_test <<"\n";
std::cout<<GridLogMessage<< "Acc. Probability = " << ((prob<1.0)? prob: 1.0)<< " ";
if((prob >1.0) || (rn_test <= prob)){ // accepted
std::cout<<GridLogMessage <<"-- ACCEPTED\n";
@ -62,8 +63,7 @@ namespace Grid{
/////////////////////////////////////////////////////////
// Evolution
/////////////////////////////////////////////////////////
RealD evolve_step(LatticeLorentzColourMatrix& U){
RealD evolve_step(LatticeGaugeField& U){
MD.init(U); // set U and initialize P and phi's
RealD H0 = MD.S(U); // initial state action
@ -73,8 +73,6 @@ namespace Grid{
RealD H1 = MD.S(U); // updated state action
std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
std::cout<<GridLogMessage<<"DeltaH is "<< H1-H0 << "\n";
return (H1-H0);
}
@ -89,12 +87,10 @@ namespace Grid{
sRNG.SeedRandomDevice();
}
~HybridMonteCarlo(){};
void evolve(LatticeLorentzColourMatrix& Uin){
void evolve(LatticeGaugeField& Uin){
Real DeltaH;
// Thermalizations
@ -102,13 +98,13 @@ namespace Grid{
std::cout<<GridLogMessage << "-- # Thermalization step = "<< iter << "\n";
DeltaH = evolve_step(Uin);
std::cout<<GridLogMessage<< " dH = "<< DeltaH << "\n";
std::cout<<GridLogMessage<< "dH = "<< DeltaH << "\n";
}
// Actual updates (evolve a copy Ucopy then copy back eventually)
LatticeLorentzColourMatrix Ucopy(Uin._grid);
for(int iter=Params.StartingConfig; iter < Params.Nsweeps+Params.StartingConfig; ++iter){
LatticeGaugeField Ucopy(Uin._grid);
for(int iter=Params.StartingConfig;
iter < Params.Nsweeps+Params.StartingConfig; ++iter){
std::cout<<GridLogMessage << "-- # Sweep = "<< iter << "\n";
Ucopy = Uin;
@ -117,8 +113,6 @@ namespace Grid{
if(metropolis_test(DeltaH)) Uin = Ucopy;
//need sync? // Query Guido on this comment
}
}
};

55
lib/qcd/hmc/integrators/Integrator.h
View File

@ -1,8 +1,9 @@
//--------------------------------------------------------------------
/*! @file Integrator.h
* @brief Declaration of classes for the Molecular Dynamics integrator
* @brief Classes for the Molecular Dynamics integrator
*
* @author Guido Cossu
* Time-stamp: <2015-07-30 16:21:29 neo>
*/
//--------------------------------------------------------------------
@ -16,8 +17,16 @@ class Observer;
namespace Grid{
namespace QCD{
typedef Action<LatticeLorentzColourMatrix>* ActPtr; // now force the same colours as the rest of the code
typedef std::vector<ActPtr> ActionLevel;
typedef Action<LatticeGaugeField>* ActPtr; // now force the same colours as the rest of the code
struct ActionLevel{
int multiplier;
public:
std::vector<ActPtr> actions;
explicit ActionLevel(int mul = 1):multiplier(mul){assert (mul > 0);};
void push_back(ActPtr ptr){
actions.push_back(ptr);
}
};
typedef std::vector<ActionLevel> ActionSet;
typedef std::vector<Observer*> ObserverList;
@ -35,8 +44,8 @@ namespace Grid{
namespace MDutils{
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
void generate_momenta(LatticeGaugeField&,GridParallelRNG&);
void generate_momenta_su3(LatticeGaugeField&,GridParallelRNG&);
}
/*! @brief Class for Molecular Dynamics management */
@ -45,8 +54,7 @@ namespace Grid{
private:
IntegratorParameters Params;
const ActionSet as;
const std::vector<int> Nrel; //relative step size per level
std::unique_ptr<LatticeLorentzColourMatrix> P;
std::unique_ptr<LatticeGaugeField> P;
GridParallelRNG pRNG;
//ObserverList observers; // not yet
@ -55,15 +63,15 @@ namespace Grid{
void register_observers();
void notify_observers();
void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
for(int a=0; a<as[level].size(); ++a){
LatticeLorentzColourMatrix force(U._grid);
as[level].at(a)->deriv(U,force);
void update_P(LatticeGaugeField&U, int level,double ep){
for(int a=0; a<as[level].actions.size(); ++a){
LatticeGaugeField force(U._grid);
as[level].actions.at(a)->deriv(U,force);
*P -= force*ep;
}
}
void update_U(LatticeLorentzColourMatrix&U, double ep){
void update_U(LatticeGaugeField&U, double ep){
//rewrite exponential to deal automatically with the lorentz index?
LatticeColourMatrix Umu(U._grid);
LatticeColourMatrix Pmu(U._grid);
@ -77,34 +85,31 @@ namespace Grid{
}
friend void IntegratorAlgorithm::step (LatticeLorentzColourMatrix& U,
friend void IntegratorAlgorithm::step (LatticeGaugeField& U,
int level, std::vector<int>& clock,
Integrator<IntegratorAlgorithm>* Integ);
public:
Integrator(GridBase* grid, IntegratorParameters Par,
ActionSet& Aset, std::vector<int> Nrel_):
Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
assert(as.size() == Nrel.size());
ActionSet& Aset):
Params(Par),as(Aset),P(new LatticeGaugeField(grid)),pRNG(grid){
pRNG.SeedRandomDevice();
};
~Integrator(){}
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U){
void init(LatticeGaugeField& U){
std::cout<<GridLogMessage<< "Integrator init\n";
MDutils::generate_momenta(*P,pRNG);
for(int level=0; level< as.size(); ++level){
for(int actionID=0; actionID<as.at(level).size(); ++actionID){
as[level].at(actionID)->init(U, pRNG);
for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
as[level].actions.at(actionID)->init(U, pRNG);
}
}
}
// Calculate action
RealD S(LatticeLorentzColourMatrix& U){
RealD S(LatticeGaugeField& U){
LatticeComplex Hloc(U._grid);
Hloc = zero;
// Momenta
@ -120,15 +125,15 @@ namespace Grid{
// Actions
for(int level=0; level<as.size(); ++level)
for(int actionID=0; actionID<as.at(level).size(); ++actionID)
H += as[level].at(actionID)->S(U);
for(int actionID=0; actionID<as[level].actions.size(); ++actionID)
H += as[level].actions.at(actionID)->S(U);
std::cout<<GridLogMessage << "Total action H = "<< H << "\n";
return H;
}
void integrate(LatticeLorentzColourMatrix& U){
void integrate(LatticeGaugeField& U){
std::vector<int> clock;
clock.resize(as.size(),0);
for(int step=0; step< Params.MDsteps; ++step) // MD step

17
lib/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -27,14 +27,13 @@ namespace Grid{
int fl = Integ->as.size() -1;
double eps = Integ->Params.stepsize;
for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
for(int l=0; l<=level; ++l) eps/= 2.0*Integ->as[l].multiplier;
int fin = Integ->Nrel[0];
for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
int fin = Integ->as[0].multiplier;
for(int l=1; l<=level; ++l) fin*= 2.0*Integ->as[l].multiplier;
fin = 3*Integ->Params.MDsteps*fin -1;
for(int e=0; e<Integ->Nrel[level]; ++e){
for(int e=0; e<Integ->as[level].multiplier; ++e){
if(clock[level] == 0){ // initial half step
Integ->update_P(U,level,lambda*eps);
++clock[level];
@ -98,13 +97,13 @@ namespace Grid{
double eps = Integ->Params.stepsize;
// Get current level step size
for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
for(int l=0; l<=level; ++l) eps/= Integ->as[l].multiplier;
int fin = 1;
for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
for(int l=0; l<=level; ++l) fin*= Integ->as[l].multiplier;
fin = 2*Integ->Params.MDsteps*fin - 1;
for(int e=0; e<Integ->Nrel[level]; ++e){
for(int e=0; e<Integ->as[level].multiplier; ++e){
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2.0);