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

Moved all the HMC tests over to using a single HmcRunner class that manages checkpoint strategies and such like

Browse Source
This commit is contained in:
paboyle 2015-12-22 11:19:25 +00:00
parent 08edbb5cbe
commit 0afcf1cf13
13 changed files with 679 additions and 652 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
lib/Grid.h
View File

@ -47,9 +47,11 @@
#include <parallelIO/BinaryIO.h>
#include <qcd/QCD.h>
#include <parallelIO/NerscIO.h>
#include <qcd/hmc/NerscCheckpointer.h>
#include <Init.h>
#include <qcd/hmc/NerscCheckpointer.h>
#include <qcd/hmc/HmcRunner.h>
#endif

113
lib/qcd/hmc/HmcRunner.h Normal file
View File

@ -0,0 +1,113 @@
#ifndef HMC_RUNNER
#define HMC_RUNNER
namespace Grid{
namespace QCD{
class NerscHmcRunner {
public:
enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
ActionSet<LatticeGaugeField> TheAction;
GridCartesian * UGrid ;
GridCartesian * FGrid ;
GridRedBlackCartesian * UrbGrid ;
GridRedBlackCartesian * FrbGrid ;
virtual void BuildTheAction (int argc, char **argv) = 0;
void Run (int argc, char **argv){
StartType_t StartType = HotStart;
std::string arg;
if( GridCmdOptionExists(argv,argv+argc,"--StartType") ){
arg = GridCmdOptionPayload(argv,argv+argc,"--StartType");
if ( arg == "HotStart" ) { StartType = HotStart; }
else if ( arg == "ColdStart" ) { StartType = ColdStart; }
else if ( arg == "TepidStart" ) { StartType = TepidStart; }
else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
else assert(0);
}
int StartTraj = 0;
if( GridCmdOptionExists(argv,argv+argc,"--StartTrajectory") ){
arg= GridCmdOptionPayload(argv,argv+argc,"--StartTrajectory");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg,ivec);
StartTraj = ivec[0];
}
int NumTraj = 1;
if( GridCmdOptionExists(argv,argv+argc,"--Trajectories") ){
arg= GridCmdOptionPayload(argv,argv+argc,"--Trajectories");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg,ivec);
NumTraj = ivec[0];
}
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorType;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorType MDynamics(UGrid,MDpar, TheAction);
// Checkpoint strategy
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
HMCparameters HMCpar;
HMCpar.StartTrajectory = StartTraj;
HMCpar.Trajectories = NumTraj;
GridSerialRNG sRNG;
GridParallelRNG pRNG(UGrid);
LatticeGaugeField U(UGrid);
std::vector<int> SerSeed({1,2,3,4,5});
std::vector<int> ParSeed({6,7,8,9,10});
if ( StartType == HotStart ) {
// Hot start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::HotConfiguration(pRNG, U);
} else if ( StartType == ColdStart ) {
// Cold start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::ColdConfiguration(pRNG, U);
} else if ( StartType == TepidStart ) {
// Tepid start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::TepidConfiguration(pRNG, U);
} else if ( StartType == CheckpointStart ) {
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
// CheckpointRestart
Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
}
HybridMonteCarlo<LatticeGaugeField,IntegratorType> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
HMC.AddObservable(&PlaqLog);
// Run it
HMC.evolve();
}
};
}}
#endif

116
tests/Test_hmc_EODWFRatio.cc
View File

