portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-17 07:17:06 +01:00
Code Releases Activity

Merge branch 'feature/gparity_HMC' into feature/ddhmc

Browse Source
This commit is contained in:
Peter Boyle
2021-05-06 20:55:03 +02:00
parent cff884929c 1c70d8c4d9
commit f776a7fe4a
51 changed files with 4291 additions and 426 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
tests/core/Test_fft.cc
View File

@ -299,12 +299,12 @@ int main (int argc, char ** argv)
SpinColourVectorD ferm; gaussian(sRNG,ferm);
pokeSite(ferm,src,point);
const int Ls=32;
const int Ls=64;
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
RealD mass=0.01;
RealD M5 =0.8;
RealD mass=1.0;
RealD M5 =0.99;
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5);
// Momentum space prop
@ -353,6 +353,12 @@ int main (int argc, char ** argv)
std::cout << " Taking difference" <<std::endl;
std::cout << "Ddwf result4 "<<norm2(result4)<<std::endl;
std::cout << "Ddwf ref "<<norm2(ref)<<std::endl;
auto twopoint = localInnerProduct(result4,result4);
std::vector<TComplex> pion_prop;
sliceSum(twopoint,pion_prop,Nd-1);
for(int t=0;t<pion_prop.size();t++){
std::cout << "Pion_prop["<<t<<"]="<<pion_prop[t]<<std::endl;
}
diff = ref - result4;
std::cout << "result - ref "<<norm2(diff)<<std::endl;
@ -383,7 +389,7 @@ int main (int argc, char ** argv)
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
RealD mass=0.01;
RealD mass=1.0;
RealD M5 =0.8;
OverlapWilsonCayleyTanhFermionD Dov(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5,1.0);

36
tests/core/Test_gparity.cc
View File

@ -55,13 +55,17 @@ static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIM
int main (int argc, char ** argv)
{
int nu = 0;
int tbc_aprd = 0; //use antiperiodic BCs in the time direction?
Grid_init(&argc,&argv);
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "--Gparity-dir"){
std::stringstream ss; ss << argv[i+1]; ss >> nu;
std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl;
}else if(std::string(argv[i]) == "--Tbc-APRD"){
tbc_aprd = 1;
std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl;
}
}
@ -155,13 +159,18 @@ int main (int argc, char ** argv)
//Coordinate grid for reference
LatticeInteger xcoor_1f5(FGrid_1f);
LatticeCoordinate(xcoor_1f5,1+nu);
LatticeCoordinate(xcoor_1f5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0
Replicate(src,src_1f);
src_1f = where( xcoor_1f5 >= Integer(L), 2.0*src_1f,src_1f );
RealD mass=0.0;
RealD M5=1.8;
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS);
//Standard Dirac op
AcceleratorVector<Complex,4> bc_std(Nd, 1.0);
if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC
StandardDiracOp::ImplParams std_params(bc_std);
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS, std_params);
StandardFermionField src_o_1f(FrbGrid_1f);
StandardFermionField result_o_1f(FrbGrid_1f);
@ -172,9 +181,11 @@ int main (int argc, char ** argv)
ConjugateGradient<StandardFermionField> CG(1.0e-8,10000);
CG(HermOpEO,src_o_1f,result_o_1f);
// const int nu = 3;
//Gparity Dirac op
std::vector<int> twists(Nd,0);
twists[nu] = 1;
if(tbc_aprd) twists[Nd-1] = 1;
GparityDiracOp::ImplParams params;
params.twists = twists;
GparityDiracOp GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5 DOP_PARAMS,params);
@ -271,8 +282,11 @@ int main (int argc, char ** argv)
std::cout << "2f cb "<<result_o_2f.Checkerboard()<<std::endl;
std::cout << "1f cb "<<result_o_1f.Checkerboard()<<std::endl;
std::cout << " result norms " <<norm2(result_o_2f)<<" " <<norm2(result_o_1f)<<std::endl;
//Compare norms
std::cout << " result norms 2f: " <<norm2(result_o_2f)<<" 1f: " <<norm2(result_o_1f)<<std::endl;
//Take the 2f solution and convert into the corresponding 1f solution (odd cb only)
StandardFermionField res0o (FrbGrid_2f);
StandardFermionField res1o (FrbGrid_2f);
StandardFermionField res0 (FGrid_2f);
@ -281,14 +295,15 @@ int main (int argc, char ** argv)
res0=Zero();
res1=Zero();
res0o = PeekIndex<0>(result_o_2f,0);
res1o = PeekIndex<0>(result_o_2f,1);
res0o = PeekIndex<0>(result_o_2f,0); //flavor 0, odd cb
res1o = PeekIndex<0>(result_o_2f,1); //flavor 1, odd cb
std::cout << "res cb "<<res0o.Checkerboard()<<std::endl;
std::cout << "res cb "<<res1o.Checkerboard()<<std::endl;
setCheckerboard(res0,res0o);
setCheckerboard(res1,res1o);
//poke odd onto non-cb field
setCheckerboard(res0,res0o);
setCheckerboard(res1,res1o);
StandardFermionField replica (FGrid_1f);
StandardFermionField replica0(FGrid_1f);
@ -296,12 +311,13 @@ int main (int argc, char ** argv)
Replicate(res0,replica0);
Replicate(res1,replica1);
//2nd half of doubled lattice has f=1
replica = where( xcoor_1f5 >= Integer(L), replica1,replica0 );
replica0 = Zero();
setCheckerboard(replica0,result_o_1f);
std::cout << "Norm2 solutions is " <<norm2(replica)<<" "<< norm2(replica0)<<std::endl;
std::cout << "Norm2 solutions 1f reconstructed from 2f: " <<norm2(replica)<<" Actual 1f: "<< norm2(replica0)<<std::endl;
replica = replica - replica0;

