1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-08 15:49:32 +00:00

Compare commits

...

16 Commits

Author SHA1 Message Date
Peter Boyle
d5c0d54f89 Gauge staples for temporal direction added (ico-T staples and T-ico staples).
Passes gauge covariance test, requiring the link x its staples = 1 on a random gauge transform.

Still to do:

temporal link double store
temporal laplacian + covariant laplacian terms

Make the "rotate" -> coalescedReadRotate as presently it is hardwired CPU domain in a couple of places in Test_icosahedron
2025-11-06 13:54:39 -05:00
Peter Boyle
d3ca16c76d Updated 2025-10-27 21:09:02 -04:00
Peter Boyle
d81d00a889 Covariance test of covariant laplacian appears to pass 2025-10-27 19:19:30 -04:00
Peter Boyle
d0ee38d1da Clean up 2025-10-22 21:44:51 -04:00
Peter Boyle
da8dc3da0d More compact 2025-10-22 21:37:40 -04:00
Peter Boyle
21514d8487 Added a free laplacian 2025-10-22 21:31:53 -04:00
Peter Boyle
77b2e9fb61 Name changes 2025-10-22 16:46:15 -04:00
Peter Boyle
a71ba05bd7 Implemented gauge transform via stencil.
Now have ability to do Vertex AND Edge grids
Should now have no barriers to
a) Double Storing links for fermion operators / laplacian
b) Laplace or Wilson operators
2025-10-22 16:27:06 -04:00
Peter Boyle
1e95e64035 Staples work in icoso-plane 2025-10-21 23:27:27 -04:00
Peter Boyle
defcac92ab Somewhat better wrapped support for Icosahedral 2025-10-20 18:09:21 -04:00
Peter Boyle
4869378f1e Now computed some plaquettes.
First cut at stencil
2025-10-20 11:15:17 -04:00
Peter Boyle
c7b74db317 Default dimensions fixed 2025-10-09 14:57:22 -04:00
Peter Boyle
0ce201efbe IcosahedralVerted() checks 2025-10-09 13:35:16 -04:00
Peter Boyle
6d8a3d8bb2 Config 2025-10-09 13:30:16 -04:00
Peter Boyle
7dfd207ebb Need to protect pole operatoins to only take place on IcosahedralVertices mesh 2025-10-08 15:18:31 -04:00
Peter Boyle
3a65a096f2 Nd verbose 2025-10-07 18:49:00 -04:00
14 changed files with 3199 additions and 69 deletions

View File

