portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Added a CshiftLink function to the GaugeImplementations and boundary condition classes that offers a boundary aware C-shift

Browse Source 
Modified gauge fixing code to use CshiftLink internally such that the steepest descent algorithm is universal
Modified gauge transformation code to use CshiftLink for a universal definition
Improved comprehensibility of Test_fft_gfix and generalized to use either periodic or charge conjugation BCs based on cmdline option
Added cmdline options to Test_fft_gfix to tune alpha and optionally disable the Fourier acceleration tests
This commit is contained in:
Christopher Kelly 2021-07-12 17:13:40 -04:00
parent 75a1f85162
commit 5b36a8af54
5 changed files with 259 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
Grid/qcd/action/gauge/GaugeImplementations.h
View File

@ -69,6 +69,11 @@ public:
return PeriodicBC::ShiftStaple(Link,mu);
}
//Same as Cshift for periodic BCs
static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
return PeriodicBC::CshiftLink(Link,mu,shift);
}
static inline bool isPeriodicGaugeField(void) { return true; }
};
@ -110,6 +115,11 @@ public:
return PeriodicBC::CovShiftBackward(Link, mu, field);
}
//If mu is a conjugate BC direction
//Out(x) = U^dag_\mu(x-mu) | x_\mu != 0
// = U^T_\mu(L-1) | x_\mu == 0
//else
//Out(x) = U^dag_\mu(x-mu mod L)
static inline GaugeLinkField
CovShiftIdentityBackward(const GaugeLinkField &Link, int mu)
{
@ -129,6 +139,13 @@ public:
return PeriodicBC::CovShiftIdentityForward(Link,mu);
}
//If mu is a conjugate BC direction
//Out(x) = S_\mu(x+mu) | x_\mu != L-1
// = S*_\mu(x+mu) | x_\mu == L-1
//else
//Out(x) = S_\mu(x+mu mod L)
//Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices
static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu)
{
assert(_conjDirs.size() == Nd);
@ -138,6 +155,27 @@ public:
return PeriodicBC::ShiftStaple(Link,mu);
}
//Boundary-aware C-shift of gauge links / gauge transformation matrices
//For conjugate BC direction
//shift = 1
//Out(x) = U_\mu(x+\hat\mu) | x_\mu != L-1
// = U*_\mu(0) | x_\mu == L-1
//shift = -1
//Out(x) = U_\mu(x-mu) | x_\mu != 0
// = U*_\mu(L-1) | x_\mu == 0
//else
//shift = 1
//Out(x) = U_\mu(x+\hat\mu mod L)
//shift = -1
//Out(x) = U_\mu(x-\hat\mu mod L)
static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
assert(_conjDirs.size() == Nd);
if(_conjDirs[mu])
return ConjugateBC::CshiftLink(Link,mu,shift);
else
return PeriodicBC::CshiftLink(Link,mu,shift);
}
static inline void setDirections(std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
static inline std::vector<int> getDirections(void) { return _conjDirs; }
static inline bool isPeriodicGaugeField(void) { return false; }

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

@ -88,6 +88,12 @@ namespace PeriodicBC {
return CovShiftBackward(Link,mu,arg);
}
//Boundary-aware C-shift of gauge links / gauge transformation matrices
template<class gauge> Lattice<gauge>
CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
{
return Cshift(Link, mu, shift);
}
}
@ -158,6 +164,9 @@ namespace ConjugateBC {
// std::cout<<"Gparity::CovCshiftBackward mu="<<mu<<std::endl;
return Cshift(tmp,mu,-1);// moves towards positive mu
}
//Out(x) = U^dag_\mu(x-mu) | x_\mu != 0
// = U^T_\mu(L-1) | x_\mu == 0
template<class gauge> Lattice<gauge>
CovShiftIdentityBackward(const Lattice<gauge> &Link, int mu) {
GridBase *grid = Link.Grid();
@ -176,6 +185,9 @@ namespace ConjugateBC {
return Link;
}
//Out(x) = S_\mu(x+\hat\mu) | x_\mu != L-1
// = S*_\mu(0) | x_\mu == L-1
//Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices
template<class gauge> Lattice<gauge>
ShiftStaple(const Lattice<gauge> &Link, int mu)
{
@ -208,6 +220,35 @@ namespace ConjugateBC {
return CovShiftBackward(Link,mu,arg);
}
//Boundary-aware C-shift of gauge links / gauge transformation matrices
//shift = 1
//Out(x) = U_\mu(x+\hat\mu) | x_\mu != L-1
// = U*_\mu(0) | x_\mu == L-1
//shift = -1
//Out(x) = U_\mu(x-mu) | x_\mu != 0
// = U*_\mu(L-1) | x_\mu == 0
template<class gauge> Lattice<gauge>
CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
{
GridBase *grid = Link.Grid();
int Lmu = grid->GlobalDimensions()[mu] - 1;
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);
return tmp;
}else if(shift == -1){
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
}
}

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