177
tests/core/Test_gparity_flavour.cc Normal file
View File

@ -0,0 +1,177 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_gparity_flavour.cc
Copyright (C) 2015-2017
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
static constexpr double tolerance = 1.0e-6;
static std::array<GparityFlavourMatrix, GparityFlavour::nSigma> testAlgebra;
void print(const GparityFlavourMatrix &g)
{
for(int i = 0; i < Ngp; i++)
{
std::cout << GridLogMessage << "(";
for(int j=0;j<Ngp;j++){
if ( abs( g(i,j)()() ) == 0 ) {
std::cout<< " 0";
} else if ( abs(g(i,j)()() - Complex(0,1)) == 0){
std::cout<< " i";
} else if ( abs(g(i,j)()() + Complex(0,1)) == 0){
std::cout<< "-i";
} else if ( abs(g(i,j)()() - Complex(1,0)) == 0){
std::cout<< " 1";
} else if ( abs(g(i,j)()() + Complex(1,0)) == 0){
std::cout<< "-1";
}
std::cout<<((j == Ngp-1) ? ")" : "," );
}
std::cout << std::endl;
}
std::cout << GridLogMessage << std::endl;
}
void createTestAlgebra(void)
{
std::array<GparityFlavourMatrix, 3> testg;
const Complex I(0., 1.), mI(0., -1.);
// 0 1
// 1 0
testg[0] = Zero();
testg[0](0, 1)()() = 1.;
testg[0](1, 0)()() = 1.;
std::cout << GridLogMessage << "test SigmaX= " << std::endl;
print(testg[0]);
// 0 -i
// i 0
testg[1] = Zero();
testg[1](0, 1)()() = mI;
testg[1](1, 0)()() = I;
std::cout << GridLogMessage << "test SigmaY= " << std::endl;
print(testg[1]);
// 1 0
// 0 -1
testg[2] = Zero();
testg[2](0, 0)()() = 1.0;
testg[2](1, 1)()() = -1.0;
std::cout << GridLogMessage << "test SigmaZ= " << std::endl;
print(testg[2]);
#define DEFINE_TEST_G(g, exp)\
testAlgebra[GparityFlavour::Algebra::g] = exp; \
testAlgebra[GparityFlavour::Algebra::Minus##g] = -exp;
DEFINE_TEST_G(SigmaX , testg[0]);
DEFINE_TEST_G(SigmaY , testg[1]);
DEFINE_TEST_G(SigmaZ , testg[2]);
DEFINE_TEST_G(Identity , 1.);
GparityFlavourMatrix pplus;
pplus = 1.0;
pplus = pplus + testg[1];
pplus = pplus * 0.5;
DEFINE_TEST_G(ProjPlus , pplus);
GparityFlavourMatrix pminus;
pminus = 1.0;
pminus = pminus - testg[1];
pminus = pminus * 0.5;
DEFINE_TEST_G(ProjMinus , pminus);
#undef DEFINE_TEST_G
}
template <typename Expr>
void test(const Expr &a, const Expr &b)
{
if (norm2(a - b) < tolerance)
{
std::cout << "[OK] ";
}
else
{
std::cout << "[fail]" << std::endl;
std::cout << GridLogError << "a= " << a << std::endl;
std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
std::cout << GridLogError << "b= " << b << std::endl;
exit(EXIT_FAILURE);
}
}
void checkSigma(const GparityFlavour::Algebra a, GridSerialRNG &rng)
{
GparityFlavourVector v;
GparityFlavourMatrix m, &testg = testAlgebra[a];
GparityFlavour g(a);
random(rng, v);
random(rng, m);
std::cout << GridLogMessage << "Checking " << GparityFlavour::name[a] << ": ";
std::cout << "vecmul ";
test(g*v, testg*v);
std::cout << "matlmul ";
test(g*m, testg*m);
std::cout << "matrmul ";
test(m*g, m*testg);
std::cout << std::endl;
}
int main(int argc, char *argv[])
{
Grid_init(&argc,&argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
std::cout << GridLogMessage << "======== Test algebra" << std::endl;
createTestAlgebra();
std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
for (int i = 0; i < GparityFlavour::nSigma; ++i)
{
checkSigma(i, sRNG);
}
std::cout << GridLogMessage << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

114
tests/core/Test_precision_change.cc Normal file
View File

@ -0,0 +1,114 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/core/Test_precision_change.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
int main (int argc, char ** argv){
Grid_init(&argc, &argv);
int Ls = 16;
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << " and Ls=" << Ls << std::endl;
GridCartesian* UGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
GridCartesian* FGrid_d = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_d);
GridRedBlackCartesian* FrbGrid_d = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_d);
GridCartesian* UGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
GridCartesian* FGrid_f = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_f);
GridRedBlackCartesian* FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_f);
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid_d);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid_d);
RNG4.SeedFixedIntegers(seeds4);
//Gauge fields
LatticeGaugeFieldD Umu_d(UGrid_d);
LatticeGaugeFieldF Umu_f(UGrid_f);
LatticeGaugeFieldD Umu_d_r(UGrid_d);
LatticeGaugeFieldD Utmp_d(UGrid_d);
for(int i=0;i<5;i++){
random(RNG4, Umu_d);
precisionChange(Umu_f, Umu_d);
std::cout << GridLogMessage << "Norm of double-prec and single-prec gauge fields (should be ~equal): " << norm2(Umu_d) << " " << norm2(Umu_f) << std::endl;
precisionChange(Umu_d_r, Umu_f);
RealD normdiff = axpy_norm(Utmp_d, -1.0, Umu_d_r, Umu_d);
std::cout << GridLogMessage << "Norm of difference of back-converted double-prec gauge fields (should be ~0) = " << normdiff << std::endl;
}
//Fermion fields
LatticeFermionD psi_d(FGrid_d);
LatticeFermionF psi_f(FGrid_f);
LatticeFermionD psi_d_r(FGrid_d);
LatticeFermionD psi_tmp_d(FGrid_d);
for(int i=0;i<5;i++){
random(RNG5, psi_d);
precisionChange(psi_f, psi_d);
std::cout << GridLogMessage << "Norm of double-prec and single-prec fermion fields (should be ~equal): " << norm2(psi_d) << " " << norm2(psi_f) << std::endl;
precisionChange(psi_d_r, psi_f);
RealD normdiff = axpy_norm(psi_tmp_d, -1.0, psi_d_r, psi_d);
std::cout << GridLogMessage << "Norm of difference of back-converted double-prec fermion fields (should be ~0)= " << normdiff << std::endl;
}
//Checkerboarded fermion fields
LatticeFermionD psi_cb_d(FrbGrid_d);
LatticeFermionF psi_cb_f(FrbGrid_f);
LatticeFermionD psi_cb_d_r(FrbGrid_d);
LatticeFermionD psi_cb_tmp_d(FrbGrid_d);
for(int i=0;i<5;i++){
random(RNG5, psi_d);
pickCheckerboard(Odd, psi_cb_d, psi_d);
precisionChange(psi_cb_f, psi_cb_d);
std::cout << GridLogMessage << "Norm of odd-cb double-prec and single-prec fermion fields (should be ~equal): " << norm2(psi_cb_d) << " " << norm2(psi_cb_f) << std::endl;
precisionChange(psi_cb_d_r, psi_cb_f);
RealD normdiff = axpy_norm(psi_cb_tmp_d, -1.0, psi_cb_d_r, psi_cb_d);
std::cout << GridLogMessage << "Norm of difference of back-converted odd-cb double-prec fermion fields (should be ~0)= " << normdiff << std::endl;
pickCheckerboard(Even, psi_cb_d, psi_d);
precisionChange(psi_cb_f, psi_cb_d);
std::cout << GridLogMessage << "Norm of even-cb double-prec and single-prec fermion fields (should be ~equal): " << norm2(psi_cb_d) << " " << norm2(psi_cb_f) << std::endl;
precisionChange(psi_cb_d_r, psi_cb_f);
normdiff = axpy_norm(psi_cb_tmp_d, -1.0, psi_cb_d_r, psi_cb_d);
std::cout << GridLogMessage << "Norm of difference of back-converted even-cb double-prec fermion fields (should be ~0)= " << normdiff << std::endl;
}
Grid_finalize();
}

