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

Merge branch 'feature/dirichlet-gparity' into feature/dirichlet

Browse Source
This commit is contained in:
Peter Boyle
2022-06-15 19:23:48 -04:00
parent 31efa5c4da 3d8146b596
commit 8208a6214f
30 changed files with 2142 additions and 188 deletions
Download Patch File Download Diff File Expand all files Collapse all files

187
tests/core/Test_fft_gfix.cc
View File

@ -29,14 +29,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/Grid.h>
using namespace Grid;
;
int main (int argc, char ** argv)
{
template<typename Gimpl>
void run(double alpha, bool do_fft_gfix){
std::vector<int> seeds({1,2,3,4});
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
Coordinate latt_size = GridDefaultLatt();
@ -55,10 +51,7 @@ int main (int argc, char ** argv)
FFT theFFT(&GRID);
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<<GridLogMessage << "Using alpha=" << alpha << std::endl;
// int coulomb_dir = -1;
int coulomb_dir = Nd-1;
@ -72,81 +65,167 @@ int main (int argc, char ** argv)
LatticeColourMatrix xform1(&GRID); // Gauge xform
LatticeColourMatrix xform2(&GRID); // Gauge xform
LatticeColourMatrix xform3(&GRID); // Gauge xform
//#########################################################################################
std::cout<< "*********************************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to Landau gauge with randomly transformed unit gauge configuration *" <<std::endl;
std::cout<< "*********************************************************************************************************" <<std::endl;
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
Uorg=Umu;
Real init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
//Apply a random gauge transformation to the unit gauge config
Urnd=Umu;
SU<Nc>::RandomGaugeTransform<Gimpl>(pRNG,Urnd,g);
SU<Nc>::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
Real alpha=0.1;
//Gauge fix the randomly transformed field
Umu = Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
// Check the gauge xform matrices
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform1);
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform1);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
Real plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
Uorg = Uorg - Umu;
std::cout << " Norm Difference "<< norm2(Uorg) << std::endl;
std::cout << " Norm "<< norm2(Umu) << std::endl;
std::cout << " Norm difference between a unit gauge configuration and the gauge fixed configuration "<< norm2(Uorg) << " (expect 0)" << std::endl;
std::cout << " Norm of gauge fixed configuration "<< norm2(Umu) << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "*************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with unit gauge configuration *" <<std::endl;
std::cout<< "*************************************************************************************" <<std::endl;
Umu=Urnd;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
Utmp=Urnd;
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform2);
Utmp = Utmp - Umu;
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
Umu=Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
//#########################################################################################
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform2);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to Landau gauge with random configuration **" <<std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing non-unit configuration *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "******************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with random configuration **" <<std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
//#########################################################################################
std::cout<< "*******************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to coulomb gauge with random configuration *" <<std::endl;
std::cout<< "*******************************************************************************************" <<std::endl;
Umu=Urnd;
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,false,coulomb_dir);
std::cout << Umu<<std::endl;
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "*******************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge with random configuration *" <<std::endl;
std::cout<< "*******************************************************************************************" <<std::endl;
Umu=Urnd;
SU<Nc>::HotConfiguration(pRNG,Umu);
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
double alpha=0.1; //step size
std::string gimpl = "periodic";
bool do_fft_gfix = true; //test fourier transformed gfix as well as steepest descent
for(int i=1;i<argc;i++){
std::string sarg(argv[i]);
if(sarg == "--gimpl"){
assert(i<argc-1 && "--gimpl option requires an argument");
gimpl = argv[i+1];
if(gimpl != "periodic" && gimpl != "conjugate")
assert(0 && "Invalid gimpl");
if(gimpl == "conjugate")
alpha = 0.025; //default alpha too large for CCBC
}else if(sarg == "--no-fft-gfix"){
std::cout << "Not doing the Fourier accelerated gauge fixing tests" << std::endl;
do_fft_gfix = false;
}else if(sarg == "--alpha"){
assert(i<argc-1 && "--alpha option requires an argument");
std::istringstream ss(argv[i+1]); ss >> alpha;
}
}
if(gimpl == "periodic"){
std::cout << GridLogMessage << "Using periodic boundary condition" << std::endl;
run<PeriodicGimplR>(alpha, do_fft_gfix);
}else{
std::vector<int> conjdirs = {1,1,0,0}; //test with 2 conjugate dirs and 2 not
std::cout << GridLogMessage << "Using complex conjugate boundary conditions in dimensions ";
for(int i=0;i<Nd;i++)
if(conjdirs[i])
std::cout << i << " ";
std::cout << std::endl;
ConjugateGimplR::setDirections(conjdirs);
run<ConjugateGimplR>(alpha, do_fft_gfix);
}
Grid_finalize();
}

60
tests/core/Test_gamma.cc
View File