@ -40,27 +40,46 @@ public:
typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
static void GaugeLinkToLieAlgebraField(const std::vector<GaugeMat> &U,std::vector<GaugeMat> &A) {
for(int mu=0;mu<Nd;mu++){
Complex cmi(0.0,-1.0);
A[mu] = Ta(U[mu]) * cmi;
}
//A_\mu(x) = -i Ta(U_\mu(x) ) where Ta(U) = 1/2( U - U^dag ) - 1/2N tr(U - U^dag) is the traceless antihermitian part. This is an O(A^3) approximation to the logarithm of U
static void GaugeLinkToLieAlgebraField(const GaugeMat &U, GaugeMat &A) {
Complex cmi(0.0,-1.0);
A = Ta(U) * cmi;
}
static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu,int orthog) {
//The derivative of the Lie algebra field
static void DmuAmu(const std::vector<GaugeMat> &U, GaugeMat &dmuAmu,int orthog) {
GridBase* grid = U[0].Grid();
GaugeMat Ax(grid);
GaugeMat Axm1(grid);
GaugeMat Utmp(grid);
dmuAmu=Zero();
for(int mu=0;mu<Nd;mu++){
if ( mu != orthog ) {
dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
//Rather than define functionality to work out how the BCs apply to A_\mu we simply use the BC-aware Cshift to the gauge links and compute A_\mu(x) and A_\mu(x-1) separately
//Ax = A_\mu(x)
GaugeLinkToLieAlgebraField(U[mu], Ax);
//Axm1 = A_\mu(x_\mu-1)
Utmp = Gimpl::CshiftLink(U[mu], mu, -1);
GaugeLinkToLieAlgebraField(Utmp, Axm1);
//Derivative
dmuAmu = dmuAmu + Ax - Axm1;
}
}
}
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
//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) {
GridBase *grid = Umu.Grid();
GaugeMat xform(grid);
SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog);
}
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
//Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform, Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
GridBase *grid = Umu.Grid();
@ -122,27 +141,24 @@ public:
}
}
assert(0 && "Gauge fixing did not converge within the specified number of iterations");
};
static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) {
GridBase *grid = U[0].Grid();
std::vector<GaugeMat> A(Nd,grid);
GaugeMat g(grid);
GaugeLinkToLieAlgebraField(U,A);
ExpiAlphaDmuAmu(A,g,alpha,dmuAmu,orthog);
ExpiAlphaDmuAmu(U,g,alpha,dmuAmu,orthog);
Real vol = grid->gSites();
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
xform = g*xform ;
SU<Nc>::GaugeTransform(U,g);
SU<Nc>::GaugeTransform<Gimpl>(U,g);
return trG;
}
static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) {
GridBase *grid = U[0].Grid();
@ -157,11 +173,7 @@ public:
GaugeMat g(grid);
GaugeMat dmuAmu_p(grid);
std::vector<GaugeMat> A(Nd,grid);
GaugeLinkToLieAlgebraField(U,A);
DmuAmu(A,dmuAmu,orthog);
DmuAmu(U,dmuAmu,orthog);
std::vector<int> mask(Nd,1);
for(int mu=0;mu<Nd;mu++) if (mu==orthog) mask[mu]=0;
@ -205,16 +217,16 @@ public:
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
xform = g*xform ;
SU<Nc>::GaugeTransform(U,g);
SU<Nc>::GaugeTransform<Gimpl>(U,g);
return trG;
}
static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu,int orthog) {
static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &U,GaugeMat &g, Real alpha, GaugeMat &dmuAmu,int orthog) {
GridBase *grid = g.Grid();
Complex cialpha(0.0,-alpha);
GaugeMat ciadmam(grid);
DmuAmu(A,dmuAmu,orthog);
DmuAmu(U,dmuAmu,orthog);
ciadmam = dmuAmu*cialpha;
SU<Nc>::taExp(ciadmam,g);
}

24
Grid/qcd/utils/SUn.h
View File