18
tests/forces/Test_dwf_gpforce.cc
View File

@ -71,26 +71,14 @@ int main (int argc, char ** argv)
////////////////////////////////////
RealD mass=0.2; //kills the diagonal term
RealD M5=1.8;
// const int nu = 3;
// std::vector<int> twists(Nd,0); // twists[nu] = 1;
// GparityDomainWallFermionR::ImplParams params; params.twists = twists;
// GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
// DomainWallFermionR Dw (U, Grid,RBGrid,mass,M5);
const int nu = 3;
const int nu = 0; //gparity direction
std::vector<int> twists(Nd,0);
twists[nu] = 1;
twists[Nd-1] = 1; //antiperiodic in time
GparityDomainWallFermionR::ImplParams params;
params.twists = twists;
/*
params.boundary_phases[0] = 1.0;
params.boundary_phases[1] = 1.0;
params.boundary_phases[2] = 1.0;
params.boundary_phases[3] =- 1.0;
*/
GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
Dw.M (phi,Mphi);

28
tests/forces/Test_gpdwf_force.cc
View File

@ -71,8 +71,10 @@ int main (int argc, char ** argv)
RealD mass=0.01;
RealD M5=1.8;
const int nu = 3;
std::vector<int> twists(Nd,0); twists[nu] = 1;
const int nu = 1;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
twists[3] = 1;
GparityDomainWallFermionR::ImplParams params; params.twists = twists;
GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
Ddwf.M (phi,Mphi);
@ -91,16 +93,28 @@ int main (int argc, char ** argv)
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
RealD dt = 0.01;
LatticeColourMatrix zz(UGrid); zz=Zero();
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
const int Lnu=latt_size[nu];
Lattice<iScalar<vInteger> > coor(UGrid);
LatticeCoordinate(coor,nu);
for(int mu=0;mu<Nd;mu++){
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
// Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
if(0){
if(mu==nu){
mommu=where(coor==Lnu-1,mommu,zz);
} else {
mommu=Zero();
}
}
PokeIndex<LorentzIndex>(mom,mommu,mu);
@ -125,6 +139,12 @@ int main (int argc, char ** argv)
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);
LatticeComplex lip(FGrid); lip=localInnerProduct(Mphi,Mphi);
LatticeComplex lipp(FGrid); lipp=localInnerProduct(MphiPrime,MphiPrime);
LatticeComplex dip(FGrid); dip = lipp - lip;
std::cout << " dip "<<dip<<std::endl;
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////

