portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Merge pull request #411 from giltirn/patch/dirichlet-fixes

Browse Source 
Various fixes / changes
This commit is contained in:
Peter Boyle 2022-10-04 17:22:01 -04:00 committed by GitHub
commit eec0c9eb7d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 248 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Grid/lattice/Lattice_reduction.h
View File

@ -91,10 +91,7 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites)
for(int i=0;i<nthread;i++){
ssum = ssum+sumarray[i];
}
typedef typename vobj::scalar_object ssobj;
ssobj ret = ssum;
return ret;
return ssum;
}
/*
Threaded max, don't use for now

43
Grid/qcd/hmc/HMC.h
View File

@ -53,6 +53,7 @@ struct HMCparameters: Serializable {
Integer, Trajectories, /* @brief Number of sweeps in this run */
bool, MetropolisTest,
Integer, NoMetropolisUntil,
bool, PerformRandomShift, /* @brief Randomly shift the gauge configuration at the start of a trajectory */
std::string, StartingType,
IntegratorParameters, MD)
@ -63,6 +64,7 @@ struct HMCparameters: Serializable {
StartTrajectory = 0;
Trajectories = 10;
StartingType = "HotStart";
PerformRandomShift = true;
/////////////////////////////////
}
@ -83,6 +85,7 @@ struct HMCparameters: Serializable {
std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
std::cout << GridLogMessage << "[HMC parameters] Doing random shift : " << std::boolalpha << PerformRandomShift << "\n";
std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
MD.print_parameters();
}
@ -95,6 +98,7 @@ private:
const HMCparameters Params;
typedef typename IntegratorType::Field Field;
typedef typename IntegratorType::FieldImplementation FieldImplementation;
typedef std::vector< HmcObservable<Field> * > 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;d<Grid->Nd();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<typename FieldImplementation::GaugeLinkField> Umu(Grid->Nd(), U.Grid());
for(int mu=0;mu<Grid->Nd();mu++) Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
RealD rn_uniform; random(sRNG, rn_uniform);
for(int d=0;d<Grid->Nd();d++) {
int shift = (int) (rn_uniform*L);
int L = Grid->GlobalDimensions()[d];
std::cout << shift;
if(d<Grid->Nd()-1) std::cout <<",";
else std::cout <<"]\n";
RealD rn_uniform; random(sRNG, rn_uniform);
int shift = (int) (rn_uniform*L);
std::cout << shift;
if(d<Grid->Nd()-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;mu<Grid->Nd();mu++) PokeIndex<LorentzIndex>(U,Umu[mu],mu);
std::cout << GridLogMessage << "--------------------------------------------------\n";
}
std::cout << GridLogMessage << "--------------------------------------------------\n";
TheIntegrator.reset_timer();

4
Grid/qcd/hmc/integrators/Integrator.h
View File

@ -63,10 +63,10 @@ public:
};
/*! @brief Class for Molecular Dynamics management */
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy>
class Integrator {
protected:
typedef FieldImplementation_ FieldImplementation;
typedef typename FieldImplementation::Field MomentaField; //for readability
typedef typename FieldImplementation::Field Field;

15
Grid/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid);
* P 1/2 P 1/2
*/
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class LeapFrog : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>
template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class LeapFrog : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{
public:
typedef FieldImplementation_ FieldImplementation;
typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy> Algorithm;
INHERIT_FIELD_TYPES(FieldImplementation);
@ -135,13 +136,14 @@ public:
}
};
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class MinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>
template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class MinimumNorm2 : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{
private:
const RealD lambda = 0.1931833275037836;
public:
typedef FieldImplementation_ FieldImplementation;
INHERIT_FIELD_TYPES(FieldImplementation);
MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
@ -192,8 +194,8 @@ public:
}
};
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>
template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class ForceGradient : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{
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.

30
Grid/qcd/utils/CovariantCshift.h
View File

@ -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<class gauge> Lattice<gauge>
CshiftLink(const Lattice<gauge> &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<iScalar<vInteger>> coor(grid);
LatticeCoordinate(coor, mu);
Lattice<gauge> 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;
}
}

4
Grid/qcd/utils/GaugeFix.h
View File

@ -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();

2
Grid/simd/Grid_sse4.h
View File

@ -35,7 +35,7 @@ Author: neo <cossu@post.kek.jp>
*/
// Time-stamp: <2015-06-16 23:27:54 neo>
//----------------------------------------------------------------------
#include <immintrin.h>
#include <pmmintrin.h>
NAMESPACE_BEGIN(Grid);

183
tests/Test_gfield_shift.cc Normal file
View File

@ -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 <ckelly@bnl.gov>
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
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;mu<Grid->Nd();mu++){
Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = Gimpl::CshiftLink(Umu,dir,shift);
PokeIndex<LorentzIndex>(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<int> 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<int>({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<ColourIndex>(link_field, lex, 0,0); //place on 0,0 element of link
for(int mu=0;mu<Nd;mu++){
link_field_2 = link_field + mu*stride; //add in lex-mapping of mu
link_field_2 = ComplexD(0,1) * link_field_2; //make imaginary
PokeIndex<LorentzIndex>(U, link_field_2, mu);
}
}
std::stringstream ss;
ss<<"error";
for(int d=0;d<Fine._ndimension;d++){
ss<<"."<<Fine._processor_coor[d];
}
ss<<"_wr_"<<Fine._processor;
std::string fname(ss.str());
std::ofstream ferr(fname);
Integer vol4d = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
bool fail = false;
typename SiteGaugeField::scalar_object um;
TComplex cm;
for(int dir=0;dir<4;dir++){
for(int shift=-latt_size[dir]+1;shift<latt_size[dir];shift++){
if ( Fine.IsBoss() )
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction "<<dir
<< " dir is conj ? " << conj_dirs[dir] << std::endl;
ShiftU = CshiftGaugeField(U,dir,shift);
Coordinate coor(4);
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
for(coor[0]=0;coor[0]<latt_size[0];coor[0]++){
peekSite(um,ShiftU,coor);
Coordinate scoor(coor);
scoor[dir] = (scoor[dir]+shift + latt_size[dir])%latt_size[dir];
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
for(int mu = 0 ; mu < 4; mu++){
Integer slex_mu = slex + vol4d*mu;
Complex scm(0,slex_mu); //imaginary
if(
( shift > 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<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
ferr<<"Got mu "<< cm_mu << " site " <<index<<" : " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Integer scm_mu = index / vol4d;
index = index % vol4d;
Lexicographic::CoorFromIndex(peer,index,latt_size);
ferr<<"Expect mu " << scm_mu << " site " <<index<<": " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
fail = true;
}
}
}}}}
}
}
if(fail) std::cout << "Test FAILED : see " << fname << " for more details" << std::endl;
else std::cout << "Test Passed" << std::endl;
Grid_finalize();
}