portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-10-13 12:44:42 +01:00
Code Releases Activity

Changes for Nd=3

Browse Source
This commit is contained in:
Peter Boyle
2025-10-03 12:17:13 -04:00
parent d418f78352
commit 35e10a1159
31 changed files with 256 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
Grid/GridQCDcore.h
View File

@@ -37,6 +37,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <Grid/qcd/QCD.h> #include <Grid/qcd/QCD.h>
#include <Grid/qcd/spin/Spin.h> #include <Grid/qcd/spin/Spin.h>
#include <Grid/qcd/gparity/Gparity.h> #include <Grid/qcd/gparity/Gparity.h>
#include <Grid/qcd/spin/Pauli.h> // depends on Gparity
#include <Grid/qcd/utils/Utils.h> #include <Grid/qcd/utils/Utils.h>
#include <Grid/qcd/representations/Representations.h> #include <Grid/qcd/representations/Representations.h>
NAMESPACE_CHECK(GridQCDCore); NAMESPACE_CHECK(GridQCDCore);

1
Grid/cartesian/Cartesian_base.h
View File

@@ -90,6 +90,7 @@ public:
// Checkerboarding interface is virtual and overridden by // Checkerboarding interface is virtual and overridden by
// GridCartesian / GridRedBlackCartesian // GridCartesian / GridRedBlackCartesian
//////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////
virtual int Icosahedral(void) { return 0;}
virtual int CheckerBoarded(int dim) =0; virtual int CheckerBoarded(int dim) =0;
virtual int CheckerBoard(const Coordinate &site)=0; virtual int CheckerBoard(const Coordinate &site)=0;
virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0; virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;

3
Grid/qcd/QCD.h
View File

@@ -49,7 +49,7 @@ static constexpr int Tm = 7;
static constexpr int Nc=Config_Nc; static constexpr int Nc=Config_Nc;
static constexpr int Ns=4; static constexpr int Ns=4;
static constexpr int Nd=4; static constexpr int Nd=Config_Nd;
static constexpr int Nhs=2; // half spinor static constexpr int Nhs=2; // half spinor
static constexpr int Nds=8; // double stored gauge field static constexpr int Nds=8; // double stored gauge field
static constexpr int Ngp=2; // gparity index range static constexpr int Ngp=2; // gparity index range
@@ -75,6 +75,7 @@ static constexpr int InverseYes=1;
//typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE; //typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE;
const int SpinorIndex = 2; const int SpinorIndex = 2;
const int PauliIndex = 2; //TensorLevel counts from the bottom!
template<typename T> struct isSpinor { template<typename T> struct isSpinor {
static constexpr bool value = (SpinorIndex==T::TensorLevel); static constexpr bool value = (SpinorIndex==T::TensorLevel);
}; };

4
Grid/qcd/action/fermion/DomainWallVec5dImpl.h
View File

@@ -123,10 +123,10 @@ public:
GaugeGrid->LocalIndexToLocalCoor(lidx, lcoor); GaugeGrid->LocalIndexToLocalCoor(lidx, lcoor);
peekLocalSite(ScalarUmu, Umu_v, lcoor); peekLocalSite(ScalarUmu, Umu_v, lcoor);
for (int mu = 0; mu < 4; mu++) ScalarUds(mu) = ScalarUmu(mu); for (int mu = 0; mu < Nd; mu++) ScalarUds(mu) = ScalarUmu(mu);
peekLocalSite(ScalarUmu, Uadj_v, lcoor); peekLocalSite(ScalarUmu, Uadj_v, lcoor);
for (int mu = 0; mu < 4; mu++) ScalarUds(mu + 4) = ScalarUmu(mu); for (int mu = 0; mu < Nd; mu++) ScalarUds(mu + Nd) = ScalarUmu(mu);
pokeLocalSite(ScalarUds, Uds_v, lcoor); pokeLocalSite(ScalarUds, Uds_v, lcoor);
}); });

9
Grid/qcd/action/fermion/Fermion.h
View File

@@ -85,6 +85,15 @@ NAMESPACE_CHECK(DomainWall);
#include <Grid/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h> #include <Grid/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h>
#include <Grid/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h> #include <Grid/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h>
NAMESPACE_CHECK(Overlap); NAMESPACE_CHECK(Overlap);
///////////////////////////////////////////////////////////////////////////////
// Two spin wilson fermion based
///////////////////////////////////////////////////////////////////////////////
#include <Grid/qcd/action/fermion/TwoSpinWilsonFermion3plus1D.h>
NAMESPACE_CHECK(TwoSpinWilson);
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////
// G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code // G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////

3
Grid/qcd/action/fermion/FermionCore.h
View File

