portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-04-13 23:46:09 +01:00
Code Releases Activity

Covariance test of covariant laplacian appears to pass

Browse Source
This commit is contained in:
Peter Boyle
2025-10-27 19:19:30 -04:00
parent d0ee38d1da
commit d81d00a889
2 changed files with 356 additions and 98 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
Grid/stencil/IcosahedralStencil.h
View File

@@ -69,6 +69,7 @@ public:
protected: protected:
GridBase * _grid; GridBase * _grid;
GridBase * _vertexgrid;
public: public:
GridBase *Grid(void) const { return _grid; } GridBase *Grid(void) const { return _grid; }
@@ -572,11 +573,13 @@ public:
std::cout << GridLogMessage<< " Icosahedral Stencil Geometry Test Complete"<<std::endl; std::cout << GridLogMessage<< " Icosahedral Stencil Geometry Test Complete"<<std::endl;
std::cout << GridLogMessage<< "*************************************"<<std::endl; std::cout << GridLogMessage<< "*************************************"<<std::endl;
} }
IcosahedralStencil(GridBase *grid) // Must be +1 or -1 IcosahedralStencil(GridBase *grid,GridBase *vertexgrid)
{ {
this->_grid = grid; this->_grid = grid;
this->_vertexgrid = vertexgrid;
// Loop over L^2 x T x npatch and the // Loop over L^2 x T x npatch and the
assert(grid->isIcosahedral()); assert(grid->isIcosahedral());
assert(grid->isIcosahedral());
} }
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// VertexInputs = true implies the neighbour has vertex support // VertexInputs = true implies the neighbour has vertex support
@@ -589,14 +592,11 @@ public:
// can apply a vertex supported link double store to edge supported gauge field // can apply a vertex supported link double store to edge supported gauge field
// can apply a vertex supported laplace or dirac operator vertex supported matter field // can apply a vertex supported laplace or dirac operator vertex supported matter field
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void NearestNeighbourStencil(int vertexOutputs) void NearestNeighbourStencil(int vertexInputs,int vertexOutputs)
{ {
GridBase * grid = this->_grid; GridBase * grid = this->_grid; // the edge grid
GridBase * vertexgrid = this->_vertexgrid;
int vertexInputs = grid->isIcosahedralVertex();
int osites = grid->oSites();
uint64_t cart_sites = grid->CartesianOsites(); uint64_t cart_sites = grid->CartesianOsites();
uint64_t Npole_sites = grid->NorthPoleOsites(); uint64_t Npole_sites = grid->NorthPoleOsites();
uint64_t Spole_sites = grid->SouthPoleOsites(); uint64_t Spole_sites = grid->SouthPoleOsites();
@@ -641,6 +641,7 @@ public:
int Patch = Coor[nd-1]; int Patch = Coor[nd-1];
int HemiPatch = Patch%HemiPatches; int HemiPatch = Patch%HemiPatches;
int Hemisphere= Patch/HemiPatches;
int north = Patch/HemiPatches; int north = Patch/HemiPatches;
int south = 1-north; int south = 1-north;
int isPoleY; int isPoleY;
@@ -681,7 +682,11 @@ public:
int YmHemiPatch = YmCoor[nd-1]%HemiPatches; int YmHemiPatch = YmCoor[nd-1]%HemiPatches;
int YmHemisphere = YmCoor[nd-1]/HemiPatches; int YmHemisphere = YmCoor[nd-1]/HemiPatches;
if ( vertexInputs ) { int DmPatch = DmCoor[nd-1];
int DmHemiPatch = DmCoor[nd-1]%HemiPatches;
int DmHemisphere = DmCoor[nd-1]/HemiPatches;
if ( vertexInputs ) {// Neighbour will live on poles and peer point
//////////////////////////////////////////////// ////////////////////////////////////////////////
// XpCoor stencil entry; consider isPole case // XpCoor stencil entry; consider isPole case
//////////////////////////////////////////////// ////////////////////////////////////////////////
@@ -689,7 +694,8 @@ public:
SE._missing_link = false; SE._missing_link = false;
SE._offset = grid->oIndex(XpCoor); SE._offset = grid->oIndex(XpCoor);
if ( isPoleX ) { if ( isPoleX ) {
SE._offset = grid->PoleSiteForOcoor(Coor); SE._offset = vertexgrid->PoleSiteForOcoor(Coor);
// std::cout << site<<" setting X-Pole site "<<SE._offset<<" for coor "<<Coor<<std::endl;
} }
SE._polarisation = IcosahedronPatchY; SE._polarisation = IcosahedronPatchY;
SE._adjoint = false; SE._adjoint = false;
@@ -701,12 +707,13 @@ public:
SE._missing_link = false; SE._missing_link = false;
SE._offset = grid->oIndex(YpCoor); SE._offset = grid->oIndex(YpCoor);
if ( isPoleY ) { if ( isPoleY ) {
SE._offset = grid->PoleSiteForOcoor(Coor); SE._offset = vertexgrid->PoleSiteForOcoor(Coor);
// std::cout << site<<" setting Y-Pole site "<<SE._offset<<" for coor "<<Coor<<std::endl;
} }
SE._polarisation = IcosahedronPatchX; SE._polarisation = IcosahedronPatchX;
SE._adjoint = false; SE._adjoint = false;
acceleratorPut(this->_entries[lexYp],SE); acceleratorPut(this->_entries[lexYp],SE);
} else { } else { // Neighbour will be a forward edge and connection may be more complicated
//////////////////////////////////////////////// ////////////////////////////////////////////////
// XpCoor stencil entry // XpCoor stencil entry
// Store in look up table // Store in look up table
@@ -718,7 +725,7 @@ public:
SE._offset = grid->oIndex(XpCoor); SE._offset = grid->oIndex(XpCoor);
SE._polarisation = IcosahedronPatchY; SE._polarisation = IcosahedronPatchY;
SE._adjoint = false; SE._adjoint = false;
if ( DpHemiPatch != HemiPatch && south ) { if ( DpHemiPatch != HemiPatch && south ) { // These are the sneaky redirect for edge / faces
SE._offset = grid->oIndex(DpCoor); SE._offset = grid->oIndex(DpCoor);
SE._polarisation = IcosahedronPatchX; SE._polarisation = IcosahedronPatchX;
SE._adjoint = true; SE._adjoint = true;
@@ -732,7 +739,7 @@ public:
SE._offset = grid->oIndex(YpCoor); SE._offset = grid->oIndex(YpCoor);
SE._polarisation = IcosahedronPatchX; SE._polarisation = IcosahedronPatchX;
SE._adjoint = false; SE._adjoint = false;
if ( YpHemiPatch != HemiPatch && north ) { if ( YpHemiPatch != HemiPatch && north ) { // These are the sneaky redirect for edge / faces
SE._offset = grid->oIndex(DpCoor); SE._offset = grid->oIndex(DpCoor);
SE._polarisation = IcosahedronPatchY; SE._polarisation = IcosahedronPatchY;
SE._adjoint = true; SE._adjoint = true;
@@ -774,23 +781,34 @@ public:
///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////
// for DmCoor ; never needed for staples, only for vertex diff ops // for DmCoor ; never needed for staples, only for vertex diff ops
// no polarisation rotation // No polarisation rotation.
// But polarisation rotation is needed for double storing.
///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////
SE._offset = grid->oIndex(DmCoor); SE._offset = grid->oIndex(DmCoor);
SE._polarisation = IcosahedronPatchDiagonal; // should ignore SE._polarisation = IcosahedronPatchDiagonal; // default
if ( (DmHemiPatch != HemiPatch) && (DmHemisphere==Hemisphere) && south ) {
SE._polarisation = IcosahedronPatchX; // Basis rotates
}
if ( DmHemiPatch != HemiPatch && (DmHemisphere==Hemisphere) && north ) {
SE._polarisation = IcosahedronPatchY; // Basis rotates
}
SE._missing_link = missingLink; SE._missing_link = missingLink;
acceleratorPut(this->_entries[lexDm],SE); acceleratorPut(this->_entries[lexDm],SE);
} }
if ( vertexOutputs ) { if ( vertexOutputs ) {
int ndm1 = grid->Nd()-1; int ndm1 = grid->Nd()-1;
if ( grid->ownsSouthPole() ) { if ( vertexgrid->ownsSouthPole() ) {
IcosahedralStencilEntry SE; IcosahedralStencilEntry SE;
for(uint64_t site=0;site<cart_sites; site ++) { for(uint64_t site=0;site<cart_sites; site ++) { // loops over volume
Coordinate Coor; Coordinate Coor;
grid->oCoorFromOindex(Coor,site); vertexgrid->oCoorFromOindex(Coor,site);
if( (Coor[0]==L)&&(Coor[1]==0) ) { int north = Coor[ndm1]/HemiPatches;
int64_t pole_site = grid->PoleSiteForOcoor(Coor); int south = 1-north;
int64_t lex = pole_site*np+Coor[ndm1]; if( (Coor[0]==(L-1))&&(Coor[1]==0) &&south ) {
int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor);
int64_t lex = pole_site*np+(Coor[ndm1]%HemiPatches);
// std::cout << "Coor "<<Coor<<" connects to south pole_site "<<pole_site<<std::endl;
SE._offset = site; SE._offset = site;
SE._is_local = true; SE._is_local = true;
SE._polarisation = IcosahedronPatchX; // ignored SE._polarisation = IcosahedronPatchX; // ignored
@@ -804,17 +822,19 @@ public:
} }
} }
} }
if ( grid->ownsNorthPole() ) { if ( vertexgrid->ownsNorthPole() ) {
IcosahedralStencilEntry SE; IcosahedralStencilEntry SE;
for(uint64_t site=0;site<cart_sites; site ++) { for(uint64_t site=0;site<cart_sites; site ++) {
Coordinate Coor; Coordinate Coor;
grid->oCoorFromOindex(Coor,site); vertexgrid->oCoorFromOindex(Coor,site);
if( (Coor[0]==0)&&(Coor[1]==L) ) { int north = Coor[ndm1]/HemiPatches;
int64_t pole_site = grid->PoleSiteForOcoor(Coor); if( (Coor[0]==0)&&(Coor[1]==(L-1))&&north ) {
int64_t lex = pole_site*np+Coor[ndm1]; int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor);
int64_t lex = pole_site*np+(Coor[ndm1]%HemiPatches);
// std::cout << "Coor "<<Coor<<" connects to north pole_site "<<pole_site<<std::endl;
SE._offset = site; SE._offset = site;
SE._is_local = true; SE._is_local = true;
SE._polarisation = IcosahedronPatchX; // ignored SE._polarisation = IcosahedronPatchY; // ignored
SE._adjoint = false; // ignored SE._adjoint = false; // ignored
SE._missing_link = false; SE._missing_link = false;
acceleratorPut(this->_entries[lex],SE); acceleratorPut(this->_entries[lex],SE);
@@ -986,3 +1006,4 @@ public:
*/ */
}; };
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

375
tests/debug/Test_icosahedron.cc
View File

@@ -1,4 +1,4 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@@ -33,30 +33,13 @@ using namespace Grid;
const int MyNd=3; const int MyNd=3;
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> > > ;
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 class IcosahedralGimpl
{ {
public: public:
typedef LatticeIcosahedralLorentzColourMatrix GaugeField; typedef LatticeLorentzColourMatrix GaugeField;
typedef LatticeIcosahedralColourMatrix GaugeLinkField; typedef LatticeColourMatrix GaugeLinkField;
typedef LatticeComplex ComplexField; typedef LatticeDoubledGaugeField DoubledGaugeField;
typedef LatticeComplex ComplexField;
}; };
template< class Gimpl> template< class Gimpl>
@@ -66,6 +49,7 @@ public:
typedef typename Gimpl::GaugeField GaugeField; typedef typename Gimpl::GaugeField GaugeField;
typedef typename Gimpl::GaugeLinkField GaugeLinkField; typedef typename Gimpl::GaugeLinkField GaugeLinkField;
typedef typename Gimpl::ComplexField ComplexField; typedef typename Gimpl::ComplexField ComplexField;
typedef typename Gimpl::DoubledGaugeField DoubledGaugeField;
// //
GridBase *VertexGrid; GridBase *VertexGrid;
GridBase *EdgeGrid; GridBase *EdgeGrid;
@@ -74,14 +58,23 @@ public:
IcosahedralStencil NNee; // edge neighbours with edge domain IcosahedralStencil NNee; // edge neighbours with edge domain
IcosahedralStencil NNev; // vertex neighbours but in edge domain IcosahedralStencil NNev; // vertex neighbours but in edge domain
IcosahedralStencil NNvv; // vertex neighbours with vertex domain IcosahedralStencil NNvv; // vertex neighbours with vertex domain
IcosahedralStencil NNve; // edge neighbours with vertex domain
IcosahedralSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid) IcosahedralSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid)
: FaceStencil (_EdgeGrid), NNee(_EdgeGrid), NNev(_VertexGrid), NNvv(_VertexGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid) : FaceStencil (_EdgeGrid,_VertexGrid),
NNee(_EdgeGrid,_VertexGrid),
NNev(_EdgeGrid,_VertexGrid),
NNve(_EdgeGrid,_VertexGrid),
NNvv(_EdgeGrid,_VertexGrid),
VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
{ {
FaceStencil.FaceStencil(); FaceStencil.FaceStencil();
NNee.NearestNeighbourStencil(false);// Edge nearest neighbour std::cout << "NNee"<<std::endl;
NNev.NearestNeighbourStencil(false);// Edge result, vertex neighbour // true/false is "vertexInput, VertexOutput"
NNvv.NearestNeighbourStencil(true); // vertex result and neighbour NNee.NearestNeighbourStencil(false,false);// edge input + output ; used by face stencil
NNev.NearestNeighbourStencil(true,false); // vertex input, edge ouput ; used by gauge transform
NNvv.NearestNeighbourStencil(true,true); // vertex input + output ; used by Laplacian
NNve.NearestNeighbourStencil(false,true); // edge input, vertex output; used by double store
} }
//////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////
@@ -253,9 +246,6 @@ public:
}); });
} }
/*
* Should be able to use the Vertex based stencil to do the GT, picking forward hops
*/
template<class MatterField> template<class MatterField>
void Laplacian(MatterField &in,MatterField &out) void Laplacian(MatterField &in,MatterField &out)
{ {
@@ -312,29 +302,27 @@ public:
coalescedWrite(out_v[ss](),o); coalescedWrite(out_v[ss](),o);
}); });
} }
void GaugeTransform(GaugeLinkField &gt, GaugeField &Umu) template<class MatterField>
void CovariantLaplacian(MatterField &in,MatterField &out,DoubledGaugeField &Uds)
{ {
autoView(Umu_v,Umu,AcceleratorWrite); autoView(out_v,out,AcceleratorWrite);
autoView(g_v,gt,AcceleratorRead); autoView(in_v,in,AcceleratorRead);
autoView(stencil_v,NNev,AcceleratorRead); autoView(U_v,Uds,AcceleratorRead);
autoView(stencil_v,NNvv,AcceleratorRead);
const int np = NNev._npoints; const int np = NNvv._npoints;
std::cout << GridLogMessage<< "GaugeTransform via STENCIL "<<std::endl;
const int ent_Xp = 0; const int ent_Xp = 0;
const int ent_Yp = 1; const int ent_Yp = 1;
const int ent_Dp = 2; const int ent_Dp = 2;
const int ent_Xm = 3;
accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{ const int ent_Ym = 4;
const int ent_Dm = 5;
// Three forward links from this site
auto Lx = Umu_v(ss)(IcosahedronPatchX); accelerator_for(ss,VertexGrid->oSites(),vComplex::Nsimd(),{
auto Ly = Umu_v(ss)(IcosahedronPatchY);
auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal);
auto SE = stencil_v.GetEntry(ent_Xp,ss); auto SE = stencil_v.GetEntry(ent_Xp,ss);
uint64_t xp_idx = SE->_offset; uint64_t xp_idx = SE->_offset;
@@ -344,26 +332,157 @@ public:
SE = stencil_v.GetEntry(ent_Dp,ss); SE = stencil_v.GetEntry(ent_Dp,ss);
uint64_t dp_idx = SE->_offset; uint64_t dp_idx = SE->_offset;
auto g = g_v(ss)(); SE = stencil_v.GetEntry(ent_Xm,ss);
auto gx = g_v(xp_idx)(); uint64_t xm_idx = SE->_offset;
auto gy = g_v(yp_idx)();
auto gd = g_v(dp_idx)();
auto lx = Umu_v(ss)(IcosahedronPatchX); SE = stencil_v.GetEntry(ent_Ym,ss);
auto ly = Umu_v(ss)(IcosahedronPatchY); uint64_t ym_idx = SE->_offset;
auto ld = Umu_v(ss)(IcosahedronPatchDiagonal);
lx = g*lx*adj(gx); SE = stencil_v.GetEntry(ent_Dm,ss);
ly = g*ly*adj(gy); uint64_t dm_idx = SE->_offset;
ld = g*ld*adj(gd); int missingLink = SE->_missing_link;
coalescedWrite(Umu_v[ss](IcosahedronPatchX),lx); auto i = in_v(ss)();
coalescedWrite(Umu_v[ss](IcosahedronPatchY),ly); auto inxp = in_v(xp_idx)();
coalescedWrite(Umu_v[ss](IcosahedronPatchDiagonal),ld); auto inyp = in_v(yp_idx)();
auto indp = in_v(dp_idx)();
auto inxm = in_v(xm_idx)();
auto inym = in_v(ym_idx)();
auto indm = in_v(dm_idx)();
inxp = U_v(ss)(0)*inxp;
inyp = U_v(ss)(1)*inyp;
indp = U_v(ss)(2)*indp;
inxm = U_v(ss)(3)*inxm;
inym = U_v(ss)(4)*inym;
auto o = i;
if ( missingLink ) {
// std::cout << " CL site "<<ss<<" "<<i<<" "<<inxp<<" "<<inyp<<" "<<indp<<" "<<inxm<<" "<<inym<<std::endl;
o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-i;
} else {
// std::cout << " CL site "<<ss<<" "<<i<<" "<<inxp<<" "<<inyp<<" "<<indp<<" "<<inxm<<" "<<inym<<std::endl;
indm = U_v(ss)(5)*indm;
o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-i;
}
coalescedWrite(out_v[ss](),o);
}); });
} }
void DoubleStore(GaugeField &U,DoubledGaugeField &Uds)
{
assert(U.Grid()==EdgeGrid);
assert(Uds.Grid()==VertexGrid);
autoView(Uds_v,Uds,AcceleratorWrite);
autoView(U_v,U,AcceleratorRead);
autoView(stencil_v,NNvv,AcceleratorRead);
// Vertex result
// Edge valued input
// Might need an extra case (?)
const int np = NNvv._npoints;
const int ent_Xm = 3;
const int ent_Ym = 4;
const int ent_Dm = 5;
accelerator_for(ss,VertexGrid->CartesianOsites(),vComplex::Nsimd(),{
/*
Coordinate sc;
Coordinate xc;
Coordinate yc;
Coordinate dc;
int isPoleY;
int isPoleX;
*/
// EdgeGrid->oCoorFromOindex(sc,ss);
// FaceStencil.GetNbrForPlusDiagonal(EdgeGrid,sc,dc);
// FaceStencil.GetNbrForPlusX(EdgeGrid,sc,xc,isPoleX);
// FaceStencil.GetNbrForPlusY(EdgeGrid,sc,yc,isPoleY);
// Three local links
{
auto Lx = U_v(ss)(IcosahedronPatchX);
auto Ly = U_v(ss)(IcosahedronPatchY);
auto Ld = U_v(ss)(IcosahedronPatchDiagonal);
coalescedWrite(Uds_v[ss](IcosahedronPatchX),Lx);
coalescedWrite(Uds_v[ss](IcosahedronPatchY),Ly);
coalescedWrite(Uds_v[ss](IcosahedronPatchDiagonal),Ld);
}
// Three backwards links
{
auto SE = stencil_v.GetEntry(ent_Xm,ss);
auto pol = SE->_polarisation;
auto s = SE->_offset;
int pol1 = IcosahedronPatchX;
//
// Xm link is given diagonal unless HemiPatch changes in Northern hemisphere
// But here for double storing we need the reverse link so
// return either the Xplus or Diag plus direction
//
if ( pol != IcosahedronPatchDiagonal ) {
pol1 = IcosahedronPatchDiagonal;
}
auto Lx_at_xm = U_v(s)(pol1);
Lx_at_xm = adj(Lx_at_xm);
coalescedWrite(Uds_v[ss](3),Lx_at_xm );
// EdgeGrid->oCoorFromOindex(xc,s);
// std::cout << " Coor "<<sc<<" Xm "<<xc<<" stencil entry "<<s<<" "<<Lx_at_xm<<std::endl;
}
{
auto SE = stencil_v.GetEntry(ent_Ym,ss);
auto pol = SE->_polarisation;
auto s = SE->_offset;
//
// Ym link is given diagonal unless HemiPatch changes in the Southern hemisphere
// But here for double storing we need the reverse link so
// return either the Yplus or Diag plus direction
//
int pol1 = IcosahedronPatchY;
if ( pol != IcosahedronPatchDiagonal ) {
pol1 = IcosahedronPatchDiagonal;
}
auto Ly_at_ym = U_v(s)(pol1);
Ly_at_ym = adj(Ly_at_ym);
coalescedWrite(Uds_v[ss](4),Ly_at_ym );
// EdgeGrid->oCoorFromOindex(yc,s);
// std::cout << " Coor "<<sc<<" Ym "<<yc<<" stencil entry "<<s<<" "<<Ly_at_ym<<std::endl;
}
int missingLink;
{
auto SE = stencil_v.GetEntry(ent_Dm,ss);
auto s = SE->_offset;
// Dm link given the correct value for this use case
auto pol= SE->_polarisation;
missingLink = SE->_missing_link;
if ( ! missingLink ) {
auto Ld_at_dm = U_v(s)(pol);
Ld_at_dm = adj(Ld_at_dm);
coalescedWrite(Uds_v[ss](5),Ld_at_dm );
// EdgeGrid->oCoorFromOindex(dc,s);
// std::cout << " Coor "<<sc<<" Dm "<<dc<<" stencil entry "<<s<<" "<<Ld_at_dm<<std::endl;
}
}
});
auto pole_sites = VertexGrid->oSites() - VertexGrid->CartesianOsites();
auto pole_offset= VertexGrid->CartesianOsites();
accelerator_for(ss,pole_sites,vComplex::Nsimd(),{
for(int p=0;p<HemiPatches;p++){ // Neighbours of each pole
auto SE = stencil_v.GetEntry(p,pole_offset+ss);
auto s = SE->_offset;
auto pol= SE->_polarisation;
auto Link = adj(U_v(s)(pol));
coalescedWrite(Uds_v[pole_offset+ss](p),Link);
// Coordinate pc;
// EdgeGrid->oCoorFromOindex(pc,s);
// std::cout << " Stencil pole neighbour hemi "<<p<<" site "<<s<<" "<<pc<<" Link "<<Link<<std::endl;
}
});
}
/* /*
* This routine is slow and single threaded on CPU
void GaugeTransformCPU(GaugeLinkField &gt, GaugeField &Umu) void GaugeTransformCPU(GaugeLinkField &gt, GaugeField &Umu)
{ {
assert(gt.Grid()==VertexGrid); assert(gt.Grid()==VertexGrid);
@@ -475,9 +594,80 @@ public:
coalescedWrite(Umu_v[site](IcosahedronPatchDiagonal),ld); coalescedWrite(Umu_v[site](IcosahedronPatchDiagonal),ld);
}; };
} }
*/
void GaugeTransform(GaugeLinkField &gt, GaugeField &Umu)
{
autoView(Umu_v,Umu,AcceleratorWrite);
autoView(g_v,gt,AcceleratorRead);
autoView(stencil_v,NNev,AcceleratorRead);
const int np = NNev._npoints;
const int ent_Xp = 0;
const int ent_Yp = 1;
const int ent_Dp = 2;
accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{
// Three forward links from this site
auto Lx = Umu_v(ss)(IcosahedronPatchX);
auto Ly = Umu_v(ss)(IcosahedronPatchY);
auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal);
auto SE = stencil_v.GetEntry(ent_Xp,ss);
uint64_t xp_idx = SE->_offset;
SE = stencil_v.GetEntry(ent_Yp,ss);
uint64_t yp_idx = SE->_offset;
SE = stencil_v.GetEntry(ent_Dp,ss);
uint64_t dp_idx = SE->_offset;
auto g = g_v(ss)();
auto gx = g_v(xp_idx)();
auto gy = g_v(yp_idx)();
auto gd = g_v(dp_idx)();
auto lx = Umu_v(ss)(IcosahedronPatchX);
auto ly = Umu_v(ss)(IcosahedronPatchY);
auto ld = Umu_v(ss)(IcosahedronPatchDiagonal);
lx = g*lx*adj(gx);
ly = g*ly*adj(gy);
ld = g*ld*adj(gd);
/*
Coordinate sc;
Coordinate xc;
Coordinate yc;
Coordinate dc;
int isPoleY;
int isPoleX;
EdgeGrid->oCoorFromOindex(sc,ss);
FaceStencil.GetNbrForPlusDiagonal(EdgeGrid,sc,dc);
FaceStencil.GetNbrForPlusX(EdgeGrid,sc,xc,isPoleX);
FaceStencil.GetNbrForPlusY(EdgeGrid,sc,yc,isPoleY);
if (isPoleX) {
std::cout << " ss "<<ss<<" "<<sc<< " POLE xp_idx "<<xp_idx<<" gxp "<<adj(gx)<<std::endl;
} else {
std::cout << " ss "<<ss<<" "<<sc<< " xp_idx "<<xp_idx<<" "<<xc<<" gxp "<<adj(gx)<<std::endl;
}
if (isPoleY) {
std::cout << " ss "<<ss<<" "<<sc<< " POLE yp_idx "<<yp_idx<<" gyp "<<adj(gy)<<std::endl;
} else {
std::cout << " ss "<<ss<<" "<<sc<< " yp_idx "<<yp_idx<<" "<<yc<<" gyp "<<adj(gy)<<std::endl;
}
std::cout << " ss "<<ss<<" "<<sc<< " dp_idx "<<dp_idx<<" "<<dc<<" gdp "<<adj(gd)<<std::endl;
*/
coalescedWrite(Umu_v[ss](IcosahedronPatchX),lx);
coalescedWrite(Umu_v[ss](IcosahedronPatchY),ly);
coalescedWrite(Umu_v[ss](IcosahedronPatchDiagonal),ld);
});
}
}; };
*/
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
@@ -499,8 +689,8 @@ int main (int argc, char ** argv)
GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices); GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
std::cout << GridLogMessage << " Created vertex grid "<<std::endl; std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
LatticeIcosahedralLorentzColourMatrix Umu(&EdgeGrid); LatticeLorentzColourMatrix Umu(&EdgeGrid);
LatticeIcosahedralLorentzColourMatrix Umuck(&EdgeGrid); LatticeLorentzColourMatrix Umuck(&EdgeGrid);
LatticeComplex Phi(&VertexGrid); LatticeComplex Phi(&VertexGrid);
std::cout << GridLogMessage << " Created two fields "<<std::endl; std::cout << GridLogMessage << " Created two fields "<<std::endl;
@@ -557,17 +747,28 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4}); std::vector<int> seeds({1,2,3,4});
GridParallelRNG vRNG(&EdgeGrid); vRNG.SeedFixedIntegers(seeds); GridParallelRNG vRNG(&EdgeGrid); vRNG.SeedFixedIntegers(seeds);
// SU<Nc>::LieRandomize(vRNG,g); LatticeColourMatrix g(&VertexGrid);
LatticeIcosahedralColourMatrix g(&VertexGrid);
LatticeReal gr(&VertexGrid); LatticeReal gr(&VertexGrid);
LatticeComplex gc(&VertexGrid); LatticeComplex gc(&VertexGrid);
gr = 1.0; // gr = Zero();
gaussian(vRNG,gr); gaussian(vRNG,gr);
Complex ci(0.0,1.0); Complex ci(0.0,1.0);
gc = toComplex(gr); gc = toComplex(gr);
g=one; g=one;
g = g * exp(ci*gc); g = g * exp(ci*gc);
/*
Complex m1(1.0,0.0);
typename LatticeComplex::scalar_object sobj; sobj = m1;
int nd = EdgeGrid.Nd();
Coordinate coor(nd,0);
coor[0]=1;
coor[1]=latt_size[1]-1;
coor[2]=2;
coor[3]=5;
// pokeSite(m1,gc,coor);
*/
// std::cout << "GT is "<<g<<std::endl;
std::cout << GridLogMessage << "****************************************"<<std::endl; std::cout << GridLogMessage << "****************************************"<<std::endl;
std::cout << GridLogMessage << " Check plaquette is gauge invariant "<<std::endl; std::cout << GridLogMessage << " Check plaquette is gauge invariant "<<std::endl;
std::cout << GridLogMessage << "****************************************"<<std::endl; std::cout << GridLogMessage << "****************************************"<<std::endl;
@@ -636,10 +837,46 @@ int main (int argc, char ** argv)
// std::cout << " YXD\n " << closure(linkY * stapleXD) <<std::endl; // std::cout << " YXD\n " << closure(linkY * stapleXD) <<std::endl;
std::cout << GridLogMessage<< "Calling Laplacian" <<std::endl; std::cout << GridLogMessage<< "Calling Laplacian" <<std::endl;
LatticeComplex in(&VertexGrid); LatticeColourVector in(&VertexGrid);
LatticeComplex out(&VertexGrid); LatticeColourVector out(&VertexGrid);
LatticeColourVector gout(&VertexGrid);
LatticeColourVector gin(&VertexGrid);
gaussian(vRNG,in); gaussian(vRNG,in);
Support.Laplacian(in,out); Support.Laplacian(in,out);
std::cout << GridLogMessage<< "Calling double storing gauge field" <<std::endl;
LatticeDoubledGaugeField Uds(&VertexGrid);
Support.DoubleStore(Umu,Uds);
// std::cout << "Uds is"<<Uds<<std::endl;
Support.CovariantLaplacian(in,out,Uds);
std::cout << GridLogMessage<< "Applied covariant laplacian !" <<std::endl;
/*
* CovariantLaplacian testing -- check the laplacian is gauge invariant
*
* D[U_gt](gF) = g D[U] g^dag g F = g D[U] F
*/
gout = g*out;
gin = g*in;
Support.GaugeTransform(g,Umu);
Support.DoubleStore(Umu,Uds);
Support.CovariantLaplacian(gin,out,Uds);
std::cout << GridLogMessage<< "Applied gauge transformed covariant laplacian to transformed vector !" <<std::endl;
std::cout << "Testing D[U_gt](gF) = g D[U] F : defect is "<<norm2(out-gout)<<std::endl;
// std::cout << "Support.CovariantLaplacian(gin,out,Uds) "<<out<<std::endl;
// std::cout << "g * out "<<gout<<std::endl;
// out = out - gout;
// std::cout << "diff "<<out<<std::endl;
Grid_finalize(); Grid_finalize();
} }
/**********************************************************************************************
* This routine is slow and single threaded on CPU
* Preserved as potentially useful in future, but only in comments
*********************************************************************************************
*/