mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-18 15:57:05 +01:00
Merge branch 'develop' into feature/hmc_generalise
This commit is contained in:
Guido Cossu
@ -1,6 +1,6 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
grid` physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_cshift.cc
|
||||
|
||||
@ -46,57 +46,79 @@ int main (int argc, char ** argv)
|
||||
for(int d=0;d<latt_size.size();d++){
|
||||
vol = vol * latt_size[d];
|
||||
}
|
||||
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
|
||||
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
|
||||
GridRedBlackCartesian RBGRID(latt_size,simd_layout,mpi_layout);
|
||||
|
||||
LatticeComplexD one(&Fine);
|
||||
LatticeComplexD zz(&Fine);
|
||||
LatticeComplexD C(&Fine);
|
||||
LatticeComplexD Ctilde(&Fine);
|
||||
LatticeComplexD coor(&Fine);
|
||||
LatticeComplexD one(&GRID);
|
||||
LatticeComplexD zz(&GRID);
|
||||
LatticeComplexD C(&GRID);
|
||||
LatticeComplexD Ctilde(&GRID);
|
||||
LatticeComplexD Cref (&GRID);
|
||||
LatticeComplexD Csav (&GRID);
|
||||
LatticeComplexD coor(&GRID);
|
||||
|
||||
LatticeSpinMatrixD S(&Fine);
|
||||
LatticeSpinMatrixD Stilde(&Fine);
|
||||
LatticeSpinMatrixD S(&GRID);
|
||||
LatticeSpinMatrixD Stilde(&GRID);
|
||||
|
||||
std::vector<int> p({1,2,3,2});
|
||||
std::vector<int> 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 = C + (TwoPiL * p[mu]) * coor;
|
||||
}
|
||||
|
||||
C = exp(C*ci);
|
||||
|
||||
Csav = C;
|
||||
S=zero;
|
||||
S = S+C;
|
||||
|
||||
FFT theFFT(&Fine);
|
||||
FFT theFFT(&GRID);
|
||||
|
||||
theFFT.FFT_dim(Ctilde,C,0,FFT::forward); C=Ctilde; std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Ctilde,C,1,FFT::forward); C=Ctilde; std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Ctilde,C,2,FFT::forward); C=Ctilde; std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Ctilde,C,3,FFT::forward); std::cout << theFFT.MFlops()<<std::endl;
|
||||
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;
|
||||
|
||||
C=zero;
|
||||
pokeSite(cVol,C,p);
|
||||
C=C-Ctilde;
|
||||
std::cout << "diff scalar "<<norm2(C) << std::endl;
|
||||
Cref=zero;
|
||||
pokeSite(cVol,Cref,p);
|
||||
// std::cout <<"Ctilde "<< Ctilde <<std::endl;
|
||||
// std::cout <<"Cref "<< Cref <<std::endl;
|
||||
|
||||
theFFT.FFT_dim(Stilde,S,0,FFT::forward); S=Stilde; std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,1,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,2,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,3,FFT::forward);std::cout << theFFT.MFlops()<<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;
|
||||
@ -107,5 +129,331 @@ int main (int argc, char ** argv)
|
||||
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::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::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT,
|
||||
Gamma::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);
|
||||
|
||||
std::vector<int> 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::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);
|
||||
|
||||
std::vector<int> 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::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();
|
||||
}
|
||||
|
||||
300
tests/core/Test_fft_gfix.cc
Normal file
300
tests/core/Test_fft_gfix.cc
Normal file
@ -0,0 +1,300 @@
|
||||
/*************************************************************************************
|
||||
|
||||
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;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
template <class Gimpl>
|
||||
class FourierAcceleratedGaugeFixer : public Gimpl {
|
||||
public:
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
typedef typename Gimpl::GaugeLinkField GaugeMat;
|
||||
typedef typename Gimpl::GaugeField GaugeLorentz;
|
||||
|
||||
static void GaugeLinkToLieAlgebraField(const std::vector<GaugeMat> &U,std::vector<GaugeMat> &A) {
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
// ImplComplex cmi(0.0,-1.0);
|
||||
Complex cmi(0.0,-1.0);
|
||||
A[mu] = Ta(U[mu]) * cmi;
|
||||
}
|
||||
}
|
||||
static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu) {
|
||||
dmuAmu=zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
|
||||
}
|
||||
}
|
||||
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol) {
|
||||
GridBase *grid = Umu._grid;
|
||||
|
||||
Real org_plaq =WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
Real org_link_trace=WilsonLoops<Gimpl>::linkTrace(Umu);
|
||||
Real old_trace = org_link_trace;
|
||||
Real trG;
|
||||
|
||||
std::vector<GaugeMat> U(Nd,grid);
|
||||
GaugeMat dmuAmu(grid);
|
||||
|
||||
for(int i=0;i<maxiter;i++){
|
||||
for(int mu=0;mu<Nd;mu++) U[mu]= PeekIndex<LorentzIndex>(Umu,mu);
|
||||
//trG = SteepestDescentStep(U,alpha,dmuAmu);
|
||||
trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu);
|
||||
for(int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(Umu,U[mu],mu);
|
||||
// Monitor progress and convergence test
|
||||
// infrequently to minimise cost overhead
|
||||
if ( i %20 == 0 ) {
|
||||
Real plaq =WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
Real link_trace=WilsonLoops<Gimpl>::linkTrace(Umu);
|
||||
|
||||
std::cout << GridLogMessage << " Iteration "<<i<< " plaq= "<<plaq<< " dmuAmu " << norm2(dmuAmu)<< std::endl;
|
||||
|
||||
Real Phi = 1.0 - old_trace / link_trace ;
|
||||
Real Omega= 1.0 - trG;
|
||||
|
||||
|
||||
std::cout << GridLogMessage << " Iteration "<<i<< " Phi= "<<Phi<< " Omega= " << Omega<< " trG " << trG <<std::endl;
|
||||
if ( (Omega < Omega_tol) && ( ::fabs(Phi) < Phi_tol) ) {
|
||||
std::cout << GridLogMessage << "Converged ! "<<std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
old_trace = link_trace;
|
||||
|
||||
}
|
||||
}
|
||||
};
|
||||
static Real SteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
|
||||
GridBase *grid = U[0]._grid;
|
||||
|
||||
std::vector<GaugeMat> A(Nd,grid);
|
||||
GaugeMat g(grid);
|
||||
|
||||
GaugeLinkToLieAlgebraField(U,A);
|
||||
ExpiAlphaDmuAmu(A,g,alpha,dmuAmu);
|
||||
|
||||
|
||||
Real vol = grid->gSites();
|
||||
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
|
||||
|
||||
SU<Nc>::GaugeTransform(U,g);
|
||||
|
||||
return trG;
|
||||
}
|
||||
|
||||
static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
|
||||
|
||||
GridBase *grid = U[0]._grid;
|
||||
|
||||
Real vol = grid->gSites();
|
||||
|
||||
FFT theFFT((GridCartesian *)grid);
|
||||
|
||||
LatticeComplex Fp(grid);
|
||||
LatticeComplex psq(grid); psq=zero;
|
||||
LatticeComplex pmu(grid);
|
||||
LatticeComplex one(grid); one = Complex(1.0,0.0);
|
||||
|
||||
GaugeMat g(grid);
|
||||
GaugeMat dmuAmu_p(grid);
|
||||
std::vector<GaugeMat> A(Nd,grid);
|
||||
|
||||
GaugeLinkToLieAlgebraField(U,A);
|
||||
|
||||
DmuAmu(A,dmuAmu);
|
||||
|
||||
theFFT.FFT_all_dim(dmuAmu_p,dmuAmu,FFT::forward);
|
||||
|
||||
//////////////////////////////////
|
||||
// Work out Fp = psq_max/ psq...
|
||||
//////////////////////////////////
|
||||
std::vector<int> latt_size = grid->GlobalDimensions();
|
||||
std::vector<int> coor(grid->_ndimension,0);
|
||||
for(int mu=0;mu<Nd;mu++) {
|
||||
|
||||
Real TwoPiL = M_PI * 2.0/ latt_size[mu];
|
||||
LatticeCoordinate(pmu,mu);
|
||||
pmu = TwoPiL * pmu ;
|
||||
psq = psq + 4.0*sin(pmu*0.5)*sin(pmu*0.5);
|
||||
}
|
||||
|
||||
Complex psqMax(16.0);
|
||||
Fp = psqMax*one/psq;
|
||||
|
||||
static int once;
|
||||
if ( once == 0 ) {
|
||||
std::cout << " Fp " << Fp <<std::endl;
|
||||
once ++;
|
||||
}
|
||||
pokeSite(TComplex(1.0),Fp,coor);
|
||||
|
||||
dmuAmu_p = dmuAmu_p * Fp;
|
||||
|
||||
theFFT.FFT_all_dim(dmuAmu,dmuAmu_p,FFT::backward);
|
||||
|
||||
GaugeMat ciadmam(grid);
|
||||
Complex cialpha(0.0,-alpha);
|
||||
ciadmam = dmuAmu*cialpha;
|
||||
SU<Nc>::taExp(ciadmam,g);
|
||||
|
||||
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
|
||||
|
||||
SU<Nc>::GaugeTransform(U,g);
|
||||
|
||||
return trG;
|
||||
}
|
||||
|
||||
static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu) {
|
||||
GridBase *grid = g._grid;
|
||||
Complex cialpha(0.0,-alpha);
|
||||
GaugeMat ciadmam(grid);
|
||||
DmuAmu(A,dmuAmu);
|
||||
ciadmam = dmuAmu*cialpha;
|
||||
SU<Nc>::taExp(ciadmam,g);
|
||||
}
|
||||
/*
|
||||
////////////////////////////////////////////////////////////////
|
||||
// NB The FT for fields living on links has an extra phase in it
|
||||
// Could add these to the FFT class as a later task since this code
|
||||
// might be reused elsewhere ????
|
||||
////////////////////////////////////////////////////////////////
|
||||
static void InverseFourierTransformAmu(FFT &theFFT,const std::vector<GaugeMat> &Ap,std::vector<GaugeMat> &Ax) {
|
||||
GridBase * grid = theFFT.Grid();
|
||||
std::vector<int> latt_size = grid->GlobalDimensions();
|
||||
|
||||
ComplexField pmu(grid);
|
||||
ComplexField pha(grid);
|
||||
GaugeMat Apha(grid);
|
||||
|
||||
Complex ci(0.0,1.0);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
Real TwoPiL = M_PI * 2.0/ latt_size[mu];
|
||||
LatticeCoordinate(pmu,mu);
|
||||
pmu = TwoPiL * pmu ;
|
||||
pha = exp(pmu * (0.5 *ci)); // e(ipmu/2) since Amu(x+mu/2)
|
||||
|
||||
Apha = Ap[mu] * pha;
|
||||
|
||||
theFFT.FFT_all_dim(Apha,Ax[mu],FFT::backward);
|
||||
}
|
||||
}
|
||||
static void FourierTransformAmu(FFT & theFFT,const std::vector<GaugeMat> &Ax,std::vector<GaugeMat> &Ap) {
|
||||
GridBase * grid = theFFT.Grid();
|
||||
std::vector<int> latt_size = grid->GlobalDimensions();
|
||||
|
||||
ComplexField pmu(grid);
|
||||
ComplexField pha(grid);
|
||||
Complex ci(0.0,1.0);
|
||||
|
||||
// Sign convention for FFTW calls:
|
||||
// A(x)= Sum_p e^ipx A(p) / V
|
||||
// A(p)= Sum_p e^-ipx A(x)
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
Real TwoPiL = M_PI * 2.0/ latt_size[mu];
|
||||
LatticeCoordinate(pmu,mu);
|
||||
pmu = TwoPiL * pmu ;
|
||||
pha = exp(-pmu * (0.5 *ci)); // e(+ipmu/2) since Amu(x+mu/2)
|
||||
|
||||
theFFT.FFT_all_dim(Ax[mu],Ap[mu],FFT::backward);
|
||||
Ap[mu] = Ap[mu] * pha;
|
||||
}
|
||||
}
|
||||
*/
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
int threads = GridThread::GetThreads();
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout( { vComplex::Nsimd(),1,1,1});
|
||||
std::vector<int> 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);
|
||||
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); // naughty seeding
|
||||
GridParallelRNG pRNG(&GRID); pRNG.SeedFixedIntegers(seeds);
|
||||
|
||||
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;
|
||||
|
||||
LatticeGaugeField Umu(&GRID);
|
||||
LatticeGaugeField Uorg(&GRID);
|
||||
LatticeColourMatrix g(&GRID); // Gauge xform
|
||||
|
||||
|
||||
SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
|
||||
Uorg=Umu;
|
||||
|
||||
SU3::RandomGaugeTransform(pRNG,Umu,g); // Unit gauge
|
||||
Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<<plaq << std::endl;
|
||||
|
||||
|
||||
|
||||
Real alpha=0.1;
|
||||
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-10, 1.0e-10);
|
||||
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << std::endl;
|
||||
|
||||
Uorg = Uorg - Umu;
|
||||
std::cout << " Norm Difference "<< norm2(Uorg) << std::endl;
|
||||
|
||||
|
||||
// std::cout<< "*****************************************************************" <<std::endl;
|
||||
// std::cout<< "* Testing Fourier accelerated fixing *" <<std::endl;
|
||||
// std::cout<< "*****************************************************************" <<std::endl;
|
||||
|
||||
// std::cout<< "*****************************************************************" <<std::endl;
|
||||
// std::cout<< "* Testing non-unit configuration *" <<std::endl;
|
||||
// std::cout<< "*****************************************************************" <<std::endl;
|
||||
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -93,10 +93,10 @@ int main (int argc, char ** argv)
|
||||
C=C-Ctilde;
|
||||
std::cout << "diff scalar "<<norm2(C) << std::endl;
|
||||
|
||||
theFFT.FFT_dim(Stilde,S,0,FFT::forward); S=Stilde; std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,1,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,2,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,3,FFT::forward);std::cout << theFFT.MFlops()<<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,0,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<< " "<<theFFT.USec() <<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,1,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<< " "<<theFFT.USec() <<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,2,FFT::forward); S=Stilde;std::cout << theFFT.MFlops()<< " "<<theFFT.USec() <<std::endl;
|
||||
theFFT.FFT_dim(Stilde,S,3,FFT::forward);std::cout << theFFT.MFlops()<<" "<<theFFT.USec() <<std::endl;
|
||||
|
||||
SpinMatrixF Sp;
|
||||
Sp = zero; Sp = Sp+cVol;
|
||||
|
||||
138
tests/core/Test_poisson_fft.cc
Normal file
138
tests/core/Test_poisson_fft.cc
Normal file
@ -0,0 +1,138 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_poisson_fft.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;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
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;
|
||||
|
||||
int N=128;
|
||||
int N2=64;
|
||||
int W=16;
|
||||
int D=8;
|
||||
std::vector<int> latt_size ({N,N});
|
||||
std::vector<int> simd_layout({vComplexD::Nsimd(),1});
|
||||
std::vector<int> mpi_layout ({1,1});
|
||||
|
||||
int vol = 1;
|
||||
int nd = latt_size.size();
|
||||
for(int d=0;d<nd;d++){
|
||||
vol = vol * latt_size[d];
|
||||
}
|
||||
|
||||
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
|
||||
|
||||
LatticeComplexD pos(&GRID);
|
||||
LatticeComplexD zz(&GRID);
|
||||
LatticeComplexD neg(&GRID);
|
||||
LatticeInteger coor(&GRID);
|
||||
LatticeComplexD Charge(&GRID);
|
||||
LatticeComplexD ChargeTilde(&GRID);
|
||||
LatticeComplexD V(&GRID);
|
||||
LatticeComplexD Vtilde(&GRID);
|
||||
|
||||
pos = ComplexD(1.0,0.0);
|
||||
neg = -pos;
|
||||
zz = ComplexD(0.0,0.0);
|
||||
|
||||
Charge=zero;
|
||||
|
||||
// Parallel plate capacitor
|
||||
{
|
||||
int mu=0;
|
||||
LatticeCoordinate(coor,mu);
|
||||
|
||||
Charge=where(coor==Integer(N2-D),pos,zz);
|
||||
Charge=where(coor==Integer(N2+D),neg,Charge);
|
||||
}
|
||||
|
||||
{
|
||||
int mu=1;
|
||||
LatticeCoordinate(coor,mu);
|
||||
Charge=where(coor<Integer(N2-W),zz,Charge);
|
||||
Charge=where(coor>Integer(N2+W),zz,Charge);
|
||||
}
|
||||
|
||||
// std::cout << Charge <<std::endl;
|
||||
|
||||
std::vector<LatticeComplexD> k(4,&GRID);
|
||||
LatticeComplexD ksq(&GRID);
|
||||
|
||||
ksq=zero;
|
||||
for(int mu=0;mu<nd;mu++) {
|
||||
|
||||
Integer L=latt_size[mu];
|
||||
|
||||
LatticeCoordinate(coor,mu);
|
||||
LatticeCoordinate(k[mu],mu);
|
||||
|
||||
k[mu] = where ( coor > (L/2), k[mu]-L, k[mu]);
|
||||
|
||||
// std::cout << k[mu]<<std::endl;
|
||||
|
||||
RealD TwoPiL = M_PI * 2.0/ L;
|
||||
|
||||
k[mu] = TwoPiL * k[mu];
|
||||
|
||||
ksq = ksq + k[mu]*k[mu];
|
||||
|
||||
}
|
||||
|
||||
// D^2 V = - rho
|
||||
// ksq Vtilde = rhoTilde
|
||||
// Vtilde = rhoTilde/Ksq
|
||||
// Fix zero of potential : Vtilde(0) = 0;
|
||||
std::vector<int> zero_mode(nd,0);
|
||||
TComplexD Tone = ComplexD(1.0,0.0);
|
||||
pokeSite(Tone,ksq,zero_mode);
|
||||
|
||||
// std::cout << "Charge\n" << Charge <<std::endl;
|
||||
|
||||
FFT theFFT(&GRID);
|
||||
theFFT.FFT_all_dim(ChargeTilde,Charge,FFT::forward);
|
||||
// std::cout << "Rhotilde\n" << ChargeTilde <<std::endl;
|
||||
|
||||
Vtilde = ChargeTilde / ksq;
|
||||
// std::cout << "Vtilde\n" << Vtilde <<std::endl;
|
||||
|
||||
TComplexD Tzero = ComplexD(0.0,0.0);
|
||||
pokeSite(Tzero,Vtilde,zero_mode);
|
||||
|
||||
theFFT.FFT_all_dim(V,Vtilde,FFT::backward);
|
||||
|
||||
std::cout << "V\n" << V <<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
Reference in New Issue
Block a user