mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 15:55:37 +00:00
460 lines
16 KiB
C++
460 lines
16 KiB
C++
/*************************************************************************************
|
|
|
|
grid` physics library, www.github.com/paboyle/Grid
|
|
|
|
Source file: ./tests/Test_cshift.cc
|
|
|
|
Copyright (C) 2015
|
|
|
|
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
|
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;
|
|
;
|
|
|
|
int main (int argc, char ** argv)
|
|
{
|
|
Grid_init(&argc,&argv);
|
|
|
|
int threads = GridThread::GetThreads();
|
|
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
|
|
|
Coordinate latt_size = GridDefaultLatt();
|
|
Coordinate simd_layout( { vComplexD::Nsimd(),1,1,1});
|
|
Coordinate mpi_layout = GridDefaultMpi();
|
|
|
|
int vol = 1;
|
|
for(int d=0;d<latt_size.size();d++){
|
|
vol = vol * latt_size[d];
|
|
}
|
|
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
|
|
GridRedBlackCartesian RBGRID(&GRID);
|
|
|
|
LatticeComplexD one(&GRID);
|
|
LatticeComplexD zz(&GRID);
|
|
LatticeComplexD C(&GRID);
|
|
LatticeComplexD Ctilde(&GRID);
|
|
LatticeComplexD Cref (&GRID);
|
|
LatticeComplexD Csav (&GRID);
|
|
LatticeComplexD coor(&GRID);
|
|
|
|
LatticeSpinMatrixD S(&GRID);
|
|
LatticeSpinMatrixD Stilde(&GRID);
|
|
|
|
Coordinate p({1,3,2,3});
|
|
|
|
one = ComplexD(1.0,0.0);
|
|
zz = ComplexD(0.0,0.0);
|
|
|
|
ComplexD ci(0.0,1.0);
|
|
|
|
|
|
std::cout<<"*************************************************"<<std::endl;
|
|
std::cout<<"Testing Fourier from of known plane wave "<<std::endl;
|
|
std::cout<<"*************************************************"<<std::endl;
|
|
C=Zero();
|
|
for(int mu=0;mu<4;mu++){
|
|
RealD TwoPiL = M_PI * 2.0/ latt_size[mu];
|
|
LatticeCoordinate(coor,mu);
|
|
C = C + (TwoPiL * p[mu]) * coor;
|
|
}
|
|
|
|
C = exp(C*ci);
|
|
Csav = C;
|
|
S=Zero();
|
|
S = S+C;
|
|
|
|
FFT theFFT(&GRID);
|
|
|
|
Ctilde=C;
|
|
std::cout<<" Benchmarking FFT of LatticeComplex "<<std::endl;
|
|
theFFT.FFT_dim(Ctilde,Ctilde,0,FFT::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
|
|
theFFT.FFT_dim(Ctilde,Ctilde,1,FFT::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
|
|
theFFT.FFT_dim(Ctilde,Ctilde,2,FFT::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
|
|
theFFT.FFT_dim(Ctilde,Ctilde,3,FFT::forward); std::cout << theFFT.MFlops()<<" Mflops "<<std::endl;
|
|
|
|
// C=Zero();
|
|
// Ctilde = where(abs(Ctilde)<1.0e-10,C,Ctilde);
|
|
TComplexD cVol;
|
|
cVol()()() = vol;
|
|
|
|
Cref=Zero();
|
|
pokeSite(cVol,Cref,p);
|
|
// std::cout <<"Ctilde "<< Ctilde <<std::endl;
|
|
// std::cout <<"Cref "<< Cref <<std::endl;
|
|
|
|
Cref=Cref-Ctilde;
|
|
std::cout << "diff scalar "<<norm2(Cref) << std::endl;
|
|
C=Csav;
|
|
theFFT.FFT_all_dim(Ctilde,C,FFT::forward);
|
|
theFFT.FFT_all_dim(Cref,Ctilde,FFT::backward);
|
|
|
|
std::cout << norm2(C) << " " << norm2(Ctilde) << " " << norm2(Cref)<< " vol " << vol<< std::endl;
|
|
|
|
Cref= Cref - C;
|
|
std::cout << " invertible check " << norm2(Cref)<<std::endl;
|
|
|
|
Stilde=S;
|
|
std::cout<<" Benchmarking FFT of LatticeSpinMatrix "<<std::endl;
|
|
theFFT.FFT_dim(Stilde,S,0,FFT::forward); std::cout << theFFT.MFlops()<<" mflops "<<std::endl;
|
|
theFFT.FFT_dim(Stilde,S,1,FFT::forward); std::cout << theFFT.MFlops()<<" mflops "<<std::endl;
|
|
theFFT.FFT_dim(Stilde,S,2,FFT::forward); std::cout << theFFT.MFlops()<<" mflops "<<std::endl;
|
|
theFFT.FFT_dim(Stilde,S,3,FFT::forward); std::cout << theFFT.MFlops()<<" mflops "<<std::endl;
|
|
|
|
SpinMatrixD Sp;
|
|
Sp = Zero(); Sp = Sp+cVol;
|
|
|
|
S=Zero();
|
|
pokeSite(Sp,S,p);
|
|
|
|
S= S-Stilde;
|
|
std::cout << "diff FT[SpinMat] "<<norm2(S) << std::endl;
|
|
|
|
/*
|
|
*/
|
|
std::vector<int> seeds({1,2,3,4});
|
|
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); // naughty seeding
|
|
GridParallelRNG pRNG(&GRID);
|
|
pRNG.SeedFixedIntegers(seeds);
|
|
|
|
LatticeGaugeFieldD Umu(&GRID);
|
|
|
|
SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
|
|
// Umu=Zero();
|
|
////////////////////////////////////////////////////
|
|
// Wilson test
|
|
////////////////////////////////////////////////////
|
|
{
|
|
LatticeFermionD src(&GRID); gaussian(pRNG,src);
|
|
LatticeFermionD tmp(&GRID);
|
|
LatticeFermionD ref(&GRID);
|
|
|
|
RealD mass=0.01;
|
|
WilsonFermionD Dw(Umu,GRID,RBGRID,mass);
|
|
|
|
Dw.M(src,tmp);
|
|
|
|
std::cout << "Dw src = " <<norm2(src)<<std::endl;
|
|
std::cout << "Dw tmp = " <<norm2(tmp)<<std::endl;
|
|
|
|
Dw.FreePropagator(tmp,ref,mass);
|
|
|
|
std::cout << "Dw ref = " <<norm2(ref)<<std::endl;
|
|
|
|
ref = ref - src;
|
|
|
|
std::cout << "Dw ref-src = " <<norm2(ref)<<std::endl;
|
|
}
|
|
|
|
////////////////////////////////////////////////////
|
|
// Dwf matrix
|
|
////////////////////////////////////////////////////
|
|
{
|
|
std::cout<<"****************************************"<<std::endl;
|
|
std::cout<<"Testing Fourier representation of Ddwf"<<std::endl;
|
|
std::cout<<"****************************************"<<std::endl;
|
|
|
|
const int Ls=16;
|
|
const int sdir=0;
|
|
RealD mass=0.01;
|
|
RealD M5 =1.0;
|
|
Gamma G5(Gamma::Algebra::Gamma5);
|
|
|
|
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
|
|
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
|
|
|
|
std::cout<<"Making Ddwf"<<std::endl;
|
|
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5);
|
|
|
|
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds);
|
|
LatticeFermionD src5(FGrid); gaussian(RNG5,src5);
|
|
LatticeFermionD src5_p(FGrid);
|
|
LatticeFermionD result5(FGrid);
|
|
LatticeFermionD ref5(FGrid);
|
|
LatticeFermionD tmp5(FGrid);
|
|
|
|
/////////////////////////////////////////////////////////////////
|
|
// result5 is the non pert operator in 4d mom space
|
|
/////////////////////////////////////////////////////////////////
|
|
Ddwf.M(src5,tmp5);
|
|
ref5 = tmp5;
|
|
|
|
FFT theFFT5(FGrid);
|
|
|
|
theFFT5.FFT_dim(result5,tmp5,1,FFT::forward); tmp5 = result5;
|
|
theFFT5.FFT_dim(result5,tmp5,2,FFT::forward); tmp5 = result5;
|
|
theFFT5.FFT_dim(result5,tmp5,3,FFT::forward); tmp5 = result5;
|
|
theFFT5.FFT_dim(result5,tmp5,4,FFT::forward); result5 = result5*ComplexD(::sqrt(1.0/vol),0.0);
|
|
|
|
std::cout<<"Fourier xformed Ddwf"<<std::endl;
|
|
|
|
tmp5 = src5;
|
|
theFFT5.FFT_dim(src5_p,tmp5,1,FFT::forward); tmp5 = src5_p;
|
|
theFFT5.FFT_dim(src5_p,tmp5,2,FFT::forward); tmp5 = src5_p;
|
|
theFFT5.FFT_dim(src5_p,tmp5,3,FFT::forward); tmp5 = src5_p;
|
|
theFFT5.FFT_dim(src5_p,tmp5,4,FFT::forward); src5_p = src5_p*ComplexD(::sqrt(1.0/vol),0.0);
|
|
|
|
std::cout<<"Fourier xformed src5"<<std::endl;
|
|
|
|
/////////////////////////////////////////////////////////////////
|
|
// work out the predicted from Fourier
|
|
/////////////////////////////////////////////////////////////////
|
|
Gamma::Algebra Gmu [] = {
|
|
Gamma::Algebra::GammaX,
|
|
Gamma::Algebra::GammaY,
|
|
Gamma::Algebra::GammaZ,
|
|
Gamma::Algebra::GammaT,
|
|
Gamma::Algebra::Gamma5
|
|
};
|
|
LatticeFermionD Kinetic(FGrid); Kinetic = Zero();
|
|
LatticeComplexD kmu(FGrid);
|
|
LatticeInteger scoor(FGrid);
|
|
LatticeComplexD sk (FGrid); sk = Zero();
|
|
LatticeComplexD sk2(FGrid); sk2= Zero();
|
|
LatticeComplexD W(FGrid); W= Zero();
|
|
// LatticeComplexD a(FGrid); a= Zero();
|
|
LatticeComplexD one(FGrid); one =ComplexD(1.0,0.0);
|
|
ComplexD ci(0.0,1.0);
|
|
|
|
for(int mu=0;mu<Nd;mu++) {
|
|
|
|
LatticeCoordinate(kmu,mu+1);
|
|
|
|
RealD TwoPiL = M_PI * 2.0/ latt_size[mu];
|
|
|
|
kmu = TwoPiL * kmu;
|
|
|
|
sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5);
|
|
sk = sk + sin(kmu) *sin(kmu);
|
|
|
|
// -1/2 Dw -> 1/2 gmu (eip - emip) = i sinp gmu
|
|
Kinetic = Kinetic + sin(kmu)*ci*(Gamma(Gmu[mu])*src5_p);
|
|
|
|
}
|
|
|
|
// NB implicit sum over mu
|
|
//
|
|
// 1-1/2 Dw = 1 - 1/2 ( eip+emip)
|
|
// = - 1/2 (ei - 2 + emi)
|
|
// = - 1/4 2 (eih - eimh)(eih - eimh)
|
|
// = 2 sink/2 ink/2 = sk2
|
|
|
|
W = one - M5 + sk2;
|
|
Kinetic = Kinetic + W * src5_p;
|
|
|
|
LatticeCoordinate(scoor,sdir);
|
|
|
|
tmp5 = Cshift(src5_p,sdir,+1);
|
|
tmp5 = (tmp5 - G5*tmp5)*0.5;
|
|
tmp5 = where(scoor==Integer(Ls-1),mass*tmp5,-tmp5);
|
|
Kinetic = Kinetic + tmp5;
|
|
|
|
tmp5 = Cshift(src5_p,sdir,-1);
|
|
tmp5 = (tmp5 + G5*tmp5)*0.5;
|
|
tmp5 = where(scoor==Integer(0),mass*tmp5,-tmp5);
|
|
Kinetic = Kinetic + tmp5;
|
|
|
|
std::cout<<"Momentum space Ddwf "<< norm2(Kinetic)<<std::endl;
|
|
std::cout<<"Stencil Ddwf "<< norm2(result5)<<std::endl;
|
|
|
|
result5 = result5 - Kinetic;
|
|
std::cout<<"diff "<< norm2(result5)<<std::endl;
|
|
|
|
}
|
|
|
|
////////////////////////////////////////////////////
|
|
// Dwf prop
|
|
////////////////////////////////////////////////////
|
|
{
|
|
std::cout<<"****************************************"<<std::endl;
|
|
std::cout << "Testing Ddwf Ht Mom space 4d propagator \n";
|
|
std::cout<<"****************************************"<<std::endl;
|
|
|
|
LatticeFermionD src(&GRID); gaussian(pRNG,src);
|
|
LatticeFermionD tmp(&GRID);
|
|
LatticeFermionD ref(&GRID);
|
|
LatticeFermionD diff(&GRID);
|
|
|
|
Coordinate point(4,0);
|
|
src=Zero();
|
|
SpinColourVectorD ferm; gaussian(sRNG,ferm);
|
|
pokeSite(ferm,src,point);
|
|
|
|
const int Ls=32;
|
|
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
|
|
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
|
|
|
|
RealD mass=0.01;
|
|
RealD M5 =0.8;
|
|
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5);
|
|
|
|
// Momentum space prop
|
|
std::cout << " Solving by FFT and Feynman rules" <<std::endl;
|
|
Ddwf.FreePropagator(src,ref,mass) ;
|
|
|
|
Gamma G5(Gamma::Algebra::Gamma5);
|
|
|
|
LatticeFermionD src5(FGrid); src5=Zero();
|
|
LatticeFermionD tmp5(FGrid);
|
|
LatticeFermionD result5(FGrid); result5=Zero();
|
|
LatticeFermionD result4(&GRID);
|
|
const int sdir=0;
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Domain wall physical field source
|
|
////////////////////////////////////////////////////////////////////////
|
|
/*
|
|
chi_5[0] = chiralProjectPlus(chi);
|
|
chi_5[Ls-1]= chiralProjectMinus(chi);
|
|
*/
|
|
tmp = (src + G5*src)*0.5; InsertSlice(tmp,src5, 0,sdir);
|
|
tmp = (src - G5*src)*0.5; InsertSlice(tmp,src5,Ls-1,sdir);
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Conjugate gradient on normal equations system
|
|
////////////////////////////////////////////////////////////////////////
|
|
std::cout << " Solving by Conjugate Gradient (CGNE)" <<std::endl;
|
|
Ddwf.Mdag(src5,tmp5);
|
|
src5=tmp5;
|
|
MdagMLinearOperator<DomainWallFermionD,LatticeFermionD> HermOp(Ddwf);
|
|
ConjugateGradient<LatticeFermionD> CG(1.0e-16,10000);
|
|
CG(HermOp,src5,result5);
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Domain wall physical field propagator
|
|
////////////////////////////////////////////////////////////////////////
|
|
/*
|
|
psi = chiralProjectMinus(psi_5[0]);
|
|
psi += chiralProjectPlus(psi_5[Ls-1]);
|
|
*/
|
|
ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5;
|
|
ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5;
|
|
|
|
std::cout << " Taking difference" <<std::endl;
|
|
std::cout << "Ddwf result4 "<<norm2(result4)<<std::endl;
|
|
std::cout << "Ddwf ref "<<norm2(ref)<<std::endl;
|
|
|
|
diff = ref - result4;
|
|
std::cout << "result - ref "<<norm2(diff)<<std::endl;
|
|
|
|
}
|
|
|
|
|
|
|
|
////////////////////////////////////////////////////
|
|
// Dwf prop
|
|
////////////////////////////////////////////////////
|
|
{
|
|
std::cout<<"****************************************"<<std::endl;
|
|
std::cout << "Testing Dov Ht Mom space 4d propagator \n";
|
|
std::cout<<"****************************************"<<std::endl;
|
|
|
|
LatticeFermionD src(&GRID); gaussian(pRNG,src);
|
|
LatticeFermionD tmp(&GRID);
|
|
LatticeFermionD ref(&GRID);
|
|
LatticeFermionD diff(&GRID);
|
|
|
|
Coordinate point(4,0);
|
|
src=Zero();
|
|
SpinColourVectorD ferm; gaussian(sRNG,ferm);
|
|
pokeSite(ferm,src,point);
|
|
|
|
const int Ls=48;
|
|
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
|
|
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
|
|
|
|
RealD mass=0.01;
|
|
RealD M5 =0.8;
|
|
|
|
OverlapWilsonCayleyTanhFermionD Dov(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5,1.0);
|
|
|
|
// Momentum space prop
|
|
std::cout << " Solving by FFT and Feynman rules" <<std::endl;
|
|
Dov.FreePropagator(src,ref,mass) ;
|
|
|
|
Gamma G5(Gamma::Algebra::Gamma5);
|
|
|
|
LatticeFermionD src5(FGrid); src5=Zero();
|
|
LatticeFermionD tmp5(FGrid);
|
|
LatticeFermionD result5(FGrid); result5=Zero();
|
|
LatticeFermionD result4(&GRID);
|
|
const int sdir=0;
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Domain wall physical field source; need D_minus
|
|
////////////////////////////////////////////////////////////////////////
|
|
/*
|
|
chi_5[0] = chiralProjectPlus(chi);
|
|
chi_5[Ls-1]= chiralProjectMinus(chi);
|
|
*/
|
|
tmp = (src + G5*src)*0.5; InsertSlice(tmp,src5, 0,sdir);
|
|
tmp = (src - G5*src)*0.5; InsertSlice(tmp,src5,Ls-1,sdir);
|
|
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Conjugate gradient on normal equations system
|
|
////////////////////////////////////////////////////////////////////////
|
|
std::cout << " Solving by Conjugate Gradient (CGNE)" <<std::endl;
|
|
Dov.Dminus(src5,tmp5);
|
|
src5=tmp5;
|
|
Dov.Mdag(src5,tmp5);
|
|
src5=tmp5;
|
|
MdagMLinearOperator<OverlapWilsonCayleyTanhFermionD,LatticeFermionD> HermOp(Dov);
|
|
ConjugateGradient<LatticeFermionD> CG(1.0e-16,10000);
|
|
CG(HermOp,src5,result5);
|
|
|
|
////////////////////////////////////////////////////////////////////////
|
|
// Domain wall physical field propagator
|
|
////////////////////////////////////////////////////////////////////////
|
|
/*
|
|
psi = chiralProjectMinus(psi_5[0]);
|
|
psi += chiralProjectPlus(psi_5[Ls-1]);
|
|
*/
|
|
ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5;
|
|
ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5;
|
|
|
|
std::cout << " Taking difference" <<std::endl;
|
|
std::cout << "Dov result4 "<<norm2(result4)<<std::endl;
|
|
std::cout << "Dov ref "<<norm2(ref)<<std::endl;
|
|
|
|
diff = ref - result4;
|
|
std::cout << "result - ref "<<norm2(diff)<<std::endl;
|
|
|
|
}
|
|
|
|
{
|
|
/*
|
|
*
|
|
typedef GaugeImplTypes<vComplexD, 1> QEDGimplTypesD;
|
|
typedef Photon<QEDGimplTypesD> QEDGaction;
|
|
|
|
QEDGaction Maxwell(QEDGaction::FEYNMAN_L);
|
|
QEDGaction::GaugeField Prop(&GRID);Prop=Zero();
|
|
QEDGaction::GaugeField Source(&GRID);Source=Zero();
|
|
|
|
Maxwell.FreePropagator (Source,Prop);
|
|
std::cout << " MaxwellFree propagator\n";
|
|
*/
|
|
}
|
|
Grid_finalize();
|
|
}
|