@@ -41,8 +41,9 @@ NAMESPACE_CHECK(Compressor);
NAMESPACE_CHECK(FermionOperatorImpl); NAMESPACE_CHECK(FermionOperatorImpl);
#include <Grid/qcd/action/fermion/FermionOperator.h> #include <Grid/qcd/action/fermion/FermionOperator.h>
NAMESPACE_CHECK(FermionOperator); NAMESPACE_CHECK(FermionOperator);
#include <Grid/qcd/action/fermion/WilsonKernels.h> //used by all wilson type fermions #include <Grid/qcd/action/fermion/WilsonKernels.h> //used by all wilson type fermions
#include <Grid/qcd/action/fermion/StaggeredKernels.h> //used by all wilson type fermions #include <Grid/qcd/action/fermion/StaggeredKernels.h> //used by all wilson type fermions
#include <Grid/qcd/action/fermion/TwoSpinWilsonKernels.h> //used for 3D fermions, pauli in place of Dirac
NAMESPACE_CHECK(Kernels); NAMESPACE_CHECK(Kernels);
#endif #endif

6
Grid/qcd/action/fermion/FermionOperatorImpl.h
View File

@@ -180,6 +180,12 @@ NAMESPACE_CHECK(ImplGparityWilson);
#include <Grid/qcd/action/fermion/StaggeredImpl.h> #include <Grid/qcd/action/fermion/StaggeredImpl.h>
NAMESPACE_CHECK(ImplStaggered); NAMESPACE_CHECK(ImplStaggered);
/////////////////////////////////////////////////////////////////////////////
// Two component spinor Wilson action for 3d / Boston
/////////////////////////////////////////////////////////////////////////////
#include <Grid/qcd/action/fermion/TwoSpinWilsonImpl.h>
NAMESPACE_CHECK(ImplTwoSpinWilson);
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////
// Single flavour one component spinors with colour index. 5d vec // Single flavour one component spinors with colour index. 5d vec
///////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////

6
Grid/qcd/action/fermion/GparityWilsonImpl.h
View File

@@ -274,7 +274,7 @@ public:
autoView( Uds_v , Uds, CpuWrite); autoView( Uds_v , Uds, CpuWrite);
autoView( Utmp_v, Utmp, CpuWrite); autoView( Utmp_v, Utmp, CpuWrite);
thread_foreach(ss,Utmp_v,{ thread_foreach(ss,Utmp_v,{
Uds_v[ss](0)(mu+4) = Utmp_v[ss](); Uds_v[ss](0)(mu+Nd) = Utmp_v[ss]();
}); });
} }
Utmp = Uconj; Utmp = Uconj;
@@ -286,7 +286,7 @@ public:
autoView( Uds_v , Uds, CpuWrite); autoView( Uds_v , Uds, CpuWrite);
autoView( Utmp_v, Utmp, CpuWrite); autoView( Utmp_v, Utmp, CpuWrite);
thread_foreach(ss,Utmp_v,{ thread_foreach(ss,Utmp_v,{
Uds_v[ss](1)(mu+4) = Utmp_v[ss](); Uds_v[ss](1)(mu+Nd) = Utmp_v[ss]();
}); });
} }
} }
@@ -320,7 +320,7 @@ public:
} }
Uconj = conjugate(*Upoke); Uconj = conjugate(*Upoke);
pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4); pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + Nd);
} }
} }

2
Grid/qcd/action/fermion/ImprovedStaggeredFermion.h
View File

@@ -36,6 +36,8 @@ public:
static const std::vector<int> directions; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static const int npoint = 16; static const int npoint = 16;
static std::vector<int> MakeDirections(void);
static std::vector<int> MakeDisplacements(void);
}; };
template <class Impl> template <class Impl>

2
Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h
View File

@@ -40,6 +40,8 @@ public:
static const std::vector<int> directions; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
const int npoint = 16; const int npoint = 16;
static std::vector<int> MakeDirections(void);
static std::vector<int> MakeDisplacements(void);
}; };
template<class Impl> template<class Impl>

2
Grid/qcd/action/fermion/NaiveStaggeredFermion.h
View File

@@ -36,6 +36,8 @@ public:
static const std::vector<int> directions; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static const int npoint = 8; static const int npoint = 8;
static std::vector<int> MakeDirections(void);
static std::vector<int> MakeDisplacements(void);
}; };
template <class Impl> template <class Impl>

6
Grid/qcd/action/fermion/StaggeredImpl.h
View File

@@ -141,9 +141,9 @@ public:
Udag = Udag *phases; Udag = Udag *phases;
InsertGaugeField(Uds,U,mu); InsertGaugeField(Uds,U,mu);
InsertGaugeField(Uds,Udag,mu+4); InsertGaugeField(Uds,Udag,mu+Nd);
// PokeIndex<LorentzIndex>(Uds, U, mu); // PokeIndex<LorentzIndex>(Uds, U, mu);
// PokeIndex<LorentzIndex>(Uds, Udag, mu + 4); // PokeIndex<LorentzIndex>(Uds, Udag, mu + Nd);
// 3 hop based on thin links. Crazy huh ? // 3 hop based on thin links. Crazy huh ?
U = PeekIndex<LorentzIndex>(Uthin, mu); U = PeekIndex<LorentzIndex>(Uthin, mu);
@@ -156,7 +156,7 @@ public:
UUUdag = UUUdag *phases; UUUdag = UUUdag *phases;
InsertGaugeField(UUUds,UUU,mu); InsertGaugeField(UUUds,UUU,mu);
InsertGaugeField(UUUds,UUUdag,mu+4); InsertGaugeField(UUUds,UUUdag,mu+Nd);
} }
} }