35
tests/forces/Test_gpwilson_force.cc
View File

@ -64,8 +64,12 @@ int main (int argc, char ** argv)
////////////////////////////////////
RealD mass=0.01;
const int nu = 3;
std::vector<int> twists(Nd,0); twists[nu] = 1;
const int nu = 1;
const int Lnu=latt_size[nu];
std::vector<int> twists(Nd,0);
twists[nu] = 1;
twists[3]=1;
GparityWilsonFermionR::ImplParams params; params.twists = twists;
GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params);
Wil.M (phi,Mphi);
@ -87,17 +91,28 @@ int main (int argc, char ** argv)
RealD dt = 0.01;
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix zz(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
Lattice<iScalar<vInteger> > coor(UGrid);
LatticeCoordinate(coor,nu);
zz=Zero();
for(int mu=0;mu<Nd;mu++){
// Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
if(0){
if(mu==nu){
mommu=where(coor==Lnu-1,mommu,zz);
} else {
mommu=Zero();
}
}
PokeIndex<LorentzIndex>(mom,mommu,mu);
// fourth order exponential approx
autoView( mom_v, mom, CpuRead);
autoView( U_v , U, CpuRead);
@ -130,6 +145,10 @@ int main (int argc, char ** argv)
mommu=Ta(mommu)*2.0;
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
LatticeComplex lip(UGrid); lip=localInnerProduct(Mphi,Mphi);
LatticeComplex lipp(UGrid); lipp=localInnerProduct(MphiPrime,MphiPrime);
LatticeComplex dip(UGrid); dip = lipp - lip;
std::cout << " dip "<<dip<<std::endl;
LatticeComplex dS(UGrid); dS = Zero();
for(int mu=0;mu<Nd;mu++){
@ -139,12 +158,14 @@ int main (int argc, char ** argv)
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
}
std::cout << "mommu"<<mommu<<std::endl;
std::cout << "dS" << dS<<std::endl;
ComplexD dSpred = sum(dS);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
assert( fabs(real(Sprime-S-dSpred)) < 2.0 ) ;

257
tests/hmc/Test_action_dwf_gparity2fvs1f.cc Normal file
View File

@ -0,0 +1,257 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hmc/Test_action_dwf_gparity2fvs1f.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
template<typename FermionField2f, typename FermionField1f>
void copy2fTo1fFermionField(FermionField1f &out, const FermionField2f &in, int gpdir){
auto f0_halfgrid = PeekIndex<GparityFlavourIndex>(in,0); //on 2f Grid
FermionField1f f0_fullgrid_dbl(out.Grid());
Replicate(f0_halfgrid, f0_fullgrid_dbl); //double it up to live on the 1f Grid
auto f1_halfgrid = PeekIndex<GparityFlavourIndex>(in,1);
FermionField1f f1_fullgrid_dbl(out.Grid());
Replicate(f1_halfgrid, f1_fullgrid_dbl);
const Coordinate &dim_2f = in.Grid()->GlobalDimensions();
const Coordinate &dim_1f = out.Grid()->GlobalDimensions();
//We have to be careful for 5d fields; the s-direction is placed before the x,y,z,t and so we need to shift gpdir by 1
std::cout << "gpdir " << gpdir << std::endl;
gpdir+=1;
std::cout << "gpdir for 5D fields " << gpdir << std::endl;
std::cout << "dim_2f " << dim_2f << std::endl;
std::cout << "dim_1f " << dim_1f << std::endl;
assert(dim_1f[gpdir] == 2*dim_2f[gpdir]);
LatticeInteger xcoor_1f(out.Grid()); //5d lattice integer
LatticeCoordinate(xcoor_1f,gpdir);
int L = dim_2f[gpdir];
out = where(xcoor_1f < L, f0_fullgrid_dbl, f1_fullgrid_dbl);
}
//Both have the same field type
void copy2fTo1fGaugeField(LatticeGaugeField &out, const LatticeGaugeField &in, int gpdir){
LatticeGaugeField U_dbl(out.Grid());
Replicate(in, U_dbl);
LatticeGaugeField Uconj_dbl = conjugate( U_dbl );
const Coordinate &dim_2f = in.Grid()->GlobalDimensions();
LatticeInteger xcoor_1f(out.Grid());
LatticeCoordinate(xcoor_1f,gpdir);
int L = dim_2f[gpdir];
out = where(xcoor_1f < L, U_dbl, Uconj_dbl);
}
std::ostream & operator<<(std::ostream &os, const Coordinate &x){
os << "(";
for(int i=0;i<x.size();i++) os << x[i] << (i<x.size()-1 ? " " : "");
os << ")";
return os;
}
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
int Ls = 16;
Coordinate latt_2f = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd, vComplexD::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
int mu = 0; //Gparity direction
Coordinate latt_1f = latt_2f;
latt_1f[mu] *= 2;
GridCartesian * UGrid_1f = SpaceTimeGrid::makeFourDimGrid(latt_1f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_1f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_1f);
GridCartesian * FGrid_1f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_1f);
GridRedBlackCartesian * FrbGrid_1f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_1f);
GridCartesian * UGrid_2f = SpaceTimeGrid::makeFourDimGrid(latt_2f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_2f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_2f);
GridCartesian * FGrid_2f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_2f);
GridRedBlackCartesian * FrbGrid_2f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_2f);
std::cout << "SIMD layout " << simd_layout << std::endl;
std::cout << "MPI layout " << mpi_layout << std::endl;
std::cout << "2f dimensions " << latt_2f << std::endl;
std::cout << "1f dimensions " << latt_1f << std::endl;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4);
std::cout << "Generating hot 2f gauge configuration" << std::endl;
LatticeGaugeField Umu_2f(UGrid_2f);
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
std::cout << "Copying 2f->1f gauge field" << std::endl;
LatticeGaugeField Umu_1f(UGrid_1f);
copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu);
typedef GparityWilsonImplR FermionImplPolicy2f;
typedef GparityDomainWallFermionR FermionAction2f;
typedef typename FermionAction2f::FermionField FermionField2f;
typedef WilsonImplR FermionImplPolicy1f;
typedef DomainWallFermionR FermionAction1f;
typedef typename FermionAction1f::FermionField FermionField1f;
std::cout << "Generating eta 2f" << std::endl;
FermionField2f eta_2f(FGrid_2f);
gaussian(RNG5_2f, eta_2f);
RealD scale = std::sqrt(0.5);
eta_2f=eta_2f*scale;
std::cout << "Copying 2f->1f eta" << std::endl;
FermionField1f eta_1f(FGrid_1f);
copy2fTo1fFermionField(eta_1f, eta_2f, mu);
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.032;
Real pv_mass = 1.0;
RealD M5 = 1.8;
//Setup the Dirac operators
std::cout << "Initializing Dirac operators" << std::endl;
FermionAction2f::ImplParams Params_2f;
Params_2f.twists[mu] = 1;
Params_2f.twists[Nd-1] = 1; //APBC in time direction
//note 'Num' and 'Den' here refer to the determinant ratio, not the operator ratio in the pseudofermion action where the two are inverted
//to my mind the Pauli Villars and 'denominator' are synonymous but the Grid convention has this as the 'Numerator' operator in the RHMC implementation
FermionAction2f NumOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, light_mass,M5,Params_2f);
FermionAction2f DenOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, pv_mass, M5,Params_2f);
FermionAction1f::ImplParams Params_1f;
Params_1f.boundary_phases[mu] = -1; //antiperiodic in doubled lattice in GP direction
Params_1f.boundary_phases[Nd-1] = -1;
FermionAction1f NumOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, light_mass,M5,Params_1f);
FermionAction1f DenOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, pv_mass, M5,Params_1f);
//Test the replication routines by running a CG on eta
double StoppingCondition = 1e-10;
double MaxCGIterations = 30000;
ConjugateGradient<FermionField2f> CG_2f(StoppingCondition,MaxCGIterations);
ConjugateGradient<FermionField1f> CG_1f(StoppingCondition,MaxCGIterations);
NumOp_1f.ImportGauge(Umu_1f);
NumOp_2f.ImportGauge(Umu_2f);
FermionField1f test_1f(FGrid_1f);
FermionField2f test_2f(FGrid_2f);
MdagMLinearOperator<FermionAction1f, FermionField1f> Linop_1f(NumOp_1f);
MdagMLinearOperator<FermionAction2f, FermionField2f> Linop_2f(NumOp_2f);
CG_1f(Linop_1f, eta_1f, test_1f);
CG_2f(Linop_2f, eta_2f, test_2f);
RealD test_1f_norm = norm2(test_1f);
RealD test_2f_norm = norm2(test_2f);
std::cout << "Verification of replication routines: " << test_1f_norm << " " << test_2f_norm << " " << test_1f_norm - test_2f_norm << std::endl;
#if 1
typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy2f> Action2f;
typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy1f> Action1f;
RationalActionParams rational_params;
rational_params.inv_pow = 2;
rational_params.lo = 1e-5;
rational_params.hi = 32;
rational_params.md_degree = 16;
rational_params.action_degree = 16;
Action2f action_2f(DenOp_2f, NumOp_2f, rational_params);
Action1f action_1f(DenOp_1f, NumOp_1f, rational_params);
#else
typedef TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy2f> Action2f;
typedef TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy1f> Action1f;
Action2f action_2f(DenOp_2f, NumOp_2f, CG_2f, CG_2f);
Action1f action_1f(DenOp_1f, NumOp_1f, CG_1f, CG_1f);
#endif
std::cout << "Action refresh" << std::endl;
action_2f.refresh(Umu_2f, eta_2f);
action_1f.refresh(Umu_1f, eta_1f);
std::cout << "Action compute post heatbath" << std::endl;
RealD S_2f = action_2f.S(Umu_2f);
RealD S_1f = action_1f.S(Umu_1f);
std::cout << "Action comparison post heatbath" << std::endl;
std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl;
//Change the gauge field between refresh and action eval else the matrix and inverse matrices all cancel and we just get |eta|^2
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu);
//Now compute the action with the new gauge field
std::cout << "Action compute post gauge field update" << std::endl;
S_2f = action_2f.S(Umu_2f);
S_1f = action_1f.S(Umu_1f);
std::cout << "Action comparison post gauge field update" << std::endl;
std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl;
Grid_finalize();
} // main

