portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Debugged vector version of ProjectOnGroup

Browse Source
This commit is contained in:
neo 2015-07-06 02:24:58 +09:00
parent b1f94fa292
commit 62d8952c0a
10 changed files with 89 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -10,7 +10,6 @@ class Action {
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?

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

@ -17,17 +17,28 @@ namespace Grid{
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {};
virtual RealD S(const GaugeField &U) {
return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
std::cout << "Plaq : "<<plaq << "\n";
double vol = U._grid->gSites();
return beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
};
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
//FIXME loop on directions
//not optimal implementation FIXME
//extend Ta to include Lorentz indexes
RealD factor = 0.5*beta/RealD(Nc);
std::cout << "Debug force Wilson "<< factor << "\n";
MatrixField dSdU_mu(U._grid);
WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,0);
};
virtual void staple(const GaugeField &stap,GaugeField & U) {
//FIXME loop on directions
MatrixField stap_mu(U._grid);
WilsonLoops<MatrixField,GaugeField>::Staple(stap_mu,U,0);
MatrixField Umu(U._grid);
for (int mu=0; mu < Nd; mu++){
Umu = PeekIndex<LorentzIndex>(U,mu);
// Staple in direction mu
WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,mu);
std::cout << "Staple norm ["<<mu<<"] = "<< norm2(dSdU_mu) <<"\n";
dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor;
std::cout << "Force norm ["<<mu<<"] = "<< norm2(dSdU_mu) << " : Umu "<< norm2(Umu)<<"\n";
pokeLorentz(dSdU, dSdU_mu, mu);
}
};
};

2
lib/qcd/hmc/HMC.cc
View File

@ -9,7 +9,7 @@ namespace Grid{
////////////////////////////// Default values
Nsweeps = 10;
TotalSweeps = 10;
ThermalizationSteps = 0;
ThermalizationSteps = 1;
StartingConfig = 0;
SaveInterval = 1;
Filename_prefix = "Conf_";

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

@ -53,13 +53,13 @@ namespace Grid{
RealD evolve_step(LatticeLorentzColourMatrix& U){
MD->init(U,pRNG); // set U and initialize P and phi's
RealD H0 = MD->S(); // current state
MD.init(U,pRNG); // set U and initialize P and phi's
RealD H0 = MD.S(U); // current state
std::cout<<"Total H_before = "<< H0 << "\n";
MD->integrate(U,0);
MD.integrate(U,0);
RealD H1 = MD->S(); // updated state
RealD H1 = MD.S(U); // updated state
std::cout<<"Total H_after = "<< H1 << "\n";
return (H1-H0);
@ -104,10 +104,11 @@ namespace Grid{
for(int iter=Params.StartingConfig;
iter < Params.Nsweeps+Params.StartingConfig; ++iter){
std::cout << "-- # Sweep = "<< iter << "\n";
Ucopy = Uin;
DeltaH = evolve_step(Ucopy);
if(metropolis_test(DeltaH)) Uin = Ucopy;
if(metropolis_test(DeltaH)) Uin = Ucopy;
// need sync?
}

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

@ -25,7 +25,7 @@ namespace Grid{
typedef std::vector<Observer*> ObserverList;
class LeapFrog;
struct IntegratorParameters{
int Nexp;
@ -77,7 +77,9 @@ namespace Grid{
for (int mu = 0; mu < Nd; mu++){
Umu=peekLorentz(U, mu);
Pmu=peekLorentz(*P, mu);
std::cout << "U norm ["<<mu<<"] = "<< norm2(Umu) << " : Pmu "<< norm2(Pmu)<<"\n";
Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
pokeLorentz(U, Umu, mu);
}
}
@ -102,10 +104,15 @@ namespace Grid{
void init(LatticeLorentzColourMatrix& U,
GridParallelRNG& pRNG){
std::cout<< "Integrator init\n";
if (!P)
P = new LatticeLorentzColourMatrix(U._grid);
if (!P){
std::unique_ptr<LatticeLorentzColourMatrix> Pnew(new LatticeLorentzColourMatrix(U._grid));
P = std::move(Pnew);
}
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);
@ -116,12 +123,19 @@ namespace Grid{
RealD S(LatticeLorentzColourMatrix& U){
LatticeComplex Hloc(U._grid);
Hloc = zero;
// Momenta
LatticeComplex Hloc = - trace((*P)*adj(*P));
for (int mu=0; mu <Nd; mu++){
LatticeColourMatrix Pmu = peekLorentz(*P, mu);
Hloc -= trace(Pmu*adj(Pmu));
}
Complex Hsum = sum(Hloc);
RealD H = Hsum.real();
std::cout << "H_p = "<< H << "\n";
// Actions
for(int level=0; level<as.size(); ++level)
for(int actionID=0; actionID<as.at(level).size(); ++actionID)
@ -136,7 +150,7 @@ namespace Grid{
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));
TheIntegrator.step(U,0,clock, (this));
}
};

5
lib/simd/Grid_vector_unops.h
View File

@ -94,6 +94,11 @@ namespace Grid {
inline Grid_simd<S,V> rsqrt(const Grid_simd<S,V> &r) {
return SimdApply(RSqrtRealFunctor<S>(),r);
}
template < class Scalar >
inline Scalar rsqrt(const Scalar &r) {
return (RSqrtRealFunctor<Scalar>(),r);
}
template < class S, class V >
inline Grid_simd<S,V> cos(const Grid_simd<S,V> &r) {
return SimdApply(CosRealFunctor<S>(),r);

6
lib/tensors/Tensor_Ta.h
View File

@ -62,14 +62,16 @@ namespace Grid {
{
// need a check for the group type?
iMatrix<vtype,N> ret(arg);
RealD nrm;
vtype nrm;
vtype inner;
for(int c1=0;c1<N;c1++){
zeroit(inner);
for(int c2=0;c2<N;c2++)
inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
nrm = 1.0/sqrt(Reduce(toReal(inner)));
//nrm = 1.0/sqrt(Reduce(toReal(inner)));
nrm = rsqrt(inner);
for(int c2=0;c2<N;c2++)
ret._internal[c1][c2]*= nrm;

2
lib/tensors/Tensor_exp.h
View File

@ -28,7 +28,9 @@ namespace Grid {
temp *= alpha/ComplexD(i);
temp = unit + temp*arg;
}
return ProjectOnGroup(temp);//maybe not strictly necessary
}

14
tests/Test_hmc.cc
View File

@ -22,7 +22,13 @@ int main (int argc, char ** argv)
double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridParallelRNG pRNG(&Fine);
pRNG.SeedRandomDevice();
LatticeLorentzColourMatrix U(&Fine);
SU3::ColdConfiguration(pRNG, U);
// simplify template?
WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(6.0);
@ -33,12 +39,14 @@ int main (int argc, char ** argv)
FullSet.push_back(Level1);
// Create integrator
IntegratorParameters MDpar(12,10,1.0);
IntegratorParameters MDpar(12,50,1.0);
std::vector<int> rel ={1};
Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);
// Create HMC
HMCparameters HMCpar;
HybridMonteCarlo<LeapFrog> HMCrun(HMCpar, MDleapfrog, &Fine);
HybridMonteCarlo<LeapFrog> HMC(HMCpar, MDleapfrog, &Fine);
HMC.evolve(U);
}

23
tests/Test_main.cc
View File

@ -220,15 +220,19 @@ int main (int argc, char ** argv)
std::cout << cm << std::endl;
cm = Exponentiate(cm, 1.0, 12);
cm = Exponentiate(cm, 2.0, 12);
std::cout << cm << " " << std::endl;
Complex det = Determinant(cm);
std::cout << "determinant: " << det << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
cm = ProjectOnGroup(cm);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
cm = ProjectOnGroup(cm);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
// det = Determinant(cm);
// std::cout << "determinant: " << det << std::endl;
@ -245,7 +249,24 @@ int main (int argc, char ** argv)
trscMat = trace(scMat); // Trace
// Exponentiate test
std::vector<int> mysite {0,0,0,0};
random(FineRNG,cMat);
cMat = Ta(cMat);
peekSite(cm, cMat, mysite);
std::cout << cm << " " << std::endl;
cm = Exponentiate(cm, 1.0, 12);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout << "norm cMmat : " << norm2(cMat) << std::endl;
cMat = expMat(cMat, ComplexD(1.0, 0.0));
std::cout << "norm expMat: " << norm2(cMat) << std::endl;
peekSite(cm, cMat, mysite);
std::cout << cm << " " << std::endl;
std::cout << "determinant: " << Determinant(cm) << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
// LatticeComplex trlcMat(&Fine);
// trlcMat = trace(lcMat); // Trace involving iVector - now generates error