mirror of
https://github.com/paboyle/Grid.git
synced 2026-01-07 02:19:34 +00:00
Somewhat better wrapped support for Icosahedral
This commit is contained in:
Peter Boyle
@@ -391,7 +391,7 @@ template<class vobj> std::ostream& operator<< (std::ostream& stream, const Latti
|
||||
stream<<"]\t";
|
||||
stream<<ss<<std::endl;
|
||||
}
|
||||
if ( o.Grid()->isIcosahedral() ) {
|
||||
if ( o.Grid()->isIcosahedralVertex() ) {
|
||||
uint64_t psites=1;
|
||||
Coordinate perpdims;
|
||||
for(int d=2;d<o.Grid()->_ndimension-1;d++){
|
||||
|
||||
@@ -281,6 +281,131 @@ void pokePole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
|
||||
return;
|
||||
};
|
||||
|
||||
|
||||
template<class vobj,class sobj>
|
||||
void peekLocalPole(sobj &s,const Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
|
||||
{
|
||||
s=Zero();
|
||||
|
||||
GridBase *grid=l.Grid();
|
||||
|
||||
assert(grid->isIcosahedral());
|
||||
assert(grid->isIcosahedralVertex());
|
||||
|
||||
int Nsimd = grid->Nsimd();
|
||||
|
||||
int rank;
|
||||
|
||||
int Ndm1 = grid->_ndimension-1;
|
||||
Coordinate pgrid = grid->ProcessorGrid();
|
||||
const int xdim=0;
|
||||
const int ydim=1;
|
||||
const int pdim=Ndm1;
|
||||
|
||||
int64_t pole_osite;
|
||||
int64_t pole_isite;
|
||||
Coordinate rdims;
|
||||
Coordinate idims;
|
||||
Coordinate ocoor;
|
||||
Coordinate icoor;
|
||||
// Coordinate pcoor(grid->_ndimension);
|
||||
for(int d=2;d<Ndm1;d++){
|
||||
int dd=d-2;
|
||||
rdims.push_back(grid->_rdimensions[d]);
|
||||
idims.push_back(grid->_simd_layout[d]);
|
||||
icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
|
||||
ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
|
||||
// pcoor[d] = orthog[dd]/grid->_ldimensions[d];
|
||||
}
|
||||
Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
|
||||
Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
|
||||
|
||||
int64_t osite;
|
||||
if(isNorth == North){
|
||||
// pcoor[xdim] = 0;
|
||||
// pcoor[ydim] = pgrid[ydim]-1;
|
||||
// pcoor[Ndm1] = pgrid[Ndm1]-1;
|
||||
osite = pole_osite + grid->NorthPoleOsite();
|
||||
assert(grid->ownsNorthPole());
|
||||
} else {
|
||||
// pcoor[xdim] = pgrid[xdim]-1;
|
||||
// pcoor[ydim] = 0;
|
||||
// pcoor[Ndm1] = 0;
|
||||
osite = pole_osite + grid->SouthPoleOsite();
|
||||
assert(grid->ownsSouthPole());
|
||||
}
|
||||
|
||||
ExtractBuffer<sobj> buf(Nsimd);
|
||||
autoView( l_v , l, CpuWrite);
|
||||
extract(l_v[osite],buf);
|
||||
s = buf[pole_isite];
|
||||
|
||||
return;
|
||||
};
|
||||
template<class vobj,class sobj>
|
||||
void pokeLocalPole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
|
||||
{
|
||||
GridBase *grid=l.Grid();
|
||||
|
||||
assert(grid->isIcosahedral());
|
||||
assert(grid->isIcosahedralVertex());
|
||||
|
||||
int Nsimd = grid->Nsimd();
|
||||
int rank;
|
||||
int Ndm1 = grid->_ndimension-1;
|
||||
|
||||
const int xdim=0;
|
||||
const int ydim=1;
|
||||
const int pdim=Ndm1;
|
||||
|
||||
int64_t pole_osite;
|
||||
int64_t pole_isite;
|
||||
Coordinate rdims;
|
||||
Coordinate idims;
|
||||
Coordinate ocoor;
|
||||
Coordinate icoor;
|
||||
// Coordinate pcoor(grid->_ndimension,0);
|
||||
for(int d=2;d<Ndm1;d++){
|
||||
int dd = d-2;
|
||||
rdims.push_back(grid->_rdimensions[d]);
|
||||
idims.push_back(grid->_simd_layout[d]);
|
||||
icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
|
||||
ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
|
||||
// pcoor[d] = orthog[dd]/grid->_ldimensions[d];
|
||||
|
||||
int o = orthog[dd];
|
||||
int r = grid->_rdimensions[d];
|
||||
int omr = o % r;
|
||||
}
|
||||
Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
|
||||
Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
|
||||
|
||||
int64_t osite;
|
||||
int insert=0;
|
||||
if(isNorth ==North){
|
||||
// pcoor[xdim] = 0;
|
||||
// pcoor[ydim] = pgrid[ydim]-1;
|
||||
// pcoor[Ndm1] = pgrid[Ndm1]-1;
|
||||
osite = pole_osite + grid->NorthPoleOsite();
|
||||
assert(grid->ownsNorthPole());
|
||||
} else {
|
||||
// pcoor[xdim] = pgrid[xdim]-1;
|
||||
// pcoor[ydim] = 0;
|
||||
// pcoor[Ndm1] = 0;
|
||||
osite = pole_osite + grid->SouthPoleOsite();
|
||||
assert(grid->ownsSouthPole());
|
||||
}
|
||||
|
||||
// extract-modify-merge cycle is easiest way and this is not perf critical
|
||||
ExtractBuffer<sobj> buf(Nsimd);
|
||||
autoView( l_v , l, CpuWrite);
|
||||
extract(l_v[osite],buf);
|
||||
buf[pole_isite] = s;
|
||||
merge(l_v[osite],buf);
|
||||
|
||||
return;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////
|
||||
// Peek a scalar object from the SIMD array
|
||||
//////////////////////////////////////////////////////////
|
||||
|
||||
@@ -48,31 +48,45 @@ NAMESPACE_BEGIN(Grid);
|
||||
//////////////////////////////////////////////////////////////
|
||||
inline int RNGfillable(GridBase *coarse,GridBase *fine)
|
||||
{
|
||||
if ( coarse == fine ) return 1;
|
||||
|
||||
int rngdims = coarse->_ndimension;
|
||||
|
||||
// trivially extended in higher dims, with locality guaranteeing RNG state is local to node
|
||||
int lowerdims = fine->_ndimension - coarse->_ndimension;
|
||||
assert(lowerdims >= 0);
|
||||
for(int d=0;d<lowerdims;d++){
|
||||
assert(fine->_simd_layout[d]==1);
|
||||
assert(fine->_processors[d]==1);
|
||||
if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
|
||||
|
||||
if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
|
||||
assert(fine->Nd()==coarse->Nd());
|
||||
for(int d=0;d<fine->Nd();d++){
|
||||
assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
{
|
||||
|
||||
int rngdims = coarse->_ndimension;
|
||||
|
||||
int multiplicity=1;
|
||||
for(int d=0;d<lowerdims;d++){
|
||||
multiplicity=multiplicity*fine->_rdimensions[d];
|
||||
}
|
||||
// local and global volumes subdivide cleanly after SIMDization
|
||||
for(int d=0;d<rngdims;d++){
|
||||
int fd= d+lowerdims;
|
||||
assert(coarse->_processors[d] == fine->_processors[fd]);
|
||||
assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
|
||||
assert(((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d])==fine->_rdimensions[fd]);
|
||||
// trivially extended in higher dims, with locality guaranteeing RNG state is local to node
|
||||
int lowerdims = fine->_ndimension - coarse->_ndimension;
|
||||
assert(lowerdims >= 0);
|
||||
for(int d=0;d<lowerdims;d++){
|
||||
assert(fine->_simd_layout[d]==1);
|
||||
assert(fine->_processors[d]==1);
|
||||
}
|
||||
|
||||
multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d];
|
||||
int multiplicity=1;
|
||||
for(int d=0;d<lowerdims;d++){
|
||||
multiplicity=multiplicity*fine->_rdimensions[d];
|
||||
}
|
||||
// local and global volumes subdivide cleanly after SIMDization
|
||||
for(int d=0;d<rngdims;d++){
|
||||
int fd= d+lowerdims;
|
||||
assert(coarse->_processors[d] == fine->_processors[fd]);
|
||||
assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
|
||||
assert(((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d])==fine->_rdimensions[fd]);
|
||||
|
||||
multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d];
|
||||
}
|
||||
return multiplicity;
|
||||
}
|
||||
return multiplicity;
|
||||
}
|
||||
|
||||
|
||||
@@ -80,6 +94,19 @@ inline int RNGfillable(GridBase *coarse,GridBase *fine)
|
||||
// this function is necessary for the LS vectorised field
|
||||
inline int RNGfillable_general(GridBase *coarse,GridBase *fine)
|
||||
{
|
||||
|
||||
if ( coarse == fine ) return 1;
|
||||
|
||||
if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
|
||||
|
||||
if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
|
||||
assert(fine->Nd()==coarse->Nd());
|
||||
for(int d=0;d<fine->Nd();d++){
|
||||
assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
int rngdims = coarse->_ndimension;
|
||||
|
||||
// trivially extended in higher dims, with locality guaranteeing RNG state is local to node
|
||||
@@ -352,12 +379,12 @@ private:
|
||||
public:
|
||||
GridBase *Grid(void) const { return _grid; }
|
||||
int generator_idx(int os,int is) {
|
||||
return is*_grid->oSites()+os;
|
||||
return (is*_grid->CartesianOsites()+os)%_grid->lSites(); // On the pole sites wrap back to normal generators; Icosahedral hack
|
||||
}
|
||||
|
||||
GridParallelRNG(GridBase *grid) : GridRNGbase() {
|
||||
_grid = grid;
|
||||
_vol =_grid->iSites()*_grid->oSites();
|
||||
_vol =_grid->lSites();
|
||||
|
||||
_generators.resize(_vol);
|
||||
_uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
|
||||
@@ -381,7 +408,7 @@ public:
|
||||
|
||||
int multiplicity = RNGfillable_general(_grid, l.Grid()); // l has finer or same grid
|
||||
int Nsimd = _grid->Nsimd(); // guaranteed to be the same for l.Grid() too
|
||||
int osites = _grid->oSites(); // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity
|
||||
int osites = _grid->CartesianOsites(); // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity, except on Icosahedral
|
||||
int words = sizeof(scalar_object) / sizeof(scalar_type);
|
||||
|
||||
autoView(l_v, l, CpuWrite);
|
||||
@@ -402,8 +429,27 @@ public:
|
||||
// merge into SIMD lanes, FIXME suboptimal implementation
|
||||
merge(l_v[sm], buf);
|
||||
}
|
||||
});
|
||||
// });
|
||||
});
|
||||
|
||||
/*
|
||||
* Fill in the poles for an Icosahedral vertex mesh
|
||||
*/
|
||||
if (l.Grid()->isIcosahedralVertex()) {
|
||||
int64_t pole_sites=l.Grid()->NorthPoleOsites()+l.Grid()->SouthPoleOsites();
|
||||
int64_t pole_base =l.Grid()->CartesianOsites();
|
||||
|
||||
ExtractBuffer<scalar_object> buf(Nsimd);
|
||||
for (int m = 0; m < pole_sites; m++) { // Draw from same generator multiplicity times
|
||||
for (int si = 0; si < Nsimd; si++) {
|
||||
int gdx = 0;
|
||||
scalar_type *pointer = (scalar_type *)&buf[si];
|
||||
dist[gdx].reset();
|
||||
for (int idx = 0; idx < words; idx++)
|
||||
fillScalar(pointer[idx], dist[gdx], _generators[gdx]);
|
||||
}
|
||||
merge(l_v[pole_base+m], buf);
|
||||
}
|
||||
}
|
||||
|
||||
_time_counter += usecond()- inner_time_counter;
|
||||
}
|
||||
|
||||
@@ -584,6 +584,8 @@ public:
|
||||
assert(grid->isIcosahedralVertex());
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* For gauge action implementation
|
||||
*/
|
||||
|
||||
@@ -1,5 +1,6 @@
|
||||
MPICXX=mpicxx CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=clang++ ../../configure \
|
||||
--enable-simd=GEN \
|
||||
--enable-Nc=1 \
|
||||
--enable-comms=mpi-auto \
|
||||
--enable-unified=yes \
|
||||
--prefix $HOME/QCD/GridInstall \
|
||||
|
||||
@@ -32,80 +32,63 @@ using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
const int MyNd=3;
|
||||
template<typename vtype> using iIcosahedralLorentzComplex = iVector<iScalar<iScalar<vtype> >, MyNd+1 > ;
|
||||
|
||||
typedef iIcosahedralLorentzComplex<ComplexD > IcosahedralLorentzComplexD;
|
||||
typedef iIcosahedralLorentzComplex<vComplexD> vIcosahedralLorentzComplexD;
|
||||
typedef Lattice<vIcosahedralLorentzComplexD> LatticeIcosahedralLorentzComplexD;
|
||||
template<typename vtype> using iIcosahedralLorentzComplex = iVector<iScalar<iScalar<vtype> >, MyNd+1 > ;
|
||||
template<typename vtype> using iIcosahedralLorentzColourMatrix = iVector<iScalar<iMatrix<vtype,Nc> >, MyNd+1 > ;
|
||||
template<typename vtype> using iIcosahedralColourMatrix = iScalar<iScalar<iMatrix<vtype,Nc> > > ;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
typedef iIcosahedralLorentzComplex<Complex > IcosahedralLorentzComplex;
|
||||
typedef iIcosahedralLorentzComplex<vComplex> vIcosahedralLorentzComplex;
|
||||
typedef Lattice<vIcosahedralLorentzComplex> LatticeIcosahedralLorentzComplex;
|
||||
|
||||
typedef iIcosahedralLorentzColourMatrix<Complex > IcosahedralLorentzColourMatrix;
|
||||
typedef iIcosahedralLorentzColourMatrix<vComplex> vIcosahedralLorentzColourMatrix;
|
||||
typedef Lattice<vIcosahedralLorentzColourMatrix> LatticeIcosahedralLorentzColourMatrix;
|
||||
|
||||
typedef iIcosahedralColourMatrix<Complex > IcosahedralColourMatrix;
|
||||
typedef iIcosahedralColourMatrix<vComplex> vIcosahedralColourMatrix;
|
||||
typedef Lattice<vIcosahedralColourMatrix> LatticeIcosahedralColourMatrix;
|
||||
|
||||
|
||||
class IcosahedralGimpl
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int L=8;
|
||||
const int Npatch = IcosahedralPatches;
|
||||
|
||||
// Put SIMD all in time direction
|
||||
Coordinate latt_size = GridDefaultLatt();
|
||||
Coordinate simd_layout({1,1,vComplexD::Nsimd(),1});
|
||||
Coordinate mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::cout << GridLogMessage << " mpi "<<mpi_layout<<std::endl;
|
||||
std::cout << GridLogMessage << " simd "<<simd_layout<<std::endl;
|
||||
std::cout << GridLogMessage << " latt "<<latt_size<<std::endl;
|
||||
GridCartesianCrossIcosahedron EdgeGrid(latt_size,simd_layout,mpi_layout,IcosahedralEdges);
|
||||
std::cout << GridLogMessage << " Created edge grid "<<std::endl;
|
||||
GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
|
||||
|
||||
std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
|
||||
LatticeIcosahedralLorentzComplexD Umu(&EdgeGrid);
|
||||
LatticeComplex Phi(&VertexGrid);
|
||||
std::cout << GridLogMessage << " Created two fields "<<std::endl;
|
||||
|
||||
Phi = Zero();
|
||||
Umu = Zero();
|
||||
std::cout << GridLogMessage << " Zeroed two fields "<<std::endl;
|
||||
|
||||
ComplexD one (1.0);
|
||||
Phi = one;
|
||||
Umu = one;
|
||||
public:
|
||||
typedef LatticeIcosahedralLorentzColourMatrix GaugeField;
|
||||
typedef LatticeIcosahedralColourMatrix GaugeLinkField;
|
||||
typedef LatticeComplex ComplexField;
|
||||
};
|
||||
|
||||
std::cout << GridLogMessage << " V = "<<norm2(Phi)<<std::endl;
|
||||
std::cout << GridLogMessage << " Expect "<<latt_size[0]*latt_size[1]*latt_size[2]*10+2*latt_size[2]<<std::endl;
|
||||
template< class Gimpl>
|
||||
class IcosahedralEdgeSupport {
|
||||
public:
|
||||
// Move to inherit types macro as before
|
||||
typedef typename Gimpl::GaugeField GaugeField;
|
||||
typedef typename Gimpl::GaugeLinkField GaugeLinkField;
|
||||
typedef typename Gimpl::ComplexField ComplexField;
|
||||
//
|
||||
GridBase *VertexGrid;
|
||||
GridBase *EdgeGrid;
|
||||
IcosahedralStencil FaceStencil;
|
||||
|
||||
std::cout << GridLogMessage << " E = "<<norm2(Umu)<<std::endl;
|
||||
std::cout << GridLogMessage << " Expect "<<latt_size[0]*latt_size[1]*latt_size[2]*10*4<<std::endl;
|
||||
|
||||
// std::cout << " Umu "<<Umu<<std::endl;
|
||||
// std::cout << " Phi "<<Phi<<std::endl;
|
||||
LatticePole(Phi,South);
|
||||
// std::cout << " Phi South Pole set\n"<<Phi<<std::endl;
|
||||
|
||||
LatticePole(Phi,North);
|
||||
// std::cout << " Phi North Pole set\n"<<Phi<<std::endl;
|
||||
|
||||
for(int mu=0;mu<VertexGrid._ndimension;mu++){
|
||||
std::cout << " Calling lattice coordinate mu="<<mu<<std::endl;
|
||||
LatticeCoordinate(Phi,mu);
|
||||
// std::cout << " Phi coor mu="<<mu<<"\n"<<Phi<<std::endl;
|
||||
}
|
||||
IcosahedralStencil stencil(&EdgeGrid);
|
||||
|
||||
stencil.TestGeometry();
|
||||
|
||||
std::cout << "Creating face stencil"<<std::endl;
|
||||
|
||||
stencil.FaceStencil();
|
||||
Umu=one;
|
||||
LatticeComplex plaq1(&EdgeGrid);
|
||||
LatticeComplex plaq2(&EdgeGrid);
|
||||
|
||||
IcosahedralEdgeSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid)
|
||||
: FaceStencil (EdgeGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
|
||||
{
|
||||
autoView(Umu_v,Umu,AcceleratorRead);
|
||||
autoView(plaq1_v,plaq1,AcceleratorWrite);
|
||||
autoView(plaq2_v,plaq2,AcceleratorWrite);
|
||||
autoView(stencil_v,stencil,AcceleratorRead);
|
||||
accelerator_for(ss,EdgeGrid.oSites(),vComplexD::Nsimd(),{
|
||||
FaceStencil.FaceStencil();
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////
|
||||
// Fixme: will need a version of "Gimpl" and a wrapper class following "WilsonLoops" style.
|
||||
// Gauge Link field GT is the gauge transform and lives on the VERTEX field
|
||||
////////////////////////////////////////////////////////////////////////////////////
|
||||
void ForwardTriangles(GaugeField &Umu,LatticeComplex &plaq1,LatticeComplex &plaq2)
|
||||
{
|
||||
{
|
||||
autoView(Umu_v,Umu,AcceleratorRead);
|
||||
autoView(plaq1_v,plaq1,AcceleratorWrite);
|
||||
autoView(plaq2_v,plaq2,AcceleratorWrite);
|
||||
autoView(stencil_v,FaceStencil,AcceleratorRead);
|
||||
|
||||
accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
|
||||
|
||||
const int x = IcosahedronPatchX;
|
||||
const int y = IcosahedronPatchY;
|
||||
@@ -141,8 +124,213 @@ int main (int argc, char ** argv)
|
||||
coalescedWrite(plaq2_v[ss](),trace(Ly*adj(Ld)*L2 ) );
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
void GaugeTransform(GaugeLinkField >, GaugeField &Umu)
|
||||
{
|
||||
assert(gt.Grid()==VertexGrid);
|
||||
assert(Umu.Grid()==EdgeGrid);
|
||||
|
||||
GridBase * vgrid = VertexGrid;
|
||||
GridBase * grid = EdgeGrid;
|
||||
|
||||
int osites = grid->oSites();
|
||||
|
||||
uint64_t cart_sites = grid->CartesianOsites();
|
||||
uint64_t Npole_sites = grid->NorthPoleOsites();
|
||||
uint64_t Spole_sites = grid->SouthPoleOsites();
|
||||
Coordinate pcoor = grid->ThisProcessorCoor();
|
||||
Coordinate pgrid = grid->ProcessorGrid();
|
||||
/*
|
||||
* resize the stencil entries array and set npoints
|
||||
*/
|
||||
autoView(g_v,gt,CpuRead);
|
||||
autoView(Umu_v,Umu,CpuWrite);
|
||||
for(uint64_t site=0;site<cart_sites; site ++) {
|
||||
|
||||
Coordinate Coor;
|
||||
Coordinate NbrCoor;
|
||||
|
||||
int nd = grid->Nd();
|
||||
int L = grid->LocalDimensions()[0];
|
||||
|
||||
////////////////////////////////////////////////
|
||||
// Outer index of neighbour Offset calculation
|
||||
////////////////////////////////////////////////
|
||||
grid->oCoorFromOindex(Coor,site);
|
||||
NbrCoor = Coor;
|
||||
assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]);
|
||||
assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims
|
||||
assert( grid->_simd_layout[1]==1);
|
||||
assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims
|
||||
assert( grid->_processors[1]==1);
|
||||
|
||||
int Patch = Coor[nd-1];
|
||||
int HemiPatch = Patch%HemiPatches;
|
||||
int north = Patch/HemiPatches;
|
||||
int south = 1-north;
|
||||
int isPoleY;
|
||||
int isPoleX;
|
||||
|
||||
assert(Patch<IcosahedralPatches);
|
||||
assert((north==1)||(south==1));
|
||||
|
||||
Coordinate XpCoor;
|
||||
Coordinate YpCoor;
|
||||
Coordinate DpCoor;
|
||||
|
||||
FaceStencil.GetNbrForPlusDiagonal(grid,Coor,DpCoor);
|
||||
FaceStencil.GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
|
||||
FaceStencil.GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
|
||||
|
||||
int XpHemiPatch = XpCoor[nd-1]%HemiPatches;
|
||||
int XpHemisphere = XpCoor[nd-1]/HemiPatches;
|
||||
|
||||
int DpPatch = DpCoor[nd-1];
|
||||
int DpHemiPatch = DpCoor[nd-1]%HemiPatches;
|
||||
int DpHemisphere = DpCoor[nd-1]/HemiPatches;
|
||||
|
||||
// Work out the pole_osite
|
||||
Coordinate rdims;
|
||||
Coordinate ocoor;
|
||||
int64_t pole_osite;
|
||||
int Ndm1 = grid->Nd()-1;
|
||||
for(int d=2;d<Ndm1;d++){
|
||||
int dd=d-2;
|
||||
rdims.push_back(grid->_rdimensions[d]);
|
||||
ocoor.push_back(Coor[d]%grid->_rdimensions[d]);
|
||||
}
|
||||
Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
|
||||
|
||||
uint64_t xp_idx;
|
||||
uint64_t yp_idx;
|
||||
uint64_t dp_idx;
|
||||
if ( isPoleX ) {
|
||||
assert(vgrid->ownsSouthPole());
|
||||
xp_idx = pole_osite + grid->SouthPoleOsite();
|
||||
} else {
|
||||
xp_idx = grid->oIndex(XpCoor);
|
||||
}
|
||||
if ( isPoleY ) {
|
||||
assert(vgrid->ownsNorthPole());
|
||||
yp_idx = pole_osite + grid->NorthPoleOsite();
|
||||
} else {
|
||||
yp_idx = grid->oIndex(YpCoor);
|
||||
}
|
||||
dp_idx = grid->oIndex(DpCoor);
|
||||
|
||||
auto g = g_v(site)();
|
||||
auto gx = g_v(xp_idx)();
|
||||
auto gy = g_v(yp_idx)();
|
||||
auto gd = g_v(dp_idx)();
|
||||
|
||||
auto lx = Umu_v(site)(IcosahedronPatchX);
|
||||
auto ly = Umu_v(site)(IcosahedronPatchY);
|
||||
auto ld = Umu_v(site)(IcosahedronPatchDiagonal);
|
||||
|
||||
lx = g*lx*adj(gx);
|
||||
ly = g*ly*adj(gy);
|
||||
ld = g*ld*adj(gd);
|
||||
|
||||
coalescedWrite(Umu_v[site](IcosahedronPatchX),lx);
|
||||
coalescedWrite(Umu_v[site](IcosahedronPatchY),ly);
|
||||
coalescedWrite(Umu_v[site](IcosahedronPatchDiagonal),ld);
|
||||
};
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int L=8;
|
||||
const int Npatch = IcosahedralPatches;
|
||||
|
||||
// Put SIMD all in time direction
|
||||
Coordinate latt_size = GridDefaultLatt();
|
||||
Coordinate simd_layout({1,1,vComplex::Nsimd(),1});
|
||||
Coordinate mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::cout << GridLogMessage << " mpi "<<mpi_layout<<std::endl;
|
||||
std::cout << GridLogMessage << " simd "<<simd_layout<<std::endl;
|
||||
std::cout << GridLogMessage << " latt "<<latt_size<<std::endl;
|
||||
GridCartesianCrossIcosahedron EdgeGrid(latt_size,simd_layout,mpi_layout,IcosahedralEdges);
|
||||
std::cout << GridLogMessage << " Created edge grid "<<std::endl;
|
||||
GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
|
||||
|
||||
std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
|
||||
LatticeIcosahedralLorentzColourMatrix Umu(&EdgeGrid);
|
||||
LatticeComplex Phi(&VertexGrid);
|
||||
std::cout << GridLogMessage << " Created two fields "<<std::endl;
|
||||
|
||||
Phi = Zero();
|
||||
Umu = Zero();
|
||||
std::cout << GridLogMessage << " Zeroed two fields "<<std::endl;
|
||||
|
||||
Complex one (1.0);
|
||||
Phi = one;
|
||||
Umu = one;
|
||||
|
||||
std::cout << GridLogMessage << " V = "<<norm2(Phi)<<std::endl;
|
||||
std::cout << GridLogMessage << " Expect "<<latt_size[0]*latt_size[1]*latt_size[2]*10+2*latt_size[2]<<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << " E = "<<norm2(Umu)<<std::endl;
|
||||
std::cout << GridLogMessage << " Expect "<<latt_size[0]*latt_size[1]*latt_size[2]*10*4<<std::endl;
|
||||
|
||||
// std::cout << " Umu "<<Umu<<std::endl;
|
||||
// std::cout << " Phi "<<Phi<<std::endl;
|
||||
LatticePole(Phi,South);
|
||||
// std::cout << " Phi South Pole set\n"<<Phi<<std::endl;
|
||||
|
||||
LatticePole(Phi,North);
|
||||
// std::cout << " Phi North Pole set\n"<<Phi<<std::endl;
|
||||
|
||||
for(int mu=0;mu<VertexGrid._ndimension;mu++){
|
||||
std::cout << " Calling lattice coordinate mu="<<mu<<std::endl;
|
||||
LatticeCoordinate(Phi,mu);
|
||||
// std::cout << " Phi coor mu="<<mu<<"\n"<<Phi<<std::endl;
|
||||
}
|
||||
|
||||
std::cout << "Creating face stencil"<<std::endl;
|
||||
IcosahedralEdgeSupport<IcosahedralGimpl> Support(&VertexGrid,&EdgeGrid);
|
||||
|
||||
std::cout << " Calling Test Geometry "<<std::endl;
|
||||
Support.FaceStencil.TestGeometry();
|
||||
|
||||
|
||||
Umu=one;
|
||||
LatticeComplex plaq1(&EdgeGrid);
|
||||
LatticeComplex plaq2(&EdgeGrid);
|
||||
|
||||
Support.ForwardTriangles(Umu,plaq1,plaq2);
|
||||
std::cout << " plaq1 "<< norm2(plaq1)<<std::endl;
|
||||
std::cout << " plaq2 "<< norm2(plaq2)<<std::endl;
|
||||
|
||||
// Random gauge xform
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
GridParallelRNG vRNG(&EdgeGrid); vRNG.SeedFixedIntegers(seeds);
|
||||
LatticeIcosahedralColourMatrix g(&VertexGrid);
|
||||
|
||||
// SU<Nc>::LieRandomize(vRNG,g);
|
||||
|
||||
LatticeReal gr(&VertexGrid);
|
||||
LatticeComplex gc(&VertexGrid);
|
||||
gaussian(vRNG,gr);
|
||||
Complex ci(0.0,1.0);
|
||||
gc = toComplex(gr);
|
||||
g=one;
|
||||
g = g * exp(ci*gc);
|
||||
|
||||
std::cout << "applying gauge transform"<<std::endl;
|
||||
Support.GaugeTransform(g,Umu);
|
||||
std::cout << "applied gauge transform "<<Umu<<std::endl;
|
||||
|
||||
Support.ForwardTriangles(Umu,plaq1,plaq2);
|
||||
std::cout << " plaq1 "<< norm2(plaq1)<<std::endl;
|
||||
std::cout << " plaq2 "<< norm2(plaq2)<<std::endl;
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user