2
Grid/qcd/action/fermion/WilsonFermion.h
View File

@@ -38,6 +38,8 @@ public:
static int MortonOrder; static int MortonOrder;
static const std::vector<int> directions; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static std::vector<int> MakeDirections(void);
static std::vector<int> MakeDisplacements(void);
static const int npoint = 8; static const int npoint = 8;
}; };

2
Grid/qcd/action/fermion/WilsonFermion5D.h
View File

@@ -62,6 +62,8 @@ public:
static const std::vector<int> directions; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static constexpr int npoint = 8; static constexpr int npoint = 8;
static std::vector<int> MakeDirections(void);
static std::vector<int> MakeDisplacements(void);
}; };
template<class Impl> template<class Impl>

2
Grid/qcd/action/fermion/WilsonImpl.h
View File

@@ -166,7 +166,7 @@ public:
U = adj(Cshift(U, mu, -1)); U = adj(Cshift(U, mu, -1));
U = where(coor == 0, conjugate(phase) * U, U); U = where(coor == 0, conjugate(phase) * U, U);
PokeIndex<LorentzIndex>(Uds, U, mu + 4); PokeIndex<LorentzIndex>(Uds, U, mu + Nd);
} }
} }

12
Grid/qcd/action/fermion/WilsonTMFermion5D.h
View File

