portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

For G-parity BCs the Nd-1 direction is now assumed to be the time direction and setting a twist in this direction will apply antiperiodic BCs

Browse Source 
Added option to run Test_gparity with antiperiodic time BCs
This commit is contained in:
Christopher Kelly 2020-12-17 14:09:00 -05:00
parent 9e7bacb5a4
commit 249b6e61ec
2 changed files with 92 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -30,6 +30,19 @@ directory
NAMESPACE_BEGIN(Grid);
/*
Policy implementation for G-parity boundary conditions
Rather than treating the gauge field as a flavored field, the Grid implementation of G-parity treats the gauge field as a regular
field with complex conjugate boundary conditions. In order to ensure the second flavor interacts with the conjugate links and the first
with the regular links we overload the functionality of doubleStore, whose purpose is to store the gauge field and the barrel-shifted gauge field
to avoid communicating links when applying the Dirac operator, such that the double-stored field contains also a flavor index which maps to
either the link or the conjugate link. This flavored field is then used by multLink to apply the correct link to a spinor.
Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction.
mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
*/
template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
public:
@ -113,7 +126,7 @@ public:
|| ((distance== 1)&&(icoor[direction]==1))
|| ((distance==-1)&&(icoor[direction]==0));
permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu]; //only if we are going around the world
permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu] && mmu < Nd-1; //only if we are going around the world in a spatial direction
//Apply the links
int f_upper = permute_lane ? 1 : 0;
@ -139,10 +152,10 @@ public:
assert((distance == 1) || (distance == -1)); // nearest neighbour stencil hard code
assert((sl == 1) || (sl == 2));
if ( SE->_around_the_world && St.parameters.twists[mmu] ) {
//If this site is an global boundary site, perform the G-parity flavor twist
if ( mmu < Nd-1 && SE->_around_the_world && St.parameters.twists[mmu] ) {
if ( sl == 2 ) {
//Only do the twist for lanes on the edge of the physical node
ExtractBuffer<sobj> vals(Nsimd);
extract(chi,vals);
@ -197,6 +210,19 @@ public:
reg = memory;
}
//Poke 'poke_f0' onto flavor 0 and 'poke_f1' onto flavor 1 in direction mu of the doubled gauge field Uds
inline void pokeGparityDoubledGaugeField(DoubledGaugeField &Uds, const GaugeLinkField &poke_f0, const GaugeLinkField &poke_f1, const int mu){
autoView(poke_f0_v, poke_f0, CpuRead);
autoView(poke_f1_v, poke_f1, CpuRead);
autoView(Uds_v, Uds, CpuWrite);
thread_foreach(ss,poke_f0_v,{
Uds_v[ss](0)(mu) = poke_f0_v[ss]();
Uds_v[ss](1)(mu) = poke_f1_v[ss]();
});
}
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
{
conformable(Uds.Grid(),GaugeGrid);
@ -207,14 +233,16 @@ public:
GaugeLinkField Uconj(GaugeGrid);
Lattice<iScalar<vInteger> > coor(GaugeGrid);
for(int mu=0;mu<Nd;mu++){
//Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction.
//mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
for(int mu=0;mu<Nd-1;mu++){
LatticeCoordinate(coor,mu);
U = PeekIndex<LorentzIndex>(Umu,mu);
Uconj = conjugate(U);
// Implement the isospin rotation sign on the boundary between f=1 and f=0
// This phase could come from a simple bc 1,1,-1,1 ..
int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
if ( Params.twists[mu] ) {
@ -260,6 +288,38 @@ public:
});
}
}
{ //periodic / antiperiodic temporal BCs
int mu = Nd-1;
int L = GaugeGrid->GlobalDimensions()[mu];
int Lmu = L - 1;
LatticeCoordinate(coor, mu);
U = PeekIndex<LorentzIndex>(Umu, mu); //Get t-directed links
GaugeLinkField *Upoke = &U;
if(Params.twists[mu]){ //antiperiodic
Utmp = where(coor == Lmu, -U, U);
Upoke = &Utmp;
}
Uconj = conjugate(*Upoke); //second flavor interacts with conjugate links
pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu);
//Get the barrel-shifted field
Utmp = adj(Cshift(U, mu, -1)); //is a forward shift!
Upoke = &Utmp;
if(Params.twists[mu]){
U = where(coor == 0, -Utmp, Utmp); //boundary phase
Upoke = &U;
}
Uconj = conjugate(*Upoke);
pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4);
}
}
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A, int mu) {

34
tests/core/Test_gparity.cc
View File

@ -55,13 +55,17 @@ static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIM
int main (int argc, char ** argv)
{
int nu = 0;
int tbc_aprd = 0; //use antiperiodic BCs in the time direction?
Grid_init(&argc,&argv);
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "--Gparity-dir"){
std::stringstream ss; ss << argv[i+1]; ss >> nu;
std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl;
}else if(std::string(argv[i]) == "--Tbc-APRD"){
tbc_aprd = 1;
std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl;
}
}
@ -161,7 +165,12 @@ int main (int argc, char ** argv)
RealD mass=0.0;
RealD M5=1.8;
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS);
//Standard Dirac op
AcceleratorVector<Complex,4> bc_std(Nd, 1.0);
if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC
StandardDiracOp::ImplParams std_params(bc_std);
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS, std_params);
StandardFermionField src_o_1f(FrbGrid_1f);
StandardFermionField result_o_1f(FrbGrid_1f);
@ -172,9 +181,11 @@ int main (int argc, char ** argv)
ConjugateGradient<StandardFermionField> CG(1.0e-8,10000);
CG(HermOpEO,src_o_1f,result_o_1f);
// const int nu = 3;
//Gparity Dirac op
std::vector<int> twists(Nd,0);
twists[nu] = 1;
if(tbc_aprd) twists[Nd-1] = 1;
GparityDiracOp::ImplParams params;
params.twists = twists;
GparityDiracOp GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5 DOP_PARAMS,params);
@ -271,8 +282,11 @@ int main (int argc, char ** argv)
std::cout << "2f cb "<<result_o_2f.Checkerboard()<<std::endl;
std::cout << "1f cb "<<result_o_1f.Checkerboard()<<std::endl;
std::cout << " result norms " <<norm2(result_o_2f)<<" " <<norm2(result_o_1f)<<std::endl;
//Compare norms
std::cout << " result norms 2f: " <<norm2(result_o_2f)<<" 1f: " <<norm2(result_o_1f)<<std::endl;
//Take the 2f solution and convert into the corresponding 1f solution (odd cb only)
StandardFermionField res0o (FrbGrid_2f);
StandardFermionField res1o (FrbGrid_2f);
StandardFermionField res0 (FGrid_2f);
@ -281,14 +295,15 @@ int main (int argc, char ** argv)
res0=Zero();
res1=Zero();
res0o = PeekIndex<0>(result_o_2f,0);
res1o = PeekIndex<0>(result_o_2f,1);
res0o = PeekIndex<0>(result_o_2f,0); //flavor 0, odd cb
res1o = PeekIndex<0>(result_o_2f,1); //flavor 1, odd cb
std::cout << "res cb "<<res0o.Checkerboard()<<std::endl;
std::cout << "res cb "<<res1o.Checkerboard()<<std::endl;
setCheckerboard(res0,res0o);
setCheckerboard(res1,res1o);
//poke odd onto non-cb field
setCheckerboard(res0,res0o);
setCheckerboard(res1,res1o);
StandardFermionField replica (FGrid_1f);
StandardFermionField replica0(FGrid_1f);
@ -296,12 +311,13 @@ int main (int argc, char ** argv)
Replicate(res0,replica0);
Replicate(res1,replica1);
//2nd half of doubled lattice has f=1
replica = where( xcoor_1f5 >= Integer(L), replica1,replica0 );
replica0 = Zero();
setCheckerboard(replica0,result_o_1f);
std::cout << "Norm2 solutions is " <<norm2(replica)<<" "<< norm2(replica0)<<std::endl;
std::cout << "Norm2 solutions 1f reconstructed from 2f: " <<norm2(replica)<<" Actual 1f: "<< norm2(replica0)<<std::endl;
replica = replica - replica0;