4
tests/hmc/Test_hmc_GparityIwasakiGauge.cc
View File

@ -58,7 +58,7 @@ int main(int argc, char **argv) {
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 5;
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@ -79,7 +79,7 @@ int main(int argc, char **argv) {
// that have a complex construction
// standard
RealD beta = 2.6 ;
const int nu = 3;
const int nu = 1;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);

139
tests/hmc/Test_rhmc_EOWilsonRatioPowQuarter.cc Normal file
View File

@ -0,0 +1,139 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_rhmc_EOWilsonRatio.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
//This test is for the Wilson action with the determinant det( M^dag M)^1/4
//testing the generic RHMC
int main(int argc, char **argv) {
using namespace Grid;
;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplR FermionImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 5.6 ;
WilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
Real mass = -0.77;
Real pv = 0.0;
// Can we define an overloaded operator that does not need U and initialises
// it with zeroes?
FermionAction DenOp(U, *GridPtr, *GridRBPtr, mass);
FermionAction NumOp(U, *GridPtr, *GridRBPtr, pv);
// 1/2+1/2 flavour
// RationalActionParams(int _inv_pow = 2,
// RealD _lo = 0.0,
// RealD _hi = 1.0,
// int _maxit = 1000,
// RealD tol = 1.0e-8,
// int _degree = 10,
// int _precision = 64,
// int _BoundsCheckFreq=20)
int inv_pow = 4;
RationalActionParams Params(inv_pow,1.0e-2,64.0,1000,1.0e-6,14,64,1);
GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> RHMC(NumOp,DenOp,Params);
// Collect actions
ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&RHMC);
ActionLevel<HMCWrapper::Field> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run();
Grid_finalize();
} // main