@@ -56,7 +56,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
Frbgrid, Frbgrid,
Ugrid, Ugrid,
Urbgrid, Urbgrid,
4.0,p) Nd*1.0,p)
{ {
update(_mass,_mu); update(_mass,_mu);
@@ -83,7 +83,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
//axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in
for (int s=0;s<(int)this->mass.size();s++) { for (int s=0;s<(int)this->mass.size();s++) {
ComplexD a = 4.0+this->mass[s]; ComplexD a = Nd*1.0+this->mass[s];
ComplexD b(0.0,this->mu[s]); ComplexD b(0.0,this->mu[s]);
axpbg5y_ssp(out,a,in,b,in,s,s); axpbg5y_ssp(out,a,in,b,in,s,s);
} }
@@ -92,7 +92,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
virtual void MooeeDag(const FermionField &in, FermionField &out) { virtual void MooeeDag(const FermionField &in, FermionField &out) {
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
for (int s=0;s<(int)this->mass.size();s++) { for (int s=0;s<(int)this->mass.size();s++) {
ComplexD a = 4.0+this->mass[s]; ComplexD a = Nd*1.0+this->mass[s];
ComplexD b(0.0,-this->mu[s]); ComplexD b(0.0,-this->mu[s]);
axpbg5y_ssp(out,a,in,b,in,s,s); axpbg5y_ssp(out,a,in,b,in,s,s);
} }
@@ -101,7 +101,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
for (int s=0;s<(int)this->mass.size();s++) { for (int s=0;s<(int)this->mass.size();s++) {
RealD m = this->mass[s]; RealD m = this->mass[s];
RealD tm = this->mu[s]; RealD tm = this->mu[s];
RealD mtil = 4.0+this->mass[s]; RealD mtil = Nd*1.0+this->mass[s];
RealD sq = mtil*mtil+tm*tm; RealD sq = mtil*mtil+tm*tm;
ComplexD a = mtil/sq; ComplexD a = mtil/sq;
ComplexD b(0.0, -tm /sq); ComplexD b(0.0, -tm /sq);
@@ -112,7 +112,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
for (int s=0;s<(int)this->mass.size();s++) { for (int s=0;s<(int)this->mass.size();s++) {
RealD m = this->mass[s]; RealD m = this->mass[s];
RealD tm = this->mu[s]; RealD tm = this->mu[s];
RealD mtil = 4.0+this->mass[s]; RealD mtil = Nd*1.0+this->mass[s];
RealD sq = mtil*mtil+tm*tm; RealD sq = mtil*mtil+tm*tm;
ComplexD a = mtil/sq; ComplexD a = mtil/sq;
ComplexD b(0.0,tm /sq); ComplexD b(0.0,tm /sq);
@@ -126,7 +126,7 @@ class WilsonTMFermion5D : public WilsonFermion5D<Impl>
this->Dhop(in, out, DaggerNo); this->Dhop(in, out, DaggerNo);
FermionField tmp(out.Grid()); FermionField tmp(out.Grid());
for (int s=0;s<(int)this->mass.size();s++) { for (int s=0;s<(int)this->mass.size();s++) {
ComplexD a = 4.0+this->mass[s]; ComplexD a = Nd*1.0+this->mass[s];
ComplexD b(0.0,this->mu[s]); ComplexD b(0.0,this->mu[s]);
axpbg5y_ssp(tmp,a,in,b,in,s,s); axpbg5y_ssp(tmp,a,in,b,in,s,s);
} }

2
Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h
View File

@@ -240,7 +240,7 @@ void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::ve
this->ceo.resize(Ls); this->ceo.resize(Ls);
for(int i=0; i<Ls; ++i){ for(int i=0; i<Ls; ++i){
this->bee[i] = 4.0 - this->M5 + 1.0; this->bee[i] = Nd*1.0 - this->M5 + 1.0;
this->cee[i] = 1.0; this->cee[i] = 1.0;
} }

14
Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h
View File

@@ -61,7 +61,7 @@ WilsonCloverFermion<Impl, CloverHelpers>::WilsonCloverFermion(GaugeField&
diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
} else { } else {
csw_r = _csw_r * 0.5; csw_r = _csw_r * 0.5;
diag_mass = 4.0 + _mass; diag_mass = Nd*1.0 + _mass;
} }
csw_t = _csw_t * 0.5; csw_t = _csw_t * 0.5;
@@ -297,9 +297,9 @@ void WilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField &force, const F
{ {
if (mu == nu) if (mu == nu)
continue; continue;
RealD factor; RealD factor;
if (nu == 4 || mu == 4) if (nu == (Nd-1) || mu == (Nd-1)) // This was a bug - surely mu/nu is NEVER 4 but rather (Nd-1)=3 ??
{ {
factor = 2.0 * csw_t; factor = 2.0 * csw_t;
} }
@@ -307,9 +307,11 @@ void WilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField &force, const F
{ {
factor = 2.0 * csw_r; factor = 2.0 * csw_r;
} }
PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked if ( mu < Nd && nu < Nd ) { // Allow to restrict range to Nd=3, but preserve orders of SigmaMuNu in table by counting ALL
Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok
force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked
}
count++; count++;
} }

36
Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h
View File

@@ -63,10 +63,10 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
Dirichlet(0) Dirichlet(0)
{ {
// some assertions // some assertions
assert(FiveDimGrid._ndimension==5); assert(FiveDimGrid._ndimension==Nd+1);
assert(FourDimGrid._ndimension==4); assert(FourDimGrid._ndimension==Nd);
assert(FourDimRedBlackGrid._ndimension==4); assert(FourDimRedBlackGrid._ndimension==Nd);
assert(FiveDimRedBlackGrid._ndimension==5); assert(FiveDimRedBlackGrid._ndimension==Nd+1);
assert(FiveDimRedBlackGrid._checker_dim==1); // Don't checker the s direction assert(FiveDimRedBlackGrid._checker_dim==1); // Don't checker the s direction
// extent of fifth dim and not spread out // extent of fifth dim and not spread out
@@ -76,7 +76,7 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
assert(FiveDimRedBlackGrid._processors[0] ==1); assert(FiveDimRedBlackGrid._processors[0] ==1);
// Other dimensions must match the decomposition of the four-D fields // Other dimensions must match the decomposition of the four-D fields
for(int d=0;d<4;d++){ for(int d=0;d<Nd;d++){
assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]); assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]);
assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]); assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]);
@@ -93,11 +93,13 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
if ( p.dirichlet.size() == Nd+1) { if ( p.dirichlet.size() == Nd+1) {
Coordinate block = p.dirichlet; Coordinate block = p.dirichlet;
if ( block[0] || block[1] || block[2] || block[3] || block[4] ){ for(int d=0;d<Nd+1;d++) {
Dirichlet = 1; if ( block[d] ){
std::cout << GridLogMessage << " WilsonFermion: non-trivial Dirichlet condition "<< block << std::endl; Dirichlet = 1;
std::cout << GridLogMessage << " WilsonFermion: partial Dirichlet "<< p.partialDirichlet << std::endl; std::cout << GridLogMessage << " WilsonFermion: non-trivial Dirichlet condition "<< block << std::endl;
Block = block; std::cout << GridLogMessage << " WilsonFermion: partial Dirichlet "<< p.partialDirichlet << std::endl;
Block = block;
}
} }
} else { } else {
Coordinate block(Nd+1,0); Coordinate block(Nd+1,0);
@@ -112,7 +114,7 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
assert(FiveDimGrid._simd_layout[0] ==nsimd); assert(FiveDimGrid._simd_layout[0] ==nsimd);
assert(FiveDimRedBlackGrid._simd_layout[0]==nsimd); assert(FiveDimRedBlackGrid._simd_layout[0]==nsimd);
for(int d=0;d<4;d++){ for(int d=0;d<Nd;d++){
assert(FourDimGrid._simd_layout[d]==1); assert(FourDimGrid._simd_layout[d]==1);
assert(FourDimRedBlackGrid._simd_layout[d]==1); assert(FourDimRedBlackGrid._simd_layout[d]==1);
assert(FiveDimRedBlackGrid._simd_layout[d+1]==1); assert(FiveDimRedBlackGrid._simd_layout[d+1]==1);
@@ -183,8 +185,8 @@ void WilsonFermion5D<Impl>::DhopDir(const FermionField &in, FermionField &out,in
// assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t; // assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t;
int skip = (disp==1) ? 0 : 1; int skip = (disp==1) ? 0 : 1;
int dirdisp = dir+skip*4; int dirdisp = dir+skip*Nd;
int gamma = dir+(1-skip)*4; int gamma = dir+(1-skip)*Nd;
Compressor compressor(DaggerNo); Compressor compressor(DaggerNo);
Stencil.HaloExchange(in,compressor); Stencil.HaloExchange(in,compressor);
@@ -483,7 +485,7 @@ void WilsonFermion5D<Impl>::DW(const FermionField &in, FermionField &out,int dag
{ {
out.Checkerboard()=in.Checkerboard(); out.Checkerboard()=in.Checkerboard();
Dhop(in,out,dag); // -0.5 is included Dhop(in,out,dag); // -0.5 is included
axpy(out,4.0-M5,in,out); axpy(out,Nd*1.0-M5,in,out);
} }
template <class Impl> template <class Impl>
void WilsonFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out) void WilsonFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out)
@@ -509,7 +511,7 @@ template <class Impl>
void WilsonFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out) void WilsonFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out)
{ {
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
typename FermionField::scalar_type scal(4.0 + M5); typename FermionField::scalar_type scal(Nd*1.0 + M5);
out = scal * in; out = scal * in;
} }
@@ -524,7 +526,7 @@ template<class Impl>
void WilsonFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out) void WilsonFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out)
{ {
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
out = (1.0/(4.0 + M5))*in; out = (1.0/(Nd*1.0 + M5))*in;
} }
template<class Impl> template<class Impl>
@@ -635,7 +637,7 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt_5d(FermionField &out,const
A = one / (abs(W) * sinha * 2.0) * one / (sinhaLs * 2.0); A = one / (abs(W) * sinha * 2.0) * one / (sinhaLs * 2.0);
F = eaLs * (one - Wea + (Wema - one) * mass*mass); F = eaLs * (one - Wea + (Wema - one) * mass*mass);
F = F + emaLs * (Wema - one + (one - Wea) * mass*mass); F = F + emaLs * (Wema - one + (one - Wea) * mass*mass);
F = F - abs(W) * sinha * 4.0 * mass; F = F - abs(W) * sinha * (Nd* 1.0) * mass;
Bpp = (A/F) * (ema2Ls - one) * (one - Wema) * (one - mass*mass * one); Bpp = (A/F) * (ema2Ls - one) * (one - Wema) * (one - mass*mass * one);
Bmm = (A/F) * (one - ea2Ls) * (one - Wea) * (one - mass*mass * one); Bmm = (A/F) * (one - ea2Ls) * (one - Wea) * (one - mass*mass * one);

10
Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h
View File

@@ -63,7 +63,7 @@ WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
if (anisotropyCoeff.isAnisotropic){ if (anisotropyCoeff.isAnisotropic){
diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0); diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0);
} else { } else {
diag_mass = 4.0 + mass; diag_mass = Nd*1.0 + mass;
} }
int vol4; int vol4;
@@ -354,8 +354,8 @@ void WilsonFermion<Impl>::DhopDir(const FermionField &in, FermionField &out, int
Stencil.HaloExchange(in, compressor); Stencil.HaloExchange(in, compressor);
int skip = (disp == 1) ? 0 : 1; int skip = (disp == 1) ? 0 : 1;
int dirdisp = dir + skip * 4; int dirdisp = dir + skip * Nd;
int gamma = dir + (1 - skip) * 4; int gamma = dir + (1 - skip) * Nd;
DhopDirCalc(in, out, dirdisp, gamma, DaggerNo); DhopDirCalc(in, out, dirdisp, gamma, DaggerNo);
}; };
@@ -370,8 +370,8 @@ void WilsonFermion<Impl>::DhopDirAll(const FermionField &in, std::vector<Fermion
for(int disp=-1;disp<=1;disp+=2){ for(int disp=-1;disp<=1;disp+=2){
int skip = (disp == 1) ? 0 : 1; int skip = (disp == 1) ? 0 : 1;
int dirdisp = dir + skip * 4; int dirdisp = dir + skip * Nd;
int gamma = dir + (1 - skip) * 4; int gamma = dir + (1 - skip) * Nd;
DhopDirCalc(in, out[dirdisp], dirdisp, gamma, DaggerNo); DhopDirCalc(in, out[dirdisp], dirdisp, gamma, DaggerNo);
} }

2
Grid/qcd/action/fermion/implementation/WilsonKernelsHandGparityImplementation.h
View File

@@ -97,7 +97,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
distance = st._distances[DIR]; \ distance = st._distances[DIR]; \
sl = st._simd_layout[direction]; \ sl = st._simd_layout[direction]; \
inplace_twist = 0; \ inplace_twist = 0; \
if(SE->_around_the_world && st.parameters.twists[DIR % 4]){ \ if(SE->_around_the_world && st.parameters.twists[DIR % Nd]){ \
if(sl == 1){ \ if(sl == 1){ \
g = (F+1) % 2; \ g = (F+1) % 2; \
}else{ \ }else{ \

26
Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermion5DInstantiation.cc
View File

@@ -32,8 +32,30 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
// S-direction is INNERMOST and takes no part in the parity. // S-direction is INNERMOST and takes no part in the parity.
const std::vector<int> ImprovedStaggeredFermion5DStatic::directions({1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4}); const std::vector<int> ImprovedStaggeredFermion5DStatic::directions(ImprovedStaggeredFermion5DStatic::MakeDirections());
const std::vector<int> ImprovedStaggeredFermion5DStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); const std::vector<int> ImprovedStaggeredFermion5DStatic::displacements(ImprovedStaggeredFermion5DStatic::MakeDisplacements());
std::vector<int> ImprovedStaggeredFermion5DStatic::MakeDirections(void)
{
std::vector<int> directions(4*Nd);
for(int d=0;d<Nd;d++){
directions[d+Nd*0] = d+1;
directions[d+Nd*1] = d+1;
directions[d+Nd*2] = d+1;
directions[d+Nd*3] = d+1;
}
return directions;
}
std::vector<int> ImprovedStaggeredFermion5DStatic::MakeDisplacements(void)
{
std::vector<int> displacements(4*Nd);
for(int d=0;d<Nd;d++){
displacements[d+Nd*0] =+1;
displacements[d+Nd*1] =-1;
displacements[d+Nd*2] =+3;
displacements[d+Nd*3] =-3;
}
return displacements;
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

23
Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermionInstantiation.cc
View File

@@ -32,5 +32,26 @@ NAMESPACE_BEGIN(Grid);
const std::vector<int> ImprovedStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3}); const std::vector<int> ImprovedStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3});
const std::vector<int> ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); const std::vector<int> ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3});
std::vector<int> ImprovedStaggeredFermionStatic::MakeDirections(void)
{
std::vector<int> directions(4*Nd);
for(int d=0;d<Nd;d++){
directions[d+Nd*0] = d;
directions[d+Nd*1] = d;
directions[d+Nd*2] = d;
directions[d+Nd*3] = d;
}
return directions;
}
std::vector<int> ImprovedStaggeredFermionStatic::MakeDisplacements(void)
{
std::vector<int> displacements(4*Nd);
for(int d=0;d<Nd;d++){
displacements[d+Nd*0] =+1;
displacements[d+Nd*1] =-1;
displacements[d+Nd*2] =+3;
displacements[d+Nd*3] =-3;
}
return displacements;
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

24
Grid/qcd/action/fermion/instantiation/NaiveStaggeredFermionInstantiation.cc
View File

@@ -30,7 +30,27 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
const std::vector<int> NaiveStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); //const std::vector<int> NaiveStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3});
const std::vector<int> NaiveStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); //const std::vector<int> NaiveStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1});
const std::vector<int> NaiveStaggeredFermionStatic::directions(NaiveStaggeredFermionStatic::MakeDirections());
const std::vector<int> NaiveStaggeredFermionStatic::displacements(NaiveStaggeredFermionStatic::MakeDisplacements());
std::vector<int> NaiveStaggeredFermionStatic::MakeDirections(void)
{
std::vector<int> directions(4*Nd);
for(int d=0;d<Nd;d++){
directions[d+Nd*0] = d;
directions[d+Nd*1] = d;
}
return directions;
}
std::vector<int> NaiveStaggeredFermionStatic::MakeDisplacements(void)
{
std::vector<int> displacements(4*Nd);
for(int d=0;d<Nd;d++){
displacements[d+Nd*0] =+1;
displacements[d+Nd*1] =-1;
}
return displacements;
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

24
Grid/qcd/action/fermion/instantiation/WilsonFermion5DInstantiation.cc
View File

@@ -34,8 +34,28 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
// S-direction is INNERMOST and takes no part in the parity. // S-direction is INNERMOST and takes no part in the parity.
const std::vector<int> WilsonFermion5DStatic::directions ({1,2,3,4, 1, 2, 3, 4});
const std::vector<int> WilsonFermion5DStatic::displacements({1,1,1,1,-1,-1,-1,-1}); const std::vector<int> WilsonFermion5DStatic::directions (WilsonFermion5DStatic::MakeDirections());
const std::vector<int> WilsonFermion5DStatic::displacements(WilsonFermion5DStatic::MakeDisplacements());
std::vector<int> WilsonFermion5DStatic::MakeDirections (void)
{
std::vector<int> directions(2*Nd);
for(int d=0;d<Nd;d++){
directions[d] = d+1;
directions[d+Nd] = d+1;
}
return directions;
}
std::vector<int> WilsonFermion5DStatic::MakeDisplacements(void)
{
std::vector<int> displacements(2*Nd);
for(int d=0;d<Nd;d++){
displacements[d] = +1;
displacements[d+Nd] = -1;
}
return displacements;
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

22
Grid/qcd/action/fermion/instantiation/WilsonFermionInstantiation.cc
View File

@@ -33,9 +33,27 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
const std::vector<int> WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); const std::vector<int> WilsonFermionStatic::directions(WilsonFermionStatic::MakeDirections());
const std::vector<int> WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); const std::vector<int> WilsonFermionStatic::displacements(WilsonFermionStatic::MakeDisplacements());
int WilsonFermionStatic::HandOptDslash; int WilsonFermionStatic::HandOptDslash;
std::vector<int> WilsonFermionStatic::MakeDirections (void)
{
std::vector<int> directions(2*Nd);
for(int d=0;d<Nd;d++){
directions[d] = d;
directions[d+Nd] = d;
}
return directions;
}
std::vector<int> WilsonFermionStatic::MakeDisplacements(void)
{
std::vector<int> displacements(2*Nd);
for(int d=0;d<Nd;d++){
displacements[d] = +1;
displacements[d+Nd] = -1;
}
return displacements;
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

18
Grid/qcd/action/fermion/instantiation/generate_instantiations.sh
View File

@@ -36,11 +36,16 @@ DWF_IMPL_LIST=" \
ZWilsonImplF \ ZWilsonImplF \
ZWilsonImplD2 " ZWilsonImplD2 "
TWOSPIN_WILSON_IMPL_LIST=" \
TwoSpinWilsonImplF \
TwoSpinWilsonImplD "
GDWF_IMPL_LIST=" \ GDWF_IMPL_LIST=" \
GparityWilsonImplF \ GparityWilsonImplF \
GparityWilsonImplD " GparityWilsonImplD "
IMPL_LIST="$STAG_IMPL_LIST $WILSON_IMPL_LIST $DWF_IMPL_LIST $GDWF_IMPL_LIST" IMPL_LIST="$STAG_IMPL_LIST $WILSON_IMPL_LIST $DWF_IMPL_LIST $GDWF_IMPL_LIST $TWOSPIN_WILSON_IMPL_LIST"
for impl in $IMPL_LIST for impl in $IMPL_LIST
do do
@@ -110,7 +115,12 @@ do
done done
done done
CC_LIST=" \ CC_LIST="TwoSpinWilsonFermion3plus1DInstantiation.cc.master TwoSpinWilsonKernelsInstantiation.cc.master"
ImprovedStaggeredFermion5DInstantiation \
StaggeredKernelsInstantiation "
for impl in $TWOSPIN_WILSON_IMPL_LIST
do
for f in $CC_LIST
do
ln -f -s ../$f.cc.master $impl/$f$impl.cc
done
done

4
Grid/qcd/smearing/WilsonFlow.h
View File

@@ -158,8 +158,8 @@ RealD WilsonFlowBase<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeF
LatticeComplexD R(U.Grid()); LatticeComplexD R(U.Grid());
R = Zero(); R = Zero();
for(int mu=0;mu<3;mu++){ for(int mu=0;mu<Nd-1;mu++){
for(int nu=mu+1;nu<4;nu++){ for(int nu=mu+1;nu<Nd;nu++){
WilsonLoops<Gimpl>::FieldStrength(F, U, mu, nu); WilsonLoops<Gimpl>::FieldStrength(F, U, mu, nu);
R = R + trace(F*F); R = R + trace(F*F);
} }

45
Grid/qcd/utils/WilsonLoops.h
View File

@@ -179,20 +179,17 @@ public:
////////////////////////////////////////////////// //////////////////////////////////////////////////
// average over all x,y,z the temporal loop // average over all x,y,z the temporal loop
////////////////////////////////////////////////// //////////////////////////////////////////////////
static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4 static ComplexD avgPolyakovLoop(const GaugeField &Umu) {
GaugeMat Ut(Umu.Grid()), P(Umu.Grid()); GaugeMat Ut(Umu.Grid()), P(Umu.Grid());
ComplexD out; ComplexD out;
int T = Umu.Grid()->GlobalDimensions()[3]; uint64_t vol = Umu.Grid()->gSites();
int X = Umu.Grid()->GlobalDimensions()[0]; int T = Umu.Grid()->GlobalDimensions()[Nd-1];
int Y = Umu.Grid()->GlobalDimensions()[1]; Ut = peekLorentz(Umu,Nd-1); //Select temporal direction
int Z = Umu.Grid()->GlobalDimensions()[2];
Ut = peekLorentz(Umu,3); //Select temporal direction
P = Ut; P = Ut;
for (int t=1;t<T;t++){ for (int t=1;t<T;t++){
P = Gimpl::CovShiftForward(Ut,3,P); P = Gimpl::CovShiftForward(Ut,Nd-1,P);
} }
RealD norm = 1.0/(Nc*X*Y*Z*T); RealD norm = 1.0/(Nc*vol);
out = sum(trace(P))*norm; out = sum(trace(P))*norm;
return out; return out;
} }
@@ -215,7 +212,7 @@ public:
double vol = Umu.Grid()->gSites(); double vol = Umu.Grid()->gSites();
return p.real() / vol / (4.0 * Nc ) ; return p.real() / vol / (Nd * Nc ) ;
}; };
////////////////////////////////////////////////// //////////////////////////////////////////////////
@@ -740,6 +737,7 @@ public:
//cf https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 6 //cf https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 6
//output is the charge by timeslice: sum over timeslices to obtain the total //output is the charge by timeslice: sum over timeslices to obtain the total
static std::vector<Real> TimesliceTopologicalChargeMxN(const GaugeLorentz &U, int M, int N){ static std::vector<Real> TimesliceTopologicalChargeMxN(const GaugeLorentz &U, int M, int N){
// Audit: 4D epsilon is hard coded
assert(Nd == 4); assert(Nd == 4);
std::vector<std::vector<GaugeMat*> > F(Nd,std::vector<GaugeMat*>(Nd,nullptr)); std::vector<std::vector<GaugeMat*> > F(Nd,std::vector<GaugeMat*>(Nd,nullptr));
//Note F_numu = - F_munu //Note F_numu = - F_munu
@@ -829,6 +827,25 @@ public:
return out; return out;
} }
//Compute the 5Li topological charge density
static std::vector<Real> TopologicalChargeDensity5Li(const GaugeLorentz &U){
static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} };
std::vector<std::vector<Real> > loops = TimesliceTopologicalCharge5LiContributions(U);
double c5=1./20.;
double c4=1./5.-2.*c5;
double c3=(-64.+640.*c5)/45.;
double c2=(1-64.*c5)/9.;
double c1=(19.-55.*c5)/9.;
int Lt = loops[0].size();
std::vector<Real> out(Lt,0.);
for(int t=0;t<Lt;t++)
out[t] += c1*loops[0][t] + c2*loops[1][t] + c3*loops[2][t] + c4*loops[3][t] + c5*loops[4][t];
return out;
}
static Real TopologicalCharge5Li(const GaugeLorentz &U){ static Real TopologicalCharge5Li(const GaugeLorentz &U){
std::vector<Real> Qt = TimesliceTopologicalCharge5Li(U); std::vector<Real> Qt = TimesliceTopologicalCharge5Li(U);
Real Q = 0.; Real Q = 0.;
@@ -1455,7 +1472,7 @@ public:
////////////////////////////////////////////////// //////////////////////////////////////////////////
static Real sumWilsonLoop(const GaugeLorentz &Umu, static Real sumWilsonLoop(const GaugeLorentz &Umu,
const int R1, const int R2) { const int R1, const int R2) {
std::vector<GaugeMat> U(4, Umu.Grid()); std::vector<GaugeMat> U(Nd, Umu.Grid());
for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu); U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
@@ -1474,7 +1491,7 @@ public:
////////////////////////////////////////////////// //////////////////////////////////////////////////
static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu, static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu,
const int R1, const int R2) { const int R1, const int R2) {
std::vector<GaugeMat> U(4, Umu.Grid()); std::vector<GaugeMat> U(Nd, Umu.Grid());
for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu); U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
@@ -1492,8 +1509,8 @@ public:
// sum over all x,y,z,t and over all planes of spatial Wilson loop // sum over all x,y,z,t and over all planes of spatial Wilson loop
////////////////////////////////////////////////// //////////////////////////////////////////////////
static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu, static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu,
const int R1, const int R2) { const int R1, const int R2) {
std::vector<GaugeMat> U(4, Umu.Grid()); std::vector<GaugeMat> U(Nd, Umu.Grid());
for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu); U[mu] = PeekIndex<LorentzIndex>(Umu, mu);

4
Grid/util/Init.cc
View File

@@ -187,8 +187,8 @@ void GridParseLayout(char **argv,int argc,
Coordinate &latt_c, Coordinate &latt_c,
Coordinate &mpi_c) Coordinate &mpi_c)
{ {
auto mpi =std::vector<int>({1,1,1,1}); auto mpi =std::vector<int>(1,Nd);
auto latt=std::vector<int>({8,8,8,8}); auto latt=std::vector<int>(8,Nd);
GridThread::SetMaxThreads(); GridThread::SetMaxThreads();

17
configure.ac
View File

@@ -198,6 +198,8 @@ AC_ARG_ENABLE([Nc],
[ac_Nc=${enable_Nc}], [ac_Nc=3]) [ac_Nc=${enable_Nc}], [ac_Nc=3])
case ${ac_Nc} in case ${ac_Nc} in
1)
AC_DEFINE([Config_Nc],[1],[Gauge group Nc]);;
2) 2)
AC_DEFINE([Config_Nc],[2],[Gauge group Nc]);; AC_DEFINE([Config_Nc],[2],[Gauge group Nc]);;
3) 3)
@@ -211,6 +213,21 @@ case ${ac_Nc} in
*) *)
AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);; AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);;
esac esac
############### Nd
AC_ARG_ENABLE([Nd],
[AS_HELP_STRING([--enable-Nd=2|3|4],[enable default LGT dimension])],
[ac_Nd=${enable_Nd}], [ac_Nd=4])
case ${ac_Nd} in
2)
AC_DEFINE([Config_Nd],[2],[Gauge field dimension Nd]);;
3)
AC_DEFINE([Config_Nd],[3],[Gauge field dimension Nd]);;
4)
AC_DEFINE([Config_Nd],[4],[Gauge field dimension Nd]);;
*)
AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);;
esac
############### Symplectic group ############### Symplectic group
AC_ARG_ENABLE([Sp], AC_ARG_ENABLE([Sp],