@ -1,4 +1,5 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
@ -6,20 +7,17 @@ using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class NerscHmcRunner {
class HmcRunner : public NerscHmcRunner {
public:
enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
ActionSet<LatticeGaugeField> TheAction;
GridCartesian * UGrid ;
GridCartesian * FGrid ;
GridRedBlackCartesian * UrbGrid ;
GridRedBlackCartesian * FrbGrid ;
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef DomainWallFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
const int Ls = 8;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@ -37,11 +35,11 @@ public:
Real mass=0.04;
Real pv =1.0;
RealD M5=1.5;
DomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
DomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
@ -51,95 +49,7 @@ public:
Run(argc,argv);
};
void Run (int argc, char **argv){
StartType_t StartType = HotStart;
std::string arg;
if( GridCmdOptionExists(argv,argv+argc,"--StartType") ){
arg = GridCmdOptionPayload(argv,argv+argc,"--StartType");
if ( arg == "HotStart" ) { StartType = HotStart; }
else if ( arg == "ColdStart" ) { StartType = ColdStart; }
else if ( arg == "TepidStart" ) { StartType = TepidStart; }
else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
else assert(0);
}
int StartTraj = 0;
if( GridCmdOptionExists(argv,argv+argc,"--StartTrajectory") ){
arg= GridCmdOptionPayload(argv,argv+argc,"--StartTrajectory");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg,ivec);
StartTraj = ivec[0];
}
int NumTraj = 1;
if( GridCmdOptionExists(argv,argv+argc,"--Trajectories") ){
arg= GridCmdOptionPayload(argv,argv+argc,"--Trajectories");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg,ivec);
NumTraj = ivec[0];
}
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorType;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorType MDynamics(UGrid,MDpar, TheAction);
// Checkpoint strategy
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
HMCparameters HMCpar;
HMCpar.StartTrajectory = StartTraj;
HMCpar.Trajectories = NumTraj;
GridSerialRNG sRNG;
GridParallelRNG pRNG(UGrid);
LatticeGaugeField U(UGrid);
std::vector<int> SerSeed({1,2,3,4,5});
std::vector<int> ParSeed({6,7,8,9,10});
if ( StartType == HotStart ) {
// Hot start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::HotConfiguration(pRNG, U);
} else if ( StartType == ColdStart ) {
// Cold start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::ColdConfiguration(pRNG, U);
} else if ( StartType == TepidStart ) {
// Tepid start
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed);
SU3::TepidConfiguration(pRNG, U);
} else if ( StartType == CheckpointStart ) {
HMCpar.NoMetropolisUntil =0;
HMCpar.MetropolisTest = true;
// CheckpointRestart
Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
}
HybridMonteCarlo<LatticeGaugeField,IntegratorType> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
HMC.AddObservable(&PlaqLog);
// Run it
HMC.evolve();
}
};
}}
@ -151,7 +61,7 @@ int main (int argc, char ** argv)
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
NerscHmcRunner TheHMC;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);

114
tests/Test_hmc_EODWFRatio_Gparity.cc
View File

@ -1,73 +1,69 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef GparityWilsonImplR ImplPolicy;
typedef GparityDomainWallFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
const int Ls = 8;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
Real mass=0.04;
Real pv =1.0;
RealD M5=1.5;
FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&Nf2);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls = 8;
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
GridSerialRNG sRNG;
GridParallelRNG pRNG(UGrid);
sRNG.SeedRandomDevice();
pRNG.SeedRandomDevice();
LatticeLorentzColourMatrix U(UGrid);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=0.04;
Real pv =1.0;
RealD M5=1.5;
typedef typename GparityDomainWallFermionR::FermionField FermionField;
const int nu = 3;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
GparityDomainWallFermionR::ImplParams params;
params.twists = twists;
GparityDomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
GparityDomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,params);
HmcRunner TheHMC;
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<GparityWilsonImplR> Nf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&Nf2);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorType;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorType MDynamics(UGrid,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorType> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

115
tests/Test_hmc_EOWilsonFermionGauge.cc
View File

@ -1,69 +1,72 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourEvenOddPseudoFermionAction<ImplPolicy> Nf2(FermOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<LatticeGaugeField> Level2(4);
Level2.push_back(&Waction);
TheAction.push_back(Level1);
TheAction.push_back(Level2);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
WilsonFermionR FermOp(U,Fine,RBFine,mass);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
TwoFlavourEvenOddPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
ActionLevel<LatticeGaugeField> Level2(4);
Level1.push_back(&WilsonNf2);
Level2.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
FullSet.push_back(Level2);
// Create integrator
//typedef LeapFrog<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
//typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
typedef ForceGradient<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(16);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

107
tests/Test_hmc_EOWilsonRatio.cc
View File

@ -1,67 +1,66 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
RealD mass=-0.77;
RealD pv =0.0;
FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&Nf2);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
Real pv =0.0;
WilsonFermionR DenOp(U,Fine,RBFine,mass);
WilsonFermionR NumOp(U,Fine,RBFine,pv);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf2);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
// typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