@ -694,32 +694,32 @@ public:
* Adjoint rep gauge xform
*/
template<typename GaugeField,typename GaugeMat>
static void GaugeTransform( GaugeField &Umu, GaugeMat &g){
template<typename Gimpl>
static void GaugeTransform(typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){
GridBase *grid = Umu.Grid();
conformable(grid,g.Grid());
GaugeMat U(grid);
GaugeMat ag(grid); ag = adj(g);
typename Gimpl::GaugeLinkField U(grid);
typename Gimpl::GaugeLinkField ag(grid); ag = adj(g);
for(int mu=0;mu<Nd;mu++){
U= PeekIndex<LorentzIndex>(Umu,mu);
U = g*U*Cshift(ag, mu, 1);
U = g*U*Gimpl::CshiftLink(ag, mu, 1); //BC-aware
PokeIndex<LorentzIndex>(Umu,U,mu);
}
}
template<typename GaugeMat>
static void GaugeTransform( std::vector<GaugeMat> &U, GaugeMat &g){
template<typename Gimpl>
static void GaugeTransform( std::vector<typename Gimpl::GaugeLinkField> &U, typename Gimpl::GaugeLinkField &g){
GridBase *grid = g.Grid();
GaugeMat ag(grid); ag = adj(g);
typename Gimpl::GaugeLinkField ag(grid); ag = adj(g);
for(int mu=0;mu<Nd;mu++){
U[mu] = g*U[mu]*Cshift(ag, mu, 1);
U[mu] = g*U[mu]*Gimpl::CshiftLink(ag, mu, 1); //BC-aware
}
}
template<typename GaugeField,typename GaugeMat>
static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){
template<typename Gimpl>
static void RandomGaugeTransform(GridParallelRNG &pRNG, typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){
LieRandomize(pRNG,g,1.0);
GaugeTransform(Umu,g);
GaugeTransform<Gimpl>(Umu,g);
}
// Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )

185
tests/core/Test_fft_gfix.cc
View File

@ -29,14 +29,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/Grid.h>
using namespace Grid;
;
int main (int argc, char ** argv)
{
template<typename Gimpl>
void run(double alpha, bool do_fft_gfix){
std::vector<int> seeds({1,2,3,4});
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
Coordinate latt_size = GridDefaultLatt();
@ -55,10 +51,7 @@ int main (int argc, char ** argv)
FFT theFFT(&GRID);
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<<GridLogMessage << "Using alpha=" << alpha << std::endl;
// int coulomb_dir = -1;
int coulomb_dir = Nd-1;
@ -72,81 +65,165 @@ int main (int argc, char ** argv)
LatticeColourMatrix xform1(&GRID); // Gauge xform
LatticeColourMatrix xform2(&GRID); // Gauge xform
LatticeColourMatrix xform3(&GRID); // Gauge xform
//#########################################################################################
std::cout<< "*********************************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to Landau gauge with randomly transformed unit gauge configuration *" <<std::endl;
std::cout<< "*********************************************************************************************************" <<std::endl;
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
Uorg=Umu;
Real init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
//Apply a random gauge transformation to the unit gauge config
Urnd=Umu;
SU<Nc>::RandomGaugeTransform<Gimpl>(pRNG,Urnd,g);
SU<Nc>::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
Real alpha=0.1;
//Gauge fix the randomly transformed field
Umu = Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
// Check the gauge xform matrices
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform1);
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform1);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
Real plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
Uorg = Uorg - Umu;
std::cout << " Norm Difference "<< norm2(Uorg) << std::endl;
std::cout << " Norm "<< norm2(Umu) << std::endl;
std::cout << " Norm difference between a unit gauge configuration and the gauge fixed configuration "<< norm2(Uorg) << " (expect 0)" << std::endl;
std::cout << " Norm of gauge fixed configuration "<< norm2(Umu) << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "*************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with unit gauge configuration *" <<std::endl;
std::cout<< "*************************************************************************************" <<std::endl;
Umu=Urnd;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
Utmp=Urnd;
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform2);
Utmp = Utmp - Umu;
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
Umu=Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
//#########################################################################################
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform2);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to Landau gauge with random configuration **" <<std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing non-unit configuration *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "******************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with random configuration **" <<std::endl;
std::cout<< "******************************************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
//#########################################################################################
std::cout<< "*******************************************************************************************" <<std::endl;
std::cout<< "* Testing steepest descent fixing to coulomb gauge with random configuration *" <<std::endl;
std::cout<< "*******************************************************************************************" <<std::endl;
Umu=Urnd;
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::HotConfiguration(pRNG,Umu);
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,false,coulomb_dir);
std::cout << Umu<<std::endl;
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
//#########################################################################################
if(do_fft_gfix){
std::cout<< "*******************************************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge with random configuration *" <<std::endl;
std::cout<< "*******************************************************************************************" <<std::endl;
Umu=Urnd;
SU<Nc>::HotConfiguration(pRNG,Umu);
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<< init_plaq << std::endl;
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
}
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
double alpha=0.1; //step size
std::string gimpl = "periodic";
bool do_fft_gfix = true; //test fourier transformed gfix as well as steepest descent
for(int i=1;i<argc;i++){
std::string sarg(argv[i]);
if(sarg == "--gimpl"){
assert(i<argc-1 && "--gimpl option requires an argument");
gimpl = argv[i+1];
if(gimpl != "periodic" && gimpl != "conjugate")
assert(0 && "Invalid gimpl");
}else if(sarg == "--no-fft-gfix"){
std::cout << "Not doing the Fourier accelerated gauge fixing tests" << std::endl;
do_fft_gfix = false;
}else if(sarg == "--alpha"){
assert(i<argc-1 && "--alpha option requires an argument");
std::istringstream ss(argv[i+1]); ss >> alpha;
}
}
if(gimpl == "periodic"){
std::cout << GridLogMessage << "Using periodic boundary condition" << std::endl;
run<PeriodicGimplR>(alpha, do_fft_gfix);
}else{
std::vector<int> conjdirs = {1,1,0,0}; //test with 2 conjugate dirs and 2 not
std::cout << GridLogMessage << "Using complex conjugate boundary conditions in dimensions ";
for(int i=0;i<Nd;i++)
if(conjdirs[i])
std::cout << i << " ";
std::cout << std::endl;
ConjugateGimplR::setDirections(conjdirs);
run<ConjugateGimplR>(alpha, do_fft_gfix);
}
Grid_finalize();
}