portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-03 21:44:33 +00:00
Code Releases Activity

HMC ready but untested

Browse Source
This commit is contained in:
neo
2015-07-04 17:47:50 +09:00
parent 59be55c0ab
commit e6087e1820
6 changed files with 225 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
lib/qcd/QCD.h
View File

@@ -252,47 +252,47 @@ namespace QCD {
//////////////////////////////////////////////////////////////////////////////
//spin
template<class vobj> auto peekSpin(const vobj &rhs,int i) -> decltype(peekIndex<SpinIndex>(rhs,0))
template<class vobj> auto peekSpin(const vobj &rhs,int i) -> decltype(PeekIndex<SpinIndex>(rhs,0))
{
return peekIndex<SpinIndex>(rhs,i);
return PeekIndex<SpinIndex>(rhs,i);
}
template<class vobj> auto peekSpin(const vobj &rhs,int i,int j) -> decltype(peekIndex<SpinIndex>(rhs,0,0))
template<class vobj> auto peekSpin(const vobj &rhs,int i,int j) -> decltype(PeekIndex<SpinIndex>(rhs,0,0))
{
return peekIndex<SpinIndex>(rhs,i,j);
return PeekIndex<SpinIndex>(rhs,i,j);
}
template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<SpinIndex>(rhs,0))
template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<SpinIndex>(rhs,0))
{
return peekIndex<SpinIndex>(rhs,i);
return PeekIndex<SpinIndex>(rhs,i);
}
template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i,int j) -> decltype(peekIndex<SpinIndex>(rhs,0,0))
template<class vobj> auto peekSpin(const Lattice<vobj> &rhs,int i,int j) -> decltype(PeekIndex<SpinIndex>(rhs,0,0))
{
return peekIndex<SpinIndex>(rhs,i,j);
return PeekIndex<SpinIndex>(rhs,i,j);
}
//colour
template<class vobj> auto peekColour(const vobj &rhs,int i) -> decltype(peekIndex<ColourIndex>(rhs,0))
template<class vobj> auto peekColour(const vobj &rhs,int i) -> decltype(PeekIndex<ColourIndex>(rhs,0))
{
return peekIndex<ColourIndex>(rhs,i);
return PeekIndex<ColourIndex>(rhs,i);
}
template<class vobj> auto peekColour(const vobj &rhs,int i,int j) -> decltype(peekIndex<ColourIndex>(rhs,0,0))
template<class vobj> auto peekColour(const vobj &rhs,int i,int j) -> decltype(PeekIndex<ColourIndex>(rhs,0,0))
{
return peekIndex<ColourIndex>(rhs,i,j);
return PeekIndex<ColourIndex>(rhs,i,j);
}
template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<ColourIndex>(rhs,0))
template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<ColourIndex>(rhs,0))
{
return peekIndex<ColourIndex>(rhs,i);
return PeekIndex<ColourIndex>(rhs,i);
}
template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i,int j) -> decltype(peekIndex<ColourIndex>(rhs,0,0))
template<class vobj> auto peekColour(const Lattice<vobj> &rhs,int i,int j) -> decltype(PeekIndex<ColourIndex>(rhs,0,0))
{
return peekIndex<ColourIndex>(rhs,i,j);
return PeekIndex<ColourIndex>(rhs,i,j);
}
//lorentz
template<class vobj> auto peekLorentz(const vobj &rhs,int i) -> decltype(peekIndex<LorentzIndex>(rhs,0))
template<class vobj> auto peekLorentz(const vobj &rhs,int i) -> decltype(PeekIndex<LorentzIndex>(rhs,0))
{
return peekIndex<LorentzIndex>(rhs,i);
return PeekIndex<LorentzIndex>(rhs,i);
}
template<class vobj> auto peekLorentz(const Lattice<vobj> &rhs,int i) -> decltype(peekIndex<LorentzIndex>(rhs,0))
template<class vobj> auto peekLorentz(const Lattice<vobj> &rhs,int i) -> decltype(PeekIndex<LorentzIndex>(rhs,0))
{
return peekIndex<LorentzIndex>(rhs,i);
return PeekIndex<LorentzIndex>(rhs,i);
}
//////////////////////////////////////////////