119
tests/hmc/Test_rhmc_EOWilsonRatio_doubleVsMixedPrec.cc Normal file
View File

@ -0,0 +1,119 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_rhmc_EOWilsonRatio_doubleVsMixedPrec.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
//This test ensures the mixed precision RHMC gives the same result as the regular double precision
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplD FermionImplPolicyD;
typedef WilsonFermionD FermionActionD;
typedef typename FermionActionD::FermionField FermionFieldD;
typedef WilsonImplF FermionImplPolicyF;
typedef WilsonFermionF FermionActionF;
typedef typename FermionActionF::FermionField FermionFieldF;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
TheHMC.Resources.AddFourDimGrid("gauge");
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
auto GridPtrD = TheHMC.Resources.GetCartesian();
auto GridRBPtrD = TheHMC.Resources.GetRBCartesian();
GridCartesian* GridPtrF = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF);
// temporarily need a gauge field
LatticeGaugeFieldD Ud(GridPtrD);
LatticeGaugeFieldF Uf(GridPtrF);
Real mass = -0.77;
Real pv = 0.0;
FermionActionD DenOpD(Ud, *GridPtrD, *GridRBPtrD, mass);
FermionActionD NumOpD(Ud, *GridPtrD, *GridRBPtrD, pv);
FermionActionF DenOpF(Uf, *GridPtrF, *GridRBPtrF, mass);
FermionActionF NumOpF(Uf, *GridPtrF, *GridRBPtrF, pv);
TheHMC.Resources.AddRNGs();
PeriodicGimplR::HotConfiguration(TheHMC.Resources.GetParallelRNG(), Ud);
std::string seed_string = "the_seed";
//Setup the pseudofermion actions
RationalActionParams GenParams;
GenParams.inv_pow = 2;
GenParams.lo = 1e-2;
GenParams.hi = 64.0;
GenParams.MaxIter = 1000;
GenParams.action_tolerance = GenParams.md_tolerance = 1e-6;
GenParams.action_degree = GenParams.md_degree = 6;
GenParams.precision = 64;
GenParams.BoundsCheckFreq = 20;
GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicyD> GenD(NumOpD,DenOpD,GenParams);
GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicyD, FermionImplPolicyF> GenFD(NumOpD, DenOpD,
NumOpF, DenOpF,
GenParams, 50);
TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
GenD.refresh(Ud, TheHMC.Resources.GetSerialRNG(), TheHMC.Resources.GetParallelRNG());
RealD Sd = GenD.S(Ud);
LatticeGaugeField derivD(Ud);
GenD.deriv(Ud,derivD);
TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
GenFD.refresh(Ud, TheHMC.Resources.GetSerialRNG(), TheHMC.Resources.GetParallelRNG());
RealD Sfd = GenFD.S(Ud);
LatticeGaugeField derivFD(Ud);
GenFD.deriv(Ud,derivFD);
//Compare
std::cout << "Action : " << Sd << " " << Sfd << " reldiff " << (Sd - Sfd)/Sd << std::endl;
LatticeGaugeField diff(Ud);
axpy(diff, -1.0, derivD, derivFD);
std::cout << "Norm of difference in deriv " << sqrt(norm2(diff)) << std::endl;
Grid_finalize();
return 0;
}

