portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Cleaner

Browse Source
This commit is contained in:
Quadro
2021-06-08 21:50:43 -04:00
parent cd5e3fbd82
commit f9cda24781
1 changed files with 4 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files

62
Grid/qcd/action/momentum/DDHMCfilter.h
View File

@ -36,71 +36,17 @@ NAMESPACE_BEGIN(Grid);
template<typename MomentaField> template<typename MomentaField>
struct DDHMCFilter: public MomentumFilterBase<MomentaField> struct DDHMCFilter: public MomentumFilterBase<MomentaField>
{ {
// typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
// typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
// typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
// typedef iScalar<iScalar<iScalar<vector_type> > > ScalarType; //complex phase for each site
// typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
// typedef Lattice<ScalarType> LatticeScalarType;
Coordinate Block; Coordinate Block;
int Width; int Width;
// Could also zero links in plane like Luscher advocates.
DDHMCFilter(const Coordinate &_Block,int _Width): Block(_Block), Width(_Width) DDHMCFilter(const Coordinate &_Block): Block(_Block) {}
{
assert( (Width==0) || (Width==1) || (Width==2) );
}
void applyFilter(MomentaField &P) const override void applyFilter(MomentaField &P) const override
{ {
GridBase *grid = P.Grid(); DomainDecomposition Domains(Block);
Domains.ProjectDDHMC(P);
LatticeColourMatrix zz(grid); zz = Zero();
MomentaField zzz(grid); zzz = Zero();
////////////////////////////////////////////////////
// Zero strictly links crossing between domains
// Luscher also zeroes links in plane of domain boundaries
// Keeping interior only. This prevents force from plaquettes
// crossing domains and keeps whole MD trajectory local.
// This is the case Width=1.
// Width = 2 should also decouple rectangle force.
////////////////////////////////////////////////////
Coordinate Global=grid->GlobalDimensions();
LatticeInteger coor(grid);
for(int mu=0;mu<Nd;mu++) {
#define PROJECT_DOMAIN
if ( (Block[mu] <= Global[mu])&&(Block[mu]>0) ) {
LatticeCoordinate(coor,mu);
#if 1
P = where(mod(coor,Block[mu])==Integer(Block[mu]-3),zzz,P); //width 4
P = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zzz,P); //width 4
P = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zzz,P); //width 2
P = where(mod(coor,Block[mu])==Integer(0),zzz,P); //width 2
P = where(mod(coor,Block[mu])==Integer(1),zzz,P); //width 4
auto P_mu = PeekIndex<LorentzIndex>(P,mu);
P = where(mod(coor,Block[mu])==Integer(2),zzz,P); //width 6
PokeIndex<LorentzIndex>(P, P_mu, mu);
#else
P = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zzz,P); //width 2
auto P_mu = PeekIndex<LorentzIndex>(P,mu);
P = where(mod(coor,Block[mu])==Integer(0),zzz,P); //width 6
PokeIndex<LorentzIndex>(P, P_mu, mu);
#endif
}
}
#ifdef PROJECT_DOMAIN
LatticeInteger domaincb(grid); domaincb=Zero();
for(int d=0;d<Nd;d++){
LatticeCoordinate(coor,d);
domaincb = domaincb + div(coor,Block[d]);
}
P = where(mod(domaincb,2)==Integer(1),P,zzz);
#endif
} }
}; };
NAMESPACE_END(Grid); NAMESPACE_END(Grid);