8
lib/qcd/action/ActionBase.h
View File

@@ -7,10 +7,10 @@ template<class GaugeField>
class Action {
public:
virtual void init(const GaugeField &U) = 0;
virtual RealD S(const GaugeField &U) = 0; // evaluate the action
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void staple(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) = 0;
virtual RealD S(const GaugeField &U) = 0; // evaluate the action
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void staple(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
//virtual void refresh(const GaugeField & ) {} ;
// Boundary conditions?
// Heatbath?

3
lib/qcd/action/gauge/WilsonGaugeAction.h
View File

@@ -14,8 +14,7 @@ namespace Grid{
public:
WilsonGaugeAction(RealD b):beta(b){};
virtual void init(const GaugeField &U) {//FIXME
};
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {};
virtual RealD S(const GaugeField &U) {
return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);

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

@@ -1,6 +1,6 @@
//--------------------------------------------------------------------
/*! @file HMC.h
* @brief Declaration of classes for HybridMonteCarlo update
* @brief Declaration of classes for Hybrid Monte Carlo update
*
* @author Guido Cossu
*/
@@ -30,7 +30,7 @@ namespace Grid{
const HMCparameters Params;
GridSerialRNG sRNG;
GridParallelRNG pRNG;
std::unique_ptr< Integrator<IntegType> > MD;
Integrator<IntegType>& MD;
bool metropolis_test(const RealD DeltaH){
@@ -51,13 +51,13 @@ namespace Grid{
}
}
RealD evolve_step(LatticeColourMatrix& Uin){
MD->init(Uin,pRNG); // set U and initialize P and phi's
RealD evolve_step(LatticeLorentzColourMatrix& U){
MD->init(U,pRNG); // set U and initialize P and phi's
RealD H0 = MD->S(); // current state
std::cout<<"Total H_before = "<< H0 << "\n";
MD->integrate(0);
MD->integrate(U,0);
RealD H1 = MD->S(); // updated state
std::cout<<"Total H_after = "<< H1 << "\n";
@@ -68,7 +68,10 @@ namespace Grid{
public:
HybridMonteCarlo(Integrator<IntegType>& MolDyn, GridBase* grid):MD(&MolDyn),pRNG(grid){
HybridMonteCarlo(HMCparameters Pms,
Integrator<IntegType>& MolDyn,
GridBase* grid):
Params(Pms),MD(MolDyn),pRNG(grid){
//FIXME
// initialize RNGs
@@ -79,7 +82,7 @@ namespace Grid{
void evolve(LatticeColourMatrix& Uin){
void evolve(LatticeLorentzColourMatrix& Uin){
Real DeltaH;
Real timer;
@@ -90,20 +93,21 @@ namespace Grid{
DeltaH = evolve_step(Uin);
std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
std::cout<< "dH = "<< DeltaH << "\n";
std::cout<< " dH = "<< DeltaH << "\n";
// Update matrix
Uin = MD->get_U(); //accept every time
//Uin = MD->get_U(); //accept every time
}
// Actual updates
// 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){
std::cout << "-- # Sweep = "<< iter << "\n";
Ucopy = Uin;
DeltaH = evolve_step(Ucopy);
DeltaH = evolve_step(Uin);
if(metropolis_test(DeltaH)) Uin = MD->get_U();
if(metropolis_test(DeltaH)) Uin = Ucopy;
// need sync?
}

201
lib/qcd/hmc/integrators/Integrator_base.h
View File

@@ -9,6 +9,8 @@
#ifndef INTEGRATOR_INCLUDED
#define INTEGRATOR_INCLUDED
#include <memory>
class Observer;
@@ -22,59 +24,188 @@ namespace Grid{
typedef std::vector<ActionLevel> ActionSet;
typedef std::vector<Observer*> ObserverList;
class Integrator2MN{
const double lambda = 0.1931833275037836;
void step (LatticeColourMatrix&, LatticeColourMatrix&,
int, std::vector<int>&);
class LeapFrog;
struct IntegratorParameters{
int Nexp;
int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize;
IntegratorParameters(int Nexp_,
int MDsteps_,
RealD trajL_):
Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
};
namespace MDutils{
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
}
class IntegratorLeapFrog{
void step (LatticeColourMatrix&, LatticeColourMatrix&,
int, std::vector<int>&);
};
template< class IntegratorPolicy >
class Integrator{
private:
int Nexp;
int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize;
const std::vector<int> Nrel; // relative steps per level
IntegratorParameters Params;
const ActionSet as;
ObserverList observers;
// LatticeColourMatrix* const U; // is shared among all actions - or use a singleton...
LatticeColourMatrix P;
const std::vector<int> Nrel; //relative step size per level
//ObserverList observers; // not yet
std::unique_ptr<LatticeLorentzColourMatrix> P;
IntegratorPolicy TheIntegrator;// contains parameters too
void update_P(int lv,double ep);
void update_U(double ep);
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);
*P -= force*ep;
}
}
void update_U(LatticeLorentzColourMatrix&U, double ep){
//rewrite exponential to deal with the lorentz index?
LatticeColourMatrix Umu(U._grid);
LatticeColourMatrix Pmu(U._grid);
for (int mu = 0; mu < Nd; mu++){
Umu=peekLorentz(U, mu);
Pmu=peekLorentz(*P, mu);
Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
}
}
void register_observers();
void notify_observers();
void integrator_step(int level ,std::vector<Integer>& clock);
friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U,
int level, std::vector<int>& clock,
Integrator<LeapFrog>* Integ);
public:
Integrator(IntegratorParameters Par,
ActionSet& Aset, std::vector<int> Nrel_):
Params(Par),as(Aset),Nrel(Nrel_){
assert(as.size() == Nrel.size());
};
~Integrator(){}
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U,
GridParallelRNG& pRNG){
std::cout<< "Integrator init\n";
if (!P)
P = new LatticeLorentzColourMatrix(U._grid);
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);
}
}
}
public:
Integrator(int Nexp_, int MDsteps_, RealD trajL_,
ActionSet& Aset, ObserverList obs):as(Aset), observers(obs){};
~Integrator(){}
void init(LatticeLorentzColourMatrix&,
GridParallelRNG&);
double S();
void integrate(int level);
LatticeColourMatrix get_U();
RealD S(LatticeLorentzColourMatrix& U){
// Momenta
LatticeComplex Hloc = - trace((*P)*adj(*P));
Complex Hsum = sum(Hloc);
RealD H = Hsum.real();
// 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);
return H;
}
void integrate(LatticeLorentzColourMatrix& U, int level){
std::vector<int> clock;
clock.resize(as.size(),0);
for(int step=0; step< Params.MDsteps; ++step) // MD step
TheIntegrator.step(U,0,clock, *(this));
}
};
namespace MDutils{
void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
}
class MinimumNorm2{
const double lambda = 0.1931833275037836;
public:
void step (LatticeLorentzColourMatrix& U, int level, std::vector<int>& clock);
};
class LeapFrog{
public:
void step (LatticeLorentzColourMatrix& U,
int level, std::vector<int>& clock,
Integrator<LeapFrog>* Integ){
// cl : current level
// fl : final level
// eps : current step size
int fl = Integ->as.size() -1;
double eps = Integ->Params.stepsize;
// Get current level step size
for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
int fin = 1;
for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
fin = 2*Integ->Params.MDsteps*fin - 1;
for(int e=0; e<Integ->Nrel[level]; ++e){
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
if(clock[level] == fin){ // final half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
}
}
}
};
}
}