@@ -41,8 +41,17 @@ enum NorthSouth {
North = 1,
South = 0
};
enum IcoshedralDirections {
IcosahedronPatchX = 0,
IcosahedronPatchY = 1,
IcosahedronPatchDiagonal=2,
NumIcosahedralPolarizations
};
const int num_icosahedron_tiles = 10;
const int IcosahedralPatches = 10;
const int HemiPatches=IcosahedralPatches/2;
const int NorthernHemisphere = HemiPatches;
const int SouthernHemisphere = 0;
class GridCartesianCrossIcosahedron: public GridCartesian {
@@ -86,9 +95,12 @@ public:
assert(simd_layout[0]==1); // Force simd into perpendicular dimensions
assert(simd_layout[1]==1); // to avoid pole storage complexity interacting with SIMD.
assert(dimensions[_ndimension-1]==num_icosahedron_tiles);
assert(dimensions[_ndimension-1]==IcosahedralPatches);
assert(processor_grid[_ndimension-1]<=2); // Keeps the patches that need a pole on the same node
// Save a copy of the basic cartesian initialisation volume
cartesianOsites = this->_osites;
// allocate the pole storage if we are seeking vertex domain data
if ( meshType == IcosahedralVertices ) {
InitPoles();
@@ -106,14 +118,46 @@ public:
int southPoleOsite;
int northPoleOsites;
int southPoleOsites;
int cartesianOsites;
virtual int Icosahedral(void) override { return 1;}
virtual int ownsNorthPole(void) const override { return hasNorthPole; };
virtual int isIcosahedral(void) override { return 1;}
virtual int isIcosahedralVertex(void) override { return meshType==IcosahedralVertices;}
virtual int isIcosahedralEdge (void) override { return meshType==IcosahedralEdges;}
virtual int NorthPoleOsite(void) const override { return northPoleOsite; };
virtual int NorthPoleOsites(void) const override { return northPoleOsites; };
virtual int ownsSouthPole(void) const override { return hasSouthPole; };
virtual int SouthPoleOsite(void) const override { return southPoleOsite; };
virtual int SouthPoleOsites(void) const override { return southPoleOsites; };
virtual int ownsNorthPole(void) const override { return hasNorthPole; };
virtual int ownsSouthPole(void) const override { return hasSouthPole; };
virtual int CartesianOsites(void) const override { return cartesianOsites; };
virtual int64_t PoleIdxForOcoor(Coordinate &Coor) override
{
// Work out the pole_osite. Pick the higher dims
Coordinate rdims;
Coordinate ocoor;
int64_t pole_idx;
int Ndm1 = this->Nd()-1;
for(int d=2;d<Ndm1;d++){
int dd=d-2;
rdims.push_back(this->_rdimensions[d]);
ocoor.push_back(Coor[d]%this->_rdimensions[d]);
}
Lexicographic::IndexFromCoor(ocoor,pole_idx,rdims);
return pole_idx;
}
virtual int64_t PoleSiteForOcoor(Coordinate &Coor) override
{
int Ndm1 = this->Nd()-1;
int64_t pole_idx = this->PoleIdxForOcoor(Coor);
int64_t pole_osite;
if ( Coor[Ndm1] >= HemiPatches ) {
pole_osite = pole_idx + this->NorthPoleOsite();
} else {
pole_osite = pole_idx + this->SouthPoleOsite();
}
return pole_osite;
}
void InitPoles(void)
{
@@ -164,10 +208,10 @@ public:
* Hence all 5 patches associated with the pole must have the
* appropriate "corner" of the patch L^2 located on the SAME rank.
*/
if( (pcoor[xdim]==pgrid[xdim]-1) && (pcoor[ydim]==0) && (pcoor[Ndm1]==0) ){
hasSouthPole =1;
southPoleOsite=this->_osites;
southPoleOsite=this->_osites;
southPoleOsites=OrthogSize;
this->_osites += OrthogSize;
} else {
@@ -185,11 +229,11 @@ public:
northPoleOsites=0;
northPoleOsite=0;
}
std::cout << "Icosahedral vertex field volume " << this->_osites<<std::endl;
std::cout << "Icosahedral south pole offset " << this->southPoleOsite<<std::endl;
std::cout << "Icosahedral north pole offset " << this->northPoleOsite<<std::endl;
std::cout << "Icosahedral south pole size " << this->southPoleOsites<<std::endl;
std::cout << "Icosahedral north pole size " << this->northPoleOsites<<std::endl;
std::cout << GridLogDebug<<"Icosahedral vertex field volume " << this->_osites<<std::endl;
std::cout << GridLogDebug<<"Icosahedral south pole offset " << this->southPoleOsite<<std::endl;
std::cout << GridLogDebug<<"Icosahedral north pole offset " << this->northPoleOsite<<std::endl;
std::cout << GridLogDebug<<"Icosahedral south pole size " << this->southPoleOsites<<std::endl;
std::cout << GridLogDebug<<"Icosahedral north pole size " << this->northPoleOsites<<std::endl;
};
};

View File

@@ -87,13 +87,18 @@ public:
public:
// Icosahedral decisions
virtual int Icosahedral(void) { return 0;}
virtual int isIcosahedral(void) { return 0;}
virtual int isIcosahedralVertex(void) { return 0;}
virtual int isIcosahedralEdge (void) { return 0;}
virtual int ownsNorthPole(void) const { return 0; };
virtual int ownsSouthPole(void) const { return 0; };
virtual int NorthPoleOsite(void) const { return 0; };
virtual int SouthPoleOsite(void) const { return 0; };
virtual int NorthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
virtual int SouthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
virtual int CartesianOsites(void) const { return this->oSites(); };
virtual int64_t PoleIdxForOcoor(Coordinate &Coor) { return 0;};
virtual int64_t PoleSiteForOcoor(Coordinate &Coor){ return 0;}
////////////////////////////////////////////////////////////////
// Checkerboarding interface is virtual and overridden by

View File

@@ -34,7 +34,7 @@ NAMESPACE_BEGIN(Grid);
const int Cshift_verbose=0;
template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
{
assert(!rhs.Grid()->Icosahedral());
assert(!rhs.Grid()->isIcosahedral());
typedef typename vobj::vector_type vector_type;
typedef typename vobj::scalar_type scalar_type;

View File

@@ -30,7 +30,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
NAMESPACE_BEGIN(Grid);
template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
{
assert(!rhs.Grid()->Icosahedral());
assert(!rhs.Grid()->isIcosahedral());
Lattice<vobj> ret(rhs.Grid());
ret.Checkerboard() = rhs.Grid()->CheckerBoardDestination(rhs.Checkerboard(),shift,dimension);
Cshift_local(ret,rhs,dimension,shift);

View File

@@ -391,7 +391,7 @@ template<class vobj> std::ostream& operator<< (std::ostream& stream, const Latti
stream<<"]\t";
stream<<ss<<std::endl;
}
if ( o.Grid()->Icosahedral() ) {
if ( o.Grid()->isIcosahedralVertex() ) {
uint64_t psites=1;
Coordinate perpdims;
for(int d=2;d<o.Grid()->_ndimension-1;d++){

View File

@@ -40,7 +40,7 @@ template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
int Nsimd = grid->iSites();
int cartesian_vol = grid->oSites();
if ( grid->Icosahedral() ) {
if ( grid->isIcosahedral() ) {
cartesian_vol = cartesian_vol - grid->NorthPoleOsites()-grid->SouthPoleOsites();
}
{
@@ -58,7 +58,7 @@ template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
});
}
if (grid->Icosahedral()) {
if (grid->isIcosahedralVertex()) {
uint64_t psites=1;
Coordinate perpdims;
typename iobj::scalar_object ss;
@@ -95,7 +95,9 @@ template<class iobj> inline void LatticePole(Lattice<iobj> &l,NorthSouth pole)
l=Zero();
if (grid->Icosahedral()) {
assert(grid->isIcosahedralVertex());
if (grid->isIcosahedralVertex()) {
uint64_t psites=1;
Coordinate perpdims;
sobj ss;

View File

@@ -159,7 +159,8 @@ void peekPole(sobj &s,const Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
GridBase *grid=l.Grid();
assert(grid->Icosahedral());
assert(grid->isIcosahedral());
assert(grid->isIcosahedralVertex());
int Nsimd = grid->Nsimd();
@@ -219,7 +220,8 @@ void pokePole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
{
GridBase *grid=l.Grid();
assert(grid->Icosahedral());
assert(grid->isIcosahedral());
assert(grid->isIcosahedralVertex());
grid->Broadcast(grid->BossRank(),s);
@@ -279,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
//////////////////////////////////////////////////////////

View File

@@ -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;
}

View File

@@ -252,7 +252,7 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexF &o){
inline std::ostream& operator<< (std::ostream& stream, const vComplexD &o){
int nn=vComplexD::Nsimd();
std::vector<ComplexD,alignedAllocator<ComplexD> > buf(nn);
std::vector<ComplexD> buf(nn);
vstore(o,&buf[0]);
stream<<"<";
for(int i=0;i<nn;i++){
@@ -272,7 +272,7 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexD2 &o){
inline std::ostream& operator<< (std::ostream& stream, const vRealF &o){
int nn=vRealF::Nsimd();
std::vector<RealF,alignedAllocator<RealF> > buf(nn);
std::vector<RealF> buf(nn);
vstore(o,&buf[0]);
stream<<"<";
for(int i=0;i<nn;i++){

File diff suppressed because it is too large Load Diff

View File

@@ -187,8 +187,9 @@ void GridParseLayout(char **argv,int argc,
Coordinate &latt_c,
Coordinate &mpi_c)
{
auto mpi =std::vector<int>(1,Nd);
auto latt=std::vector<int>(8,Nd);
auto mpi =std::vector<int>(Nd,1);
auto latt=std::vector<int>(Nd,8);
GridThread::SetMaxThreads();
@@ -228,6 +229,9 @@ void GridParseLayout(char **argv,int argc,
}
// Copy back into coordinate format
int nd = mpi.size();
// std::cout << "mpi.size() "<<nd<<std::endl;
// std::cout << "latt.size() "<<latt.size()<<std::endl;
// std::cout << "Nd "<<Nd<<std::endl;
assert(latt.size()==nd);
latt_c.resize(nd);
mpi_c.resize(nd);

View File

@@ -226,7 +226,7 @@ case ${ac_Nd} in
4)
AC_DEFINE([Config_Nd],[4],[Gauge field dimension Nd]);;
*)
AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);;
AC_MSG_ERROR(["Unsupport dimension Nd = ${ac_Nd}"]);;
esac
############### Symplectic group
@@ -835,6 +835,7 @@ os (target) : $target_os
compiler vendor : ${ax_cv_cxx_compiler_vendor}
compiler version : ${ax_cv_gxx_version}
----- BUILD OPTIONS -----------------------------------
Nd : ${ac_Nd}
Nc : ${ac_Nc}
SIMD : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
Threading : ${ac_openmp}

View File

@@ -0,0 +1,12 @@
MPICXX=mpicxx CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=clang++ ../../configure \
--enable-simd=GEN \
--enable-Nc=1 \
--enable-debug \
--enable-unified=yes \
--prefix $HOME/QCD/GridInstall \
--with-lime=/Users/peterboyle/QCD/SciDAC/install/ \
--with-openssl=$BREW \
--disable-fermion-reps \
--disable-gparity \
--enable-debug

File diff suppressed because it is too large Load Diff