122
tests/hmc/Test_rhmc_EOWilsonRatio_genericVsOneFlavor.cc Normal file
View File

@ -0,0 +1,122 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_rhmc_EOWilsonRatio_genericVsOneFlavor.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Christopher Kelly <ckelly@bnl.gov>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
//This test ensures that the OneFlavourEvenOddRatioRationalPseudoFermionAction and GeneralEvenOddRatioRationalPseudoFermionAction action (with parameters set appropriately0
//give the same results
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplR FermionImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
TheHMC.Resources.AddFourDimGrid("gauge");
// // Checkpointer definition
// CheckpointerParameters CPparams;
// CPparams.config_prefix = "ckpoint_lat";
// CPparams.rng_prefix = "ckpoint_rng";
// CPparams.saveInterval = 5;
// CPparams.format = "IEEE64BIG";
// TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
Real mass = -0.77;
Real pv = 0.0;
FermionAction DenOp(U, *GridPtr, *GridRBPtr, mass);
FermionAction NumOp(U, *GridPtr, *GridRBPtr, pv);
TheHMC.Resources.AddRNGs();
PeriodicGimplR::HotConfiguration(TheHMC.Resources.GetParallelRNG(), U);
std::string seed_string = "the_seed";
//1f action
OneFlavourRationalParams OneFParams(1.0e-2,64.0,1000,1.0e-6,6);
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> OneF(NumOp,DenOp,OneFParams);
TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
OneF.refresh(U, TheHMC.Resources.GetParallelRNG());
RealD OneFS = OneF.S(U);
LatticeGaugeField OneFderiv(U);
OneF.deriv(U,OneFderiv);
//general action
RationalActionParams GenParams;
GenParams.inv_pow = 2;
GenParams.lo = OneFParams.lo;
GenParams.hi = OneFParams.hi;
GenParams.MaxIter = OneFParams.MaxIter;
GenParams.action_tolerance = GenParams.md_tolerance = OneFParams.tolerance;
GenParams.action_degree = GenParams.md_degree = OneFParams.degree;
GenParams.precision = OneFParams.precision;
GenParams.BoundsCheckFreq = OneFParams.BoundsCheckFreq;
GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> Gen(NumOp,DenOp,GenParams);
TheHMC.Resources.GetParallelRNG().SeedUniqueString(seed_string);
Gen.refresh(U, TheHMC.Resources.GetParallelRNG());
RealD GenS = Gen.S(U);
LatticeGaugeField Genderiv(U);
Gen.deriv(U,Genderiv);
//Compare
std::cout << "Action : " << OneFS << " " << GenS << " reldiff " << (OneFS - GenS)/OneFS << std::endl;
LatticeGaugeField diff(U);
axpy(diff, -1.0, Genderiv, OneFderiv);
std::cout << "Norm of difference in deriv " << sqrt(norm2(diff)) << std::endl;
Grid_finalize();
return 0;
}

184
tests/solver/Test_dwf_multishift_mixedprec.cc Normal file
View File