112
tests/Test_hmc_WilsonFermionGauge.cc
View File

@ -1,68 +1,70 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<LatticeGaugeField> Level2(4);
Level2.push_back(&Waction);
TheAction.push_back(Level1);
TheAction.push_back(Level2);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
WilsonFermionR FermOp(U,Fine,RBFine,mass);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
TwoFlavourPseudoFermionAction<WilsonImplR> WilsonNf2(FermOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
ActionLevel<LatticeGaugeField> Level2(4);
Level1.push_back(&WilsonNf2);
Level1.push_back(&Waction);
// Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// FullSet.push_back(Level2);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
// typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

91
tests/Test_hmc_WilsonGauge.cc
View File

@ -1,62 +1,57 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
latt_size.resize(4);
latt_size[0] = 8;
latt_size[1] = 8;
latt_size[2] = 8;
latt_size[3] = 8;
double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeGaugeField U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(6.0);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to modify the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

107
tests/Test_hmc_WilsonRatio.cc
View File

@ -1,66 +1,67 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
RealD mass=-0.77;
RealD pv =0.0;
FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
TwoFlavourRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&Nf2);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
Real pv =0.0;
WilsonFermionR FermOp(U,Fine,RBFine,mass);
WilsonFermionR NumOp(U,Fine,RBFine,pv);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
TwoFlavourRatioPseudoFermionAction<WilsonImplR> WilsonNf2(NumOp,FermOp,CG,CG);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf2);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
// typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

110
tests/Test_rhmc_EOWilson1p1.cc
View File

@ -1,67 +1,69 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
// 1+1 flavour
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
WilsonFermionR FermOp(U,Fine,RBFine,mass);
// 1+1 flavour
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

116
tests/Test_rhmc_EOWilsonRatio.cc
View File

@ -1,70 +1,72 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
RealD mass=-0.77;
RealD pv =0.0;
FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
// erange,maxiter,resid,npoly
OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
Real pv =0.0;
WilsonFermionR DenOp(U,Fine,RBFine,mass);
WilsonFermionR NumOp(U,Fine,RBFine,pv);
// erange,maxiter,resid,npoly
OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
OneFlavourEvenOddRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
// typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

110
tests/Test_rhmc_Wilson1p1.cc
View File

@ -1,68 +1,68 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
FermionAction FermOp(U,*FGrid,*FrbGrid,mass);
// 1+1 flavour
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
WilsonFermionR FermOp(U,Fine,RBFine,mass);
// 1+1 flavour
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6);
OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params);
OneFlavourRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}

114
tests/Test_rhmc_WilsonRatio.cc
View File

@ -1,69 +1,71 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
WilsonGaugeActionR Waction(5.6);
RealD mass=-0.77;
RealD pv =0.0;
FermionAction DenOp(U,*FGrid,*FrbGrid,mass);
FermionAction NumOp(U,*FGrid,*FrbGrid,pv);
// erange,maxiter,resid,npoly
OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({6,7,8,80});
GridParallelRNG pRNG(&Fine);
pRNG.SeedFixedIntegers(seeds);
std::vector<int> seedsS({1,2,3,4});
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(seedsS);
LatticeLorentzColourMatrix U(&Fine);
SU3::HotConfiguration(pRNG, U);
// simplify template declaration? Strip the lorentz from the second template
WilsonGaugeActionR Waction(5.6);
Real mass=-0.77;
Real pv =0.0;
WilsonFermionR DenOp(U,Fine,RBFine,mass);
WilsonFermionR NumOp(U,Fine,RBFine,pv);
// erange,maxiter,resid,npoly
OneFlavourRationalParams Params(1.0e-2,64.0,1000,1.0e-6,6);
OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(NumOp,DenOp,Params);
OneFlavourRatioRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(NumOp,DenOp,Params);
//Collect actions
ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction);
ActionSet<LatticeGaugeField> FullSet;
FullSet.push_back(Level1);
// Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorAlgorithm;// change here to change the algorithm
// typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm
IntegratorParameters MDpar(20);
IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
// Create HMC
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
HMCparameters HMCpar;
HybridMonteCarlo<LatticeGaugeField,IntegratorAlgorithm> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint);
// Run it
HMC.evolve();
TheHMC.BuildTheAction(argc,argv);
}