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

General detection for arbitrary domains.

Browse Source 
Can simplify and make specific if performance matters
This commit is contained in:
Quadro 2021-06-09 12:54:43 -04:00
parent f9cda24781
commit 6dcaed621c
1 changed files with 47 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -32,37 +32,67 @@ directory
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
#pragma once #pragma once
#include <Grid/qcd/action/momentum/Domains.h>
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
template<typename MomentaField> template<typename MomentaField>
struct DirichletFilter: public MomentumFilterBase<MomentaField> struct DirichletFilter: public MomentumFilterBase<MomentaField>
{ {
Coordinate Block; Coordinate Block;
DirichletFilter(const Coordinate &_Block): Block(_Block){} DirichletFilter(const Coordinate &_Block): Block(_Block) {}
void applyFilter(MomentaField &P) const override // Edge detect using domain projectors
void applyFilter (MomentaField &U) const override
{ {
GridBase *grid = P.Grid(); DomainDecomposition Domains(Block);
GridBase *grid = U.Grid();
LatticeInteger coor(grid);
LatticeInteger face(grid);
LatticeInteger one(grid); one = 1;
LatticeInteger zero(grid); zero = 0;
LatticeInteger omega(grid);
LatticeInteger omegabar(grid);
LatticeInteger tmp(grid);
//////////////////////////////////////////////////// omega=one; Domains.ProjectDomain(omega,0);
// Zero strictly links crossing between domains omegabar=one; Domains.ProjectDomain(omegabar,1);
////////////////////////////////////////////////////
LatticeInteger coor(grid); LatticeInteger nface(grid); nface=Zero();
typedef decltype(PeekIndex<LorentzIndex>(P,0)) MatrixType;
MatrixType zz(grid); zz = Zero(); MomentaField projected(grid); projected=Zero();
typedef decltype(PeekIndex<LorentzIndex>(U,0)) MomentaLinkField;
MomentaLinkField Umu(grid);
MomentaLinkField zz(grid); zz=Zero();
int dims = grid->Nd();
Coordinate Global=grid->GlobalDimensions(); Coordinate Global=grid->GlobalDimensions();
for(int mu=0;mu<Nd;mu++) { assert(dims==Nd);
if ( (Block[mu] <= Global[mu]) && (Block[mu]>1) ) {
// If costly could provide Grid earlier and precompute masks
LatticeCoordinate(coor,mu);
auto P_mu = PeekIndex<LorentzIndex>(P, mu); for(int mu=0;mu<Nd;mu++){
P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_mu);
PokeIndex<LorentzIndex>(P, P_mu, mu); if ( Block[mu]!=0 ) {
Umu = PeekIndex<LorentzIndex>(U,mu);
// Upper face
tmp = Cshift(omegabar,mu,1);
tmp = tmp + omega;
face = where(tmp == Integer(2),one,zero );
tmp = Cshift(omega,mu,1);
tmp = tmp + omegabar;
face = where(tmp == Integer(2),one,face );
Umu = where(face,zz,Umu);
PokeIndex<LorentzIndex>(U, Umu, mu);
} }
} }
} }
}; };
NAMESPACE_END(Grid); NAMESPACE_END(Grid);