mirror of
https://github.com/paboyle/Grid.git
synced 2025-11-04 05:54:32 +00:00
HMC for Wilson Gauge action works
Fixed bug in momenta generation
This commit is contained in:
neo
@@ -27,16 +27,13 @@ namespace Grid{
|
||||
//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);
|
||||
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);
|
||||
}
|
||||
};
|
||||
|
||||
@@ -7,9 +7,9 @@ namespace Grid{
|
||||
// FIXME fill this constructor now just default values
|
||||
|
||||
////////////////////////////// Default values
|
||||
Nsweeps = 10;
|
||||
TotalSweeps = 10;
|
||||
ThermalizationSteps = 1;
|
||||
Nsweeps = 100;
|
||||
TotalSweeps = 20;
|
||||
ThermalizationSteps = 20;
|
||||
StartingConfig = 0;
|
||||
SaveInterval = 1;
|
||||
Filename_prefix = "Conf_";
|
||||
|
||||
@@ -10,6 +10,7 @@
|
||||
|
||||
#include <string>
|
||||
#include <memory>
|
||||
#include <unistd.h>
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
@@ -55,12 +56,12 @@ namespace Grid{
|
||||
|
||||
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";
|
||||
std::cout<<"Total H before = "<< H0 << "\n";
|
||||
|
||||
MD.integrate(U,0);
|
||||
|
||||
RealD H1 = MD.S(U); // updated state
|
||||
std::cout<<"Total H_after = "<< H1 << "\n";
|
||||
std::cout<<"Total H after = "<< H1 << "\n";
|
||||
|
||||
return (H1-H0);
|
||||
}
|
||||
@@ -84,19 +85,14 @@ namespace Grid{
|
||||
|
||||
void evolve(LatticeLorentzColourMatrix& Uin){
|
||||
Real DeltaH;
|
||||
Real timer;
|
||||
|
||||
// Thermalizations
|
||||
for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
|
||||
std::cout << "-- # Thermalization step = "<< iter << "\n";
|
||||
|
||||
DeltaH = evolve_step(Uin);
|
||||
|
||||
std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
|
||||
std::cout<< " dH = "<< DeltaH << "\n";
|
||||
|
||||
// Update matrix
|
||||
//Uin = MD->get_U(); //accept every time
|
||||
}
|
||||
|
||||
// Actual updates (evolve a copy Ucopy then copy back eventually)
|
||||
@@ -111,11 +107,7 @@ namespace Grid{
|
||||
if(metropolis_test(DeltaH)) Uin = Ucopy;
|
||||
//need sync?
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
}// QCD
|
||||
|
||||
@@ -15,10 +15,12 @@ namespace Grid{
|
||||
|
||||
void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
|
||||
LatticeColourMatrix Pmu(P._grid);
|
||||
Pmu = zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
|
||||
pokeLorentz(P, Pmu, mu);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -71,14 +71,13 @@ namespace Grid{
|
||||
|
||||
|
||||
void update_U(LatticeLorentzColourMatrix&U, double ep){
|
||||
//rewrite exponential to deal with the lorentz index?
|
||||
//rewrite exponential to deal automatically 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);
|
||||
std::cout << "U norm ["<<mu<<"] = "<< norm2(Umu) << " : Pmu "<< norm2(Pmu)<<"\n";
|
||||
Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
|
||||
Umu = expMat(Pmu, ep)*Umu;
|
||||
pokeLorentz(U, Umu, mu);
|
||||
}
|
||||
|
||||
@@ -107,12 +106,8 @@ namespace 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);
|
||||
@@ -128,7 +123,7 @@ namespace Grid{
|
||||
// Momenta
|
||||
for (int mu=0; mu <Nd; mu++){
|
||||
LatticeColourMatrix Pmu = peekLorentz(*P, mu);
|
||||
Hloc -= trace(Pmu*adj(Pmu));
|
||||
Hloc -= trace(Pmu*Pmu);
|
||||
}
|
||||
Complex Hsum = sum(Hloc);
|
||||
|
||||
|
||||
@@ -545,21 +545,21 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
|
||||
|
||||
static void GaussianLieAlgebraMatrix(GridParallelRNG &pRNG,LatticeMatrix &out,double scale=1.0){
|
||||
GridBase *grid = out._grid;
|
||||
LatticeComplex ca (grid);
|
||||
LatticeMatrix lie(grid);
|
||||
LatticeReal ca (grid);
|
||||
LatticeMatrix la (grid);
|
||||
Complex ci(0.0,scale);
|
||||
Matrix ta;
|
||||
|
||||
out=zero;
|
||||
for(int a=0;a<generators();a++){
|
||||
|
||||
gaussian(pRNG,ca);
|
||||
generator(a,ta);
|
||||
|
||||
la=ci*ca*ta;
|
||||
out = lie+la; // e^{i la ta}
|
||||
la=toComplex(ca)*ci*ta;
|
||||
|
||||
out += la;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -69,9 +69,7 @@ namespace Grid {
|
||||
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 = rsqrt(inner);
|
||||
|
||||
for(int c2=0;c2<N;c2++)
|
||||
ret._internal[c1][c2]*= nrm;
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
|
||||
bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
|
||||
bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
|
||||
|
||||
|
||||
Test_cayley_cg_SOURCES=Test_cayley_cg.cc
|
||||
@@ -78,8 +78,8 @@ Test_GaugeAction_SOURCES=Test_GaugeAction.cc
|
||||
Test_GaugeAction_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_hmc_SOURCES=Test_hmc.cc
|
||||
Test_hmc_LDADD=-lGrid
|
||||
Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
|
||||
Test_hmc_WilsonGauge_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_lie_generators_SOURCES=Test_lie_generators.cc
|
||||
|
||||
@@ -26,7 +26,7 @@ int main (int argc, char ** argv)
|
||||
pRNG.SeedRandomDevice();
|
||||
LatticeLorentzColourMatrix U(&Fine);
|
||||
|
||||
SU3::ColdConfiguration(pRNG, U);
|
||||
SU3::HotConfiguration(pRNG, U);
|
||||
|
||||
|
||||
// simplify template?
|
||||
@@ -39,7 +39,7 @@ int main (int argc, char ** argv)
|
||||
FullSet.push_back(Level1);
|
||||
|
||||
// Create integrator
|
||||
IntegratorParameters MDpar(12,50,1.0);
|
||||
IntegratorParameters MDpar(12,30,1.0);
|
||||
std::vector<int> rel ={1};
|
||||
Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);
|
||||
|
||||
Reference in New Issue
Block a user