mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
e0e42873c1
Adapted GeneralEvenOddRationalRatio and Test_rhmc_EOWilsonRatio_doubleVsMixedPrec to recent changes that require passing in serial RNG For GeneralEvenOddRationalRatio and TwoFlavourEvenOddRatio, broke refresh into two stages, the first of which generates the random field and the second that computes the pseudofermion field. This allows derived classes to override the generation of the random field, for example in testing. Test_dwf_gpforce now uses Gparity in x-direction and APBC in time as opposed to G-parity in time Added Test_action_dwf_gparity2fvs1f that compares the DWF fermion action with the 2f and the 1f (doubled-lattice) implementations of Gparity
329 lines
13 KiB
C++
329 lines
13 KiB
C++
/*************************************************************************************
|
|
|
|
Grid physics library, www.github.com/paboyle/Grid
|
|
|
|
Source file: ./tests/Test_gparity.cc
|
|
|
|
Copyright (C) 2015
|
|
|
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
|
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
|
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 std;
|
|
using namespace Grid;
|
|
;
|
|
|
|
//typedef GparityDomainWallFermionD GparityDiracOp;
|
|
//typedef DomainWallFermionD StandardDiracOp;
|
|
//#define DOP_PARAMS
|
|
|
|
typedef GparityMobiusFermionD GparityDiracOp;
|
|
typedef MobiusFermionD StandardDiracOp;
|
|
#define DOP_PARAMS ,1.5, 0.5
|
|
|
|
|
|
typedef typename GparityDiracOp::FermionField GparityFermionField;
|
|
typedef typename GparityDiracOp::GaugeField GparityGaugeField;
|
|
typedef typename GparityFermionField::vector_type vComplexType;
|
|
|
|
typedef typename StandardDiracOp::FermionField StandardFermionField;
|
|
typedef typename StandardDiracOp::GaugeField StandardGaugeField;
|
|
|
|
enum{ same_vComplex = std::is_same<vComplexType, typename StandardFermionField::vector_type>::value };
|
|
static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIMD complex type");
|
|
|
|
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;
|
|
}
|
|
}
|
|
|
|
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
|
|
std::cout << GridLogMessage<< "* Kernel options --dslash-generic, --dslash-unroll, --dslash-asm" <<std::endl;
|
|
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
|
|
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
|
|
std::cout << GridLogMessage<< "* Testing Gparity Dirac operator "<<std::endl;
|
|
std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplexType::Nsimd()<<std::endl;
|
|
#ifdef GRID_OMP
|
|
if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
|
|
if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential comms compute" <<std::endl;
|
|
#endif
|
|
if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric ) std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
|
|
if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) std::cout << GridLogMessage<< "* Using UNROLLED Nc=3 WilsonKernels" <<std::endl;
|
|
if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3 WilsonKernels" <<std::endl;
|
|
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
|
|
|
|
const int Ls=4;
|
|
//const int L =4;
|
|
//std::vector<int> latt_2f(Nd,L);
|
|
|
|
Coordinate latt_2f = GridDefaultLatt();
|
|
Coordinate latt_1f(latt_2f); latt_1f[nu] = 2*latt_2f[nu];
|
|
int L = latt_2f[nu];
|
|
|
|
|
|
Coordinate simd_layout = GridDefaultSimd(Nd,vComplexType::Nsimd());
|
|
|
|
std::cout << GridLogMessage << "SIMD layout: ";
|
|
for(int i=0;i<simd_layout.size();i++) std::cout << simd_layout[i] << " ";
|
|
std::cout << std::endl;
|
|
|
|
Coordinate mpi_layout = GridDefaultMpi(); //node layout
|
|
|
|
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::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);
|
|
|
|
GparityGaugeField Umu_2f(UGrid_2f);
|
|
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
|
|
|
|
StandardFermionField src (FGrid_2f);
|
|
StandardFermionField tmpsrc(FGrid_2f);
|
|
GparityFermionField src_2f(FGrid_2f);
|
|
StandardFermionField src_1f(FGrid_1f);
|
|
|
|
// Replicate fermion source
|
|
random(RNG5_2f,src);
|
|
PokeIndex<0>(src_2f,src,0);
|
|
tmpsrc=src*2.0;
|
|
PokeIndex<0>(src_2f,tmpsrc,1);
|
|
|
|
StandardFermionField result_1f(FGrid_1f); result_1f=Zero();
|
|
StandardGaugeField Umu_1f(UGrid_1f);
|
|
Replicate(Umu_2f,Umu_1f);
|
|
|
|
//Coordinate grid for reference
|
|
LatticeInteger xcoor_1f(UGrid_1f);
|
|
LatticeCoordinate(xcoor_1f,nu);
|
|
|
|
//Copy-conjugate the gauge field
|
|
//First C-shift the lattice by Lx/2
|
|
{
|
|
StandardGaugeField Umu_shift = conjugate( Cshift(Umu_1f,nu,L) );
|
|
Umu_1f = where( xcoor_1f >= Integer(L), Umu_shift, Umu_1f );
|
|
|
|
// hack test to check the same
|
|
Replicate(Umu_2f,Umu_shift);
|
|
Umu_shift=Umu_shift-Umu_1f;
|
|
cout << GridLogMessage << "Umu diff " << norm2(Umu_shift)<<std::endl;
|
|
|
|
//Make the gauge field antiperiodic in nu-direction
|
|
decltype(PeekIndex<LorentzIndex>(Umu_1f,nu)) Unu(UGrid_1f);
|
|
Unu = PeekIndex<LorentzIndex>(Umu_1f,nu);
|
|
Unu = where(xcoor_1f == Integer(2*L-1), -Unu, Unu);
|
|
PokeIndex<LorentzIndex>(Umu_1f,Unu,nu);
|
|
}
|
|
|
|
//Coordinate grid for reference
|
|
LatticeInteger xcoor_1f5(FGrid_1f);
|
|
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;
|
|
|
|
//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);
|
|
pickCheckerboard(Odd,src_o_1f,src_1f);
|
|
result_o_1f=Zero();
|
|
|
|
SchurDiagMooeeOperator<StandardDiracOp,StandardFermionField> HermOpEO(Ddwf);
|
|
ConjugateGradient<StandardFermionField> CG(1.0e-8,10000);
|
|
CG(HermOpEO,src_o_1f,result_o_1f);
|
|
|
|
//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);
|
|
|
|
for(int disp=-1;disp<=1;disp+=2)
|
|
for(int mu=0;mu<5;mu++)
|
|
{
|
|
GparityFermionField Dsrc_2f(FGrid_2f);
|
|
|
|
StandardFermionField Dsrc_1f(FGrid_1f);
|
|
StandardFermionField Dsrc_2freplica(FGrid_1f);
|
|
StandardFermionField Dsrc_2freplica0(FGrid_1f);
|
|
StandardFermionField Dsrc_2freplica1(FGrid_1f);
|
|
|
|
if ( mu ==0 ) {
|
|
std::cout << GridLogMessage<< " Cross checking entire hopping term"<<std::endl;
|
|
GPDdwf.Dhop(src_2f,Dsrc_2f,DaggerNo);
|
|
Ddwf.Dhop(src_1f,Dsrc_1f,DaggerNo);
|
|
} else {
|
|
std::cout << GridLogMessage<< " Cross checking mu="<<mu<< " disp="<< disp<<std::endl;
|
|
GPDdwf.DhopDir(src_2f,Dsrc_2f,mu,disp);
|
|
Ddwf.DhopDir(src_1f,Dsrc_1f,mu,disp);
|
|
}
|
|
|
|
std::cout << GridLogMessage << "S norms "<< norm2(src_2f) << " " << norm2(src_1f) <<std::endl;
|
|
std::cout << GridLogMessage << "D norms "<< norm2(Dsrc_2f)<< " " << norm2(Dsrc_1f) <<std::endl;
|
|
|
|
StandardFermionField Dsrc_2f0(FGrid_2f); Dsrc_2f0 = PeekIndex<0>(Dsrc_2f,0);
|
|
StandardFermionField Dsrc_2f1(FGrid_2f); Dsrc_2f1 = PeekIndex<0>(Dsrc_2f,1);
|
|
|
|
// Dsrc_2f1 = Dsrc_2f1 - Dsrc_2f0;
|
|
// std::cout << GridLogMessage << " Cross check two halves " <<norm2(Dsrc_2f1)<<std::endl;
|
|
|
|
Replicate(Dsrc_2f0,Dsrc_2freplica0);
|
|
Replicate(Dsrc_2f1,Dsrc_2freplica1);
|
|
|
|
Dsrc_2freplica = where( xcoor_1f5 >= Integer(L), Dsrc_2freplica1,Dsrc_2freplica0 );
|
|
Dsrc_2freplica = Dsrc_2freplica - Dsrc_1f ;
|
|
std::cout << GridLogMessage << " Cross check against doubled latt " <<norm2(Dsrc_2freplica)<<std::endl;
|
|
|
|
// std::cout << Dsrc_2f <<std::endl;
|
|
|
|
}
|
|
|
|
{
|
|
GparityFermionField chi (FGrid_2f); gaussian(RNG5_2f,chi);
|
|
GparityFermionField phi (FGrid_2f); gaussian(RNG5_2f,phi);
|
|
|
|
GparityFermionField chi_e (FrbGrid_2f);
|
|
GparityFermionField chi_o (FrbGrid_2f);
|
|
|
|
GparityFermionField dchi_e (FrbGrid_2f);
|
|
GparityFermionField dchi_o (FrbGrid_2f);
|
|
|
|
GparityFermionField phi_e (FrbGrid_2f);
|
|
GparityFermionField phi_o (FrbGrid_2f);
|
|
|
|
GparityFermionField dphi_e (FrbGrid_2f);
|
|
GparityFermionField dphi_o (FrbGrid_2f);
|
|
|
|
pickCheckerboard(Even,chi_e,chi);
|
|
pickCheckerboard(Odd ,chi_o,chi);
|
|
pickCheckerboard(Even,phi_e,phi);
|
|
pickCheckerboard(Odd ,phi_o,phi);
|
|
|
|
GPDdwf.Meooe(chi_e,dchi_o);
|
|
GPDdwf.Meooe(chi_o,dchi_e);
|
|
GPDdwf.MeooeDag(phi_e,dphi_o);
|
|
GPDdwf.MeooeDag(phi_o,dphi_e);
|
|
|
|
ComplexD pDce = innerProduct(phi_e,dchi_e);
|
|
ComplexD pDco = innerProduct(phi_o,dchi_o);
|
|
ComplexD cDpe = innerProduct(chi_e,dphi_e);
|
|
ComplexD cDpo = innerProduct(chi_o,dphi_o);
|
|
|
|
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
|
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
|
|
|
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
|
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
|
|
|
}
|
|
|
|
GparityFermionField result_2f(FGrid_2f); result_2f=Zero();
|
|
GparityFermionField src_o_2f(FrbGrid_2f);
|
|
GparityFermionField result_o_2f(FrbGrid_2f);
|
|
pickCheckerboard(Odd,src_o_2f,src_2f);
|
|
result_o_2f=Zero();
|
|
|
|
ConjugateGradient<GparityFermionField> CG2f(1.0e-8,10000);
|
|
SchurDiagMooeeOperator<GparityDiracOp,GparityFermionField> HermOpEO2f(GPDdwf);
|
|
CG2f(HermOpEO2f,src_o_2f,result_o_2f);
|
|
|
|
std::cout << "2f cb "<<result_o_2f.Checkerboard()<<std::endl;
|
|
std::cout << "1f cb "<<result_o_1f.Checkerboard()<<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);
|
|
StandardFermionField res1 (FGrid_2f);
|
|
|
|
res0=Zero();
|
|
res1=Zero();
|
|
|
|
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;
|
|
|
|
//poke odd onto non-cb field
|
|
setCheckerboard(res0,res0o);
|
|
setCheckerboard(res1,res1o);
|
|
|
|
StandardFermionField replica (FGrid_1f);
|
|
StandardFermionField replica0(FGrid_1f);
|
|
StandardFermionField replica1(FGrid_1f);
|
|
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 1f reconstructed from 2f: " <<norm2(replica)<<" Actual 1f: "<< norm2(replica0)<<std::endl;
|
|
|
|
replica = replica - replica0;
|
|
|
|
std::cout << "Norm2 of difference in solutions is " <<norm2(replica)<<std::endl;
|
|
|
|
|
|
Grid_finalize();
|
|
}
|