@ -228,6 +228,59 @@ void checkGammaL(const Gamma::Algebra a, GridSerialRNG &rng)
std::cout << std::endl;
}
void checkChargeConjMatrix(){
//Check the properties of the charge conjugation matrix
//In the Grid basis C = -\gamma^2 \gamma^4
SpinMatrix C = testAlgebra[Gamma::Algebra::MinusGammaY] * testAlgebra[Gamma::Algebra::GammaT];
SpinMatrix mC = -C;
SpinMatrix one = testAlgebra[Gamma::Algebra::Identity];
std::cout << "Testing properties of charge conjugation matrix C = -\\gamma^2 \\gamma^4 (in Grid's basis)" << std::endl;
//C^T = -C
SpinMatrix Ct = transpose(C);
std::cout << GridLogMessage << "C^T=-C ";
test(Ct, mC);
std::cout << std::endl;
//C^\dagger = -C
SpinMatrix Cdag = adj(C);
std::cout << GridLogMessage << "C^dag=-C ";
test(Cdag, mC);
std::cout << std::endl;
//C^* = C
SpinMatrix Cstar = conjugate(C);
std::cout << GridLogMessage << "C^*=C ";
test(Cstar, C);
std::cout << std::endl;
//C^{-1} = -C
SpinMatrix CinvC = mC * C;
std::cout << GridLogMessage << "C^{-1}=-C ";
test(CinvC, one);
std::cout << std::endl;
// C^{-1} \gamma^\mu C = -[\gamma^\mu]^T
Gamma::Algebra gmu_a[4] = { Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ, Gamma::Algebra::GammaT };
for(int mu=0;mu<4;mu++){
SpinMatrix gmu = testAlgebra[gmu_a[mu]];
SpinMatrix Cinv_gmu_C = mC * gmu * C;
SpinMatrix mgmu_T = -transpose(gmu);
std::cout << GridLogMessage << "C^{-1} \\gamma^" << mu << " C = -[\\gamma^" << mu << "]^T ";
test(Cinv_gmu_C, mgmu_T);
std::cout << std::endl;
}
//[C, \gamma^5] = 0
SpinMatrix Cg5 = C * testAlgebra[Gamma::Algebra::Gamma5];
SpinMatrix g5C = testAlgebra[Gamma::Algebra::Gamma5] * C;
std::cout << GridLogMessage << "C \\gamma^5 = \\gamma^5 C";
test(Cg5, g5C);
std::cout << std::endl;
}
int main(int argc, char *argv[])
{
Grid_init(&argc,&argv);
@ -270,6 +323,13 @@ int main(int argc, char *argv[])
{
checkGammaL(i, sRNG);
}
std::cout << GridLogMessage << "======== Charge conjugation matrix check" << std::endl;
checkChargeConjMatrix();
std::cout << GridLogMessage << std::endl;
Grid_finalize();

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;
}

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);

6
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);

8
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);

81
tests/lanczos/Test_dwf_lanczos.cc
View File

@ -31,14 +31,38 @@ using namespace std;
using namespace Grid;
;
typedef typename GparityDomainWallFermionR::FermionField FermionField;
template<typename Action>
struct Setup{};
RealD AllZero(RealD x){ return 0.;}
template<>
struct Setup<GparityMobiusFermionR>{
static GparityMobiusFermionR* getAction(LatticeGaugeField &Umu,
GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){
RealD mass=0.01;
RealD M5=1.8;
RealD mob_b=1.5;
GparityMobiusFermionD ::ImplParams params;
std::vector<int> twists({1,1,1,0});
params.twists = twists;
return new GparityMobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
template<>
struct Setup<DomainWallFermionR>{
static DomainWallFermionR* getAction(LatticeGaugeField &Umu,
GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){
RealD mass=0.01;
RealD M5=1.8;
return new DomainWallFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
}
};
template<typename Action>
void run(){
typedef typename Action::FermionField FermionField;
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@ -56,24 +80,10 @@ int main (int argc, char ** argv)
LatticeGaugeField Umu(UGrid);
SU<Nc>::HotConfiguration(RNG4, Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
RealD mass=0.01;
RealD M5=1.8;
RealD mob_b=1.5;
// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
GparityMobiusFermionD ::ImplParams params;
std::vector<int> twists({1,1,1,0});
params.twists = twists;
GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
// MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
// SchurDiagTwoOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
// SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
Action *action = Setup<Action>::getAction(Umu,FGrid,FrbGrid,UGrid,UrbGrid);
//MdagMLinearOperator<Action,FermionField> HermOp(Ddwf);
SchurDiagTwoOperator<Action,FermionField> HermOp(*action);
const int Nstop = 30;
const int Nk = 40;
@ -90,8 +100,7 @@ int main (int argc, char ** argv)
PlainHermOp<FermionField> Op (HermOp);
ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby,Op,Nstop,Nk,Nm,resid,MaxIt);
std::vector<RealD> eval(Nm);
FermionField src(FrbGrid);
gaussian(RNG5rb,src);
@ -103,6 +112,28 @@ int main (int argc, char ** argv)
int Nconv;
IRL.calc(eval,evec,src,Nconv);
delete action;
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::string action = "GparityMobius";
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "-action"){
action = argv[i+1];
}
}
if(action == "GparityMobius"){
run<GparityMobiusFermionR>();
}else if(action == "DWF"){
run<DomainWallFermionR>();
}else{
std::cout << "Unknown action" << std::endl;
exit(1);
}
Grid_finalize();
}

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();
}