diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h index bcd09c04..af6f8c00 100644 --- a/Grid/lattice/Lattice_reduction.h +++ b/Grid/lattice/Lattice_reduction.h @@ -91,10 +91,7 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites) for(int i=0;i * > ObsListType; //pass these from the resource manager @@ -138,26 +142,37 @@ private: GridBase *Grid = U.Grid(); - ////////////////////////////////////////////////////////////////////////////////////////////////////// - // Mainly for DDHMC perform a random translation of U modulo volume - ////////////////////////////////////////////////////////////////////////////////////////////////////// - std::cout << GridLogMessage << "--------------------------------------------------\n"; - std::cout << GridLogMessage << "Random shifting gauge field by ["; - for(int d=0;dNd();d++) { + if(Params.PerformRandomShift){ + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // Mainly for DDHMC perform a random translation of U modulo volume + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Random shifting gauge field by ["; - int L = Grid->GlobalDimensions()[d]; + std::vector Umu(Grid->Nd(), U.Grid()); + for(int mu=0;muNd();mu++) Umu[mu] = PeekIndex(U, mu); - RealD rn_uniform; random(sRNG, rn_uniform); + for(int d=0;dNd();d++) { - int shift = (int) (rn_uniform*L); + int L = Grid->GlobalDimensions()[d]; - std::cout << shift; - if(dNd()-1) std::cout <<","; - else std::cout <<"]\n"; + RealD rn_uniform; random(sRNG, rn_uniform); + + int shift = (int) (rn_uniform*L); + + std::cout << shift; + if(dNd()-1) std::cout <<","; + else std::cout <<"]\n"; - U = Cshift(U,d,shift); + //shift all fields together in a way that respects the gauge BCs + for(int mu=0; mu < Grid->Nd(); mu++) + Umu[mu] = FieldImplementation::CshiftLink(Umu[mu],d,shift); + } + + for(int mu=0;muNd();mu++) PokeIndex(U,Umu[mu],mu); + + std::cout << GridLogMessage << "--------------------------------------------------\n"; } - std::cout << GridLogMessage << "--------------------------------------------------\n"; TheIntegrator.reset_timer(); diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 7f676ce7..f36e1f64 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -63,10 +63,10 @@ public: }; /*! @brief Class for Molecular Dynamics management */ -template +template class Integrator { protected: - + typedef FieldImplementation_ FieldImplementation; typedef typename FieldImplementation::Field MomentaField; //for readability typedef typename FieldImplementation::Field Field; diff --git a/Grid/qcd/hmc/integrators/Integrator_algorithm.h b/Grid/qcd/hmc/integrators/Integrator_algorithm.h index b05c4ea8..d8b34571 100644 --- a/Grid/qcd/hmc/integrators/Integrator_algorithm.h +++ b/Grid/qcd/hmc/integrators/Integrator_algorithm.h @@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid); * P 1/2 P 1/2 */ -template > -class LeapFrog : public Integrator +template > +class LeapFrog : public Integrator { public: + typedef FieldImplementation_ FieldImplementation; typedef LeapFrog Algorithm; INHERIT_FIELD_TYPES(FieldImplementation); @@ -135,13 +136,14 @@ public: } }; -template > -class MinimumNorm2 : public Integrator +template > +class MinimumNorm2 : public Integrator { private: const RealD lambda = 0.1931833275037836; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) @@ -192,8 +194,8 @@ public: } }; -template > -class ForceGradient : public Integrator +template > +class ForceGradient : public Integrator { private: const RealD lambda = 1.0 / 6.0; @@ -202,6 +204,7 @@ private: const RealD theta = 0.0; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); // Looks like dH scales as dt^4. tested wilson/wilson 2 level. diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h index 79cf8e0f..23984145 100644 --- a/Grid/qcd/utils/CovariantCshift.h +++ b/Grid/qcd/utils/CovariantCshift.h @@ -227,26 +227,38 @@ namespace ConjugateBC { //shift = -1 //Out(x) = U_\mu(x-mu) | x_\mu != 0 // = U*_\mu(L-1) | x_\mu == 0 + //shift = 2 + //Out(x) = U_\mu(x+2\hat\mu) | x_\mu < L-2 + // = U*_\mu(1) | x_\mu == L-1 + // = U*_\mu(0) | x_\mu == L-2 + //shift = -2 + //Out(x) = U_\mu(x-2mu) | x_\mu > 1 + // = U*_\mu(L-2) | x_\mu == 0 + // = U*_\mu(L-1) | x_\mu == 1 + //etc template Lattice CshiftLink(const Lattice &Link, int mu, int shift) { GridBase *grid = Link.Grid(); - int Lmu = grid->GlobalDimensions()[mu] - 1; + int Lmu = grid->GlobalDimensions()[mu]; + assert(abs(shift) < Lmu && "Invalid shift value"); Lattice> coor(grid); LatticeCoordinate(coor, mu); Lattice tmp(grid); - if(shift == 1){ - tmp = Cshift(Link, mu, 1); - tmp = where(coor == Lmu, conjugate(tmp), tmp); + if(shift > 0){ + tmp = Cshift(Link, mu, shift); + tmp = where(coor >= Lmu-shift, conjugate(tmp), tmp); return tmp; - }else if(shift == -1){ + }else if(shift < 0){ tmp = Link; - tmp = where(coor == Lmu, conjugate(tmp), tmp); - return Cshift(tmp, mu, -1); - }else assert(0 && "Invalid shift value"); - return tmp; //shuts up the compiler fussing about the return type + tmp = where(coor >= Lmu+shift, conjugate(tmp), tmp); + return Cshift(tmp, mu, shift); + } + + //shift == 0 + return Link; } } diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index d9d03c54..fc723fe3 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -72,12 +72,12 @@ public: //Fix the gauge field Umu //0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { GridBase *grid = Umu.Grid(); GaugeMat xform(grid); SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge); } - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { //Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform GridBase *grid = Umu.Grid(); diff --git a/Grid/simd/Grid_sse4.h b/Grid/simd/Grid_sse4.h index eb76427e..2b2a80fd 100644 --- a/Grid/simd/Grid_sse4.h +++ b/Grid/simd/Grid_sse4.h @@ -35,7 +35,7 @@ Author: neo */ // Time-stamp: <2015-06-16 23:27:54 neo> //---------------------------------------------------------------------- - +#include #include NAMESPACE_BEGIN(Grid); diff --git a/tests/Test_gfield_shift.cc b/tests/Test_gfield_shift.cc new file mode 100644 index 00000000..cf450d3c --- /dev/null +++ b/tests/Test_gfield_shift.cc @@ -0,0 +1,183 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_gfield_shift.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Azusa Yamaguchi +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +//Test the shifting of the gauge field that respects the boundary conditions +#include + +using namespace Grid; + ; + +typedef ConjugateGimplR Gimpl; //can choose periodic / charge conjugate directions at wil +typedef Gimpl::GaugeField GaugeField; +typedef Gimpl::GaugeLinkField GaugeLinkField; +typedef Gimpl::SiteGaugeField SiteGaugeField; +typedef Gimpl::SiteGaugeLink SiteGaugeLink; + +GaugeField CshiftGaugeField(const GaugeField &U, const int dir, const int shift){ + GridBase *Grid = U.Grid(); + + GaugeField out(Grid); + GaugeLinkField Umu(Grid); + for(int mu=0;muNd();mu++){ + Umu = PeekIndex(U, mu); + Umu = Gimpl::CshiftLink(Umu,dir,shift); + PokeIndex(out,Umu,mu); + } + return out; +} + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + + std::vector conj_dirs = {1,1,0,0}; + Gimpl::setDirections(conj_dirs); + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + + GridParallelRNG FineRNG(&Fine); FineRNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + + GaugeField U(&Fine); + GaugeField ShiftU(&Fine); + + GaugeLinkField link_field(&Fine), link_field_2(&Fine); + + //Like Test_cshift we put the lex coordinate index on each site but make it imaginary + //so we can tell when it was complex conjugated + LatticeComplex lex(&Fine); + lex=Zero(); + U = Zero(); + { + LatticeComplex coor(&Fine); + Integer stride =1; + for(int d=0;d<4;d++){ + LatticeCoordinate(coor,d); + lex = lex + coor*stride; + stride=stride*latt_size[d]; + } + PokeIndex(link_field, lex, 0,0); //place on 0,0 element of link + + for(int mu=0;mu(U, link_field_2, mu); + } + } + + std::stringstream ss; + ss<<"error"; + for(int d=0;d 0 && coor[dir] >= latt_size[dir]-shift && conj_dirs[dir] ) + || + ( shift < 0 && coor[dir] <= -shift-1 && conj_dirs[dir] ) + ) + scm = conjugate(scm); //CC if pulled over boundary + + cm = um(mu)()(0,0); + + RealD nrm = abs(scm-cm()()()); + //std::cout << cm << " " << scm << std::endl; + + Coordinate peer(4); + Complex tmp =cm; + Integer index=real(tmp); + + Integer cm_mu = index / vol4d; + index = index % vol4d; + Lexicographic::CoorFromIndex(peer,index,latt_size); + + if (nrm > 0){ + ferr<<"FAIL mu " << mu << " shift "<< shift<<" in dir "<< dir<<" ["<