portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-24 12:45:56 +01:00
Code Releases Activity

Starting the DDHMC filter. This will be seriously subject to tuning as force in

Browse Source 
Omega and force in Omega bar appears to have different structure
This commit is contained in:
Peter Boyle 2021-05-14 11:26:47 -04:00
parent 5b52f29b2f
commit def51267e9
1 changed files with 36 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -28,6 +28,7 @@ directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
NAMESPACE_BEGIN(Grid);
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
// DDHMC filter with sub-block size B[mu] // DDHMC filter with sub-block size B[mu]
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
@ -43,37 +44,60 @@ struct DDHMCFilter: public MomentumFilterBase<MomentaField>
typedef Lattice<ScalarType> LatticeScalarType; typedef Lattice<ScalarType> LatticeScalarType;
Coordinate Block; Coordinate Block;
int Width;
// Could also zero links in plane like Luscher advocates. // Could also zero links in plane like Luscher advocates.
DDHMCFilter(const Coordinate &_Block): Block(_Block){} DDHMCFilter(const Coordinate &_Block,int _Width): Block(_Block), Width(_Width)
{
assert( (Width==0) || (Width==1) || (Width==2) );
}
void applyFilter(MomentaField &P) const override void applyFilter(MomentaField &P) const override
{ {
GridBase *grid = P.Grid(); GridBase *grid = P.Grid();
LatticeScalarType mask_mu(grid);
LatticeLorentzScalarType mask(grid); LatticeColourMatrix zz(grid); zz = Zero();
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
// Zero strictly links crossing between domains // Zero strictly links crossing between domains
// Luscher also zeroes links in plane of domain boundaries // Luscher also zeroes links in plane of domain boundaries
// Keeping interior only. This prevents force from plaquettes // Keeping interior only. This prevents force from plaquettes
// crossing domains and keeps whole MD trajectory local. // crossing domains and keeps whole MD trajectory local.
// Step further than I want to go. // This is the case Width=1.
// Width = 2 should also decouple rectangle force.
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
ScalarType zz = scalar_type(0.0); Coordinate Global=grid->GlobalDimensions();
ScalarType one= scalar_type(1.0);
LatticeInteger coor(grid); LatticeInteger coor(grid);
for(int mu=0;mu<Nd;mu++) { for(int mu=0;mu<Nd;mu++) {
if ( Block[mu] <= Global[mu] ) {
LatticeCoordinate(coor,mu); LatticeCoordinate(coor,mu);
auto P_mu = PeekIndex<LorentzIndex>(P, mu);
P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_mu);
/* if(Width>=2) {
P_mu = where(mod(coor,Block[mu])==Integer(0),zz,P_mu);
P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zz,P_mu);
}*/
PokeIndex<LorentzIndex>(P, P_mu, mu);
mask_mu = where(mod(coor,Block[mu])==Block[mu]-1,zz,one); for(int nu=0;nu<Nd;nu++){
if ( (mu!=nu) && Block[mu] <= Global[mu]){
PokeIndex<LorentzIndex>(mask, mask_mu, mu); auto P_nu = PeekIndex<LorentzIndex>(P, nu);
if ( Width>=1 ){
P_nu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_nu);
P_nu = where(mod(coor,Block[mu])==Integer(0) ,zz,P_nu);
}
/* if ( Width>=2 ){
P_nu = where(mod(coor,Block[mu])==Integer(Block[mu]-2),zz,P_nu);
P_nu = where(mod(coor,Block[mu])==Integer(1) ,zz,P_nu);
}*/
PokeIndex<LorentzIndex>(P, P_nu, nu);
}
}
}
} }
} }
}; };
NAMESPACE_END(Grid);