@ -0,0 +1,184 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_multishift_mixedprec.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
template<typename SpeciesD, typename SpeciesF, typename GaugeStatisticsType>
void run_test(int argc, char ** argv, const typename SpeciesD::ImplParams &params){
const int Ls = 16;
GridCartesian* UGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* UrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_d);
GridCartesian* FGrid_d = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_d);
GridRedBlackCartesian* FrbGrid_d = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_d);
GridCartesian* UGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
GridCartesian* FGrid_f = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_f);
GridRedBlackCartesian* FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_f);
typedef typename SpeciesD::FermionField FermionFieldD;
typedef typename SpeciesF::FermionField FermionFieldF;
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid_d);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid_d);
RNG4.SeedFixedIntegers(seeds4);
FermionFieldD src_d(FGrid_d);
random(RNG5, src_d);
LatticeGaugeFieldD Umu_d(UGrid_d);
//CPS-created G-parity ensembles have a factor of 2 error in the plaquette that causes the read to fail unless we workaround it
bool gparity_plaquette_fix = false;
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "--gparity_plaquette_fix"){
gparity_plaquette_fix=true;
break;
}
}
bool cfg_loaded=false;
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "--load_config"){
assert(i != argc-1);
std::string file = argv[i+1];
NerscIO io;
FieldMetaData metadata;
if(gparity_plaquette_fix) NerscIO::exitOnReadPlaquetteMismatch() = false;
io.readConfiguration<GaugeStatisticsType>(Umu_d, metadata, file);
if(gparity_plaquette_fix){
metadata.plaquette *= 2.; //correct header value
//Get the true plaquette
FieldMetaData tmp;
GaugeStatisticsType gs; gs(Umu_d, tmp);
std::cout << "After correction: plaqs " << tmp.plaquette << " " << metadata.plaquette << std::endl;
assert(fabs(tmp.plaquette -metadata.plaquette ) < 1.0e-5 );
}
cfg_loaded=true;
break;
}
}
if(!cfg_loaded)
SU<Nc>::HotConfiguration(RNG4, Umu_d);
LatticeGaugeFieldF Umu_f(UGrid_f);
precisionChange(Umu_f, Umu_d);
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << " Ls: " << Ls << std::endl;
RealD mass = 0.01;
RealD M5 = 1.8;
SpeciesD Ddwf_d(Umu_d, *FGrid_d, *FrbGrid_d, *UGrid_d, *UrbGrid_d, mass, M5, params);
SpeciesF Ddwf_f(Umu_f, *FGrid_f, *FrbGrid_f, *UGrid_f, *UrbGrid_f, mass, M5, params);
FermionFieldD src_o_d(FrbGrid_d);
pickCheckerboard(Odd, src_o_d, src_d);
SchurDiagMooeeOperator<SpeciesD, FermionFieldD> HermOpEO_d(Ddwf_d);
SchurDiagMooeeOperator<SpeciesF, FermionFieldF> HermOpEO_f(Ddwf_f);
AlgRemez remez(1e-4, 64, 50);
int order = 15;
remez.generateApprox(order, 1, 2); //sqrt
MultiShiftFunction shifts(remez, 1e-10, false);
int relup_freq = 50;
double t1=usecond();
ConjugateGradientMultiShiftMixedPrec<FermionFieldD,FermionFieldF> mcg(10000, shifts, FrbGrid_f, HermOpEO_f, relup_freq);
std::vector<FermionFieldD> results_o_d(order, FrbGrid_d);
mcg(HermOpEO_d, src_o_d, results_o_d);
double t2=usecond();
//Crosscheck double and mixed prec results
ConjugateGradientMultiShift<FermionFieldD> dmcg(10000, shifts);
std::vector<FermionFieldD> results_o_d_2(order, FrbGrid_d);
dmcg(HermOpEO_d, src_o_d, results_o_d_2);
double t3=usecond();
std::cout << GridLogMessage << "Comparison of mixed prec results to double prec results |mixed - double|^2 :" << std::endl;
FermionFieldD tmp(FrbGrid_d);
for(int i=0;i<order;i++){
RealD ndiff = axpy_norm(tmp, -1., results_o_d[i], results_o_d_2[i]);
std::cout << i << " " << ndiff << std::endl;
}
std::cout<<GridLogMessage << "Mixed precision algorithm: Total usec = "<< (t2-t1)<<std::endl;
std::cout<<GridLogMessage << "Double precision algorithm: Total usec = "<< (t3-t2)<<std::endl;
}
int main (int argc, char ** argv)
{
Grid_init(&argc, &argv);
bool gparity = false;
int gpdir;
for(int i=1;i<argc;i++){
std::string arg(argv[i]);
if(arg == "--Gparity"){
assert(i!=argc-1);
gpdir = std::stoi(argv[i+1]);
assert(gpdir >= 0 && gpdir <= 2); //spatial!
gparity = true;
}
}
if(gparity){
std::cout << "Running test with G-parity BCs in " << gpdir << " direction" << std::endl;
GparityWilsonImplParams params;
params.twists[gpdir] = 1;
std::vector<int> conj_dirs(Nd,0);
conj_dirs[gpdir] = 1;
ConjugateGimplD::setDirections(conj_dirs);
run_test<GparityDomainWallFermionD, GparityDomainWallFermionF, ConjugateGaugeStatistics>(argc,argv,params);
}else{
std::cout << "Running test with periodic BCs" << std::endl;
WilsonImplParams params;
run_test<DomainWallFermionD, DomainWallFermionF, PeriodicGaugeStatistics>(argc,argv,params);
}
Grid_finalize();
}