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

Debugged vector version of ProjectOnGroup

Browse Source
This commit is contained in:
neo
2015-07-06 02:24:58 +09:00
parent e6087e1820
commit 0ffcdf6204
10 changed files with 89 additions and 26 deletions
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
}