mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-04 19:25:56 +01:00
Merge branch 'develop' into mbruno-eclover
This commit is contained in:
commit
77aa147ce5
@ -31,6 +31,7 @@ directory
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <map>
|
||||
|
||||
#include <pwd.h>
|
||||
@ -654,7 +655,8 @@ class IldgWriter : public ScidacWriter {
|
||||
// Fill ILDG header data struct
|
||||
//////////////////////////////////////////////////////
|
||||
ildgFormat ildgfmt ;
|
||||
ildgfmt.field = std::string("su3gauge");
|
||||
const std::string stNC = std::to_string( Nc ) ;
|
||||
ildgfmt.field = std::string("su"+stNC+"gauge");
|
||||
|
||||
if ( format == std::string("IEEE32BIG") ) {
|
||||
ildgfmt.precision = 32;
|
||||
@ -871,7 +873,8 @@ class IldgReader : public GridLimeReader {
|
||||
} else {
|
||||
|
||||
assert(found_ildgFormat);
|
||||
assert ( ildgFormat_.field == std::string("su3gauge") );
|
||||
const std::string stNC = std::to_string( Nc ) ;
|
||||
assert ( ildgFormat_.field == std::string("su"+stNC+"gauge") );
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////
|
||||
// Populate our Grid metadata as best we can
|
||||
@ -879,7 +882,7 @@ class IldgReader : public GridLimeReader {
|
||||
|
||||
std::ostringstream vers; vers << ildgFormat_.version;
|
||||
FieldMetaData_.hdr_version = vers.str();
|
||||
FieldMetaData_.data_type = std::string("4D_SU3_GAUGE_3X3");
|
||||
FieldMetaData_.data_type = std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC);
|
||||
|
||||
FieldMetaData_.nd=4;
|
||||
FieldMetaData_.dimension.resize(4);
|
||||
|
@ -6,8 +6,8 @@
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
|
||||
|
||||
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
|
||||
@ -182,8 +182,8 @@ class GaugeStatistics
|
||||
public:
|
||||
void operator()(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
|
||||
{
|
||||
header.link_trace=WilsonLoops<Impl>::linkTrace(data);
|
||||
header.plaquette =WilsonLoops<Impl>::avgPlaquette(data);
|
||||
header.link_trace = WilsonLoops<Impl>::linkTrace(data);
|
||||
header.plaquette = WilsonLoops<Impl>::avgPlaquette(data);
|
||||
}
|
||||
};
|
||||
typedef GaugeStatistics<PeriodicGimplD> PeriodicGaugeStatistics;
|
||||
@ -203,20 +203,24 @@ template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzCo
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
inline void reconstruct3(LorentzColourMatrix & cm)
|
||||
{
|
||||
const int x=0;
|
||||
const int y=1;
|
||||
const int z=2;
|
||||
assert( Nc < 4 && Nc > 1 ) ;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
|
||||
cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
|
||||
cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
|
||||
#if Nc == 2
|
||||
cm(mu)()(1,0) = -adj(cm(mu)()(0,y)) ;
|
||||
cm(mu)()(1,1) = adj(cm(mu)()(0,x)) ;
|
||||
#else
|
||||
const int x=0 , y=1 , z=2 ; // a little disinenuous labelling
|
||||
cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
|
||||
cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
|
||||
cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
// Some data types for intermediate storage
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, Nd >;
|
||||
template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, Nc-1>, Nd >;
|
||||
|
||||
typedef iLorentzColour2x3<Complex> LorentzColour2x3;
|
||||
typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
|
||||
@ -278,7 +282,6 @@ struct GaugeSimpleMunger{
|
||||
|
||||
template <class fobj, class sobj>
|
||||
struct GaugeSimpleUnmunger {
|
||||
|
||||
void operator()(sobj &in, fobj &out) {
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
for (int i = 0; i < Nc; i++) {
|
||||
@ -317,8 +320,8 @@ template<class fobj,class sobj>
|
||||
struct Gauge3x2munger{
|
||||
void operator() (fobj &in,sobj &out){
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
for(int i=0;i<2;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
for(int i=0;i<Nc-1;i++){
|
||||
for(int j=0;j<Nc;j++){
|
||||
out(mu)()(i,j) = in(mu)(i)(j);
|
||||
}}
|
||||
}
|
||||
@ -330,8 +333,8 @@ template<class fobj,class sobj>
|
||||
struct Gauge3x2unmunger{
|
||||
void operator() (sobj &in,fobj &out){
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
for(int i=0;i<2;i++){
|
||||
for(int j=0;j<3;j++){
|
||||
for(int i=0;i<Nc-1;i++){
|
||||
for(int j=0;j<Nc;j++){
|
||||
out(mu)(i)(j) = in(mu)()(i,j);
|
||||
}}
|
||||
}
|
||||
|
@ -9,6 +9,7 @@
|
||||
Author: Matt Spraggs <matthew.spraggs@gmail.com>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
|
||||
|
||||
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
|
||||
@ -30,6 +31,8 @@
|
||||
#ifndef GRID_NERSC_IO_H
|
||||
#define GRID_NERSC_IO_H
|
||||
|
||||
#include <string>
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
using namespace Grid;
|
||||
@ -145,15 +148,17 @@ public:
|
||||
|
||||
std::string format(header.floating_point);
|
||||
|
||||
int ieee32big = (format == std::string("IEEE32BIG"));
|
||||
int ieee32 = (format == std::string("IEEE32"));
|
||||
int ieee64big = (format == std::string("IEEE64BIG"));
|
||||
int ieee64 = (format == std::string("IEEE64") || format == std::string("IEEE64LITTLE"));
|
||||
const int ieee32big = (format == std::string("IEEE32BIG"));
|
||||
const int ieee32 = (format == std::string("IEEE32"));
|
||||
const int ieee64big = (format == std::string("IEEE64BIG"));
|
||||
const int ieee64 = (format == std::string("IEEE64") || \
|
||||
format == std::string("IEEE64LITTLE"));
|
||||
|
||||
uint32_t nersc_csum,scidac_csuma,scidac_csumb;
|
||||
// depending on datatype, set up munger;
|
||||
// munger is a function of <floating point, Real, data_type>
|
||||
if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
|
||||
const std::string stNC = std::to_string( Nc ) ;
|
||||
if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE") ) {
|
||||
if ( ieee32 || ieee32big ) {
|
||||
BinaryIO::readLatticeObject<vLorentzColourMatrixD, LorentzColour2x3F>
|
||||
(Umu,file,Gauge3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format,
|
||||
@ -164,7 +169,7 @@ public:
|
||||
(Umu,file,Gauge3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format,
|
||||
nersc_csum,scidac_csuma,scidac_csumb);
|
||||
}
|
||||
} else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) {
|
||||
} else if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC) ) {
|
||||
if ( ieee32 || ieee32big ) {
|
||||
BinaryIO::readLatticeObject<vLorentzColourMatrixD,LorentzColourMatrixF>
|
||||
(Umu,file,GaugeSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format,
|
||||
@ -209,27 +214,29 @@ public:
|
||||
template<class GaugeStats=PeriodicGaugeStatistics>
|
||||
static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
|
||||
std::string file,
|
||||
std::string ens_label = std::string("DWF"))
|
||||
std::string ens_label = std::string("DWF"),
|
||||
std::string ens_id = std::string("UKQCD"),
|
||||
unsigned int sequence_number = 1)
|
||||
{
|
||||
writeConfiguration(Umu,file,0,1,ens_label);
|
||||
writeConfiguration(Umu,file,0,1,ens_label,ens_id,sequence_number);
|
||||
}
|
||||
template<class GaugeStats=PeriodicGaugeStatistics>
|
||||
static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
|
||||
std::string file,
|
||||
int two_row,
|
||||
int bits32,
|
||||
std::string ens_label = std::string("DWF"))
|
||||
std::string ens_label = std::string("DWF"),
|
||||
std::string ens_id = std::string("UKQCD"),
|
||||
unsigned int sequence_number = 1)
|
||||
{
|
||||
typedef vLorentzColourMatrixD vobj;
|
||||
typedef typename vobj::scalar_object sobj;
|
||||
|
||||
FieldMetaData header;
|
||||
///////////////////////////////////////////
|
||||
// Following should become arguments
|
||||
///////////////////////////////////////////
|
||||
header.sequence_number = 1;
|
||||
header.ensemble_id = std::string("UKQCD");
|
||||
header.sequence_number = sequence_number;
|
||||
header.ensemble_id = ens_id;
|
||||
header.ensemble_label = ens_label;
|
||||
header.hdr_version = "1.0" ;
|
||||
|
||||
typedef LorentzColourMatrixD fobj3D;
|
||||
typedef LorentzColour2x3D fobj2D;
|
||||
@ -243,10 +250,14 @@ public:
|
||||
|
||||
uint64_t offset;
|
||||
|
||||
// Sod it -- always write 3x3 double
|
||||
header.floating_point = std::string("IEEE64BIG");
|
||||
header.data_type = std::string("4D_SU3_GAUGE_3x3");
|
||||
GaugeSimpleUnmunger<fobj3D,sobj> munge;
|
||||
// Sod it -- always write NcxNc double
|
||||
header.floating_point = std::string("IEEE64BIG");
|
||||
const std::string stNC = std::to_string( Nc ) ;
|
||||
if( two_row ) {
|
||||
header.data_type = std::string("4D_SU" + stNC + "_GAUGE" );
|
||||
} else {
|
||||
header.data_type = std::string("4D_SU" + stNC + "_GAUGE_" + stNC + "x" + stNC );
|
||||
}
|
||||
if ( grid->IsBoss() ) {
|
||||
truncate(file);
|
||||
offset = writeHeader(header,file);
|
||||
@ -254,8 +265,15 @@ public:
|
||||
grid->Broadcast(0,(void *)&offset,sizeof(offset));
|
||||
|
||||
uint32_t nersc_csum,scidac_csuma,scidac_csumb;
|
||||
BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
|
||||
nersc_csum,scidac_csuma,scidac_csumb);
|
||||
if( two_row ) {
|
||||
Gauge3x2unmunger<fobj2D,sobj> munge;
|
||||
BinaryIO::writeLatticeObject<vobj,fobj2D>(Umu,file,munge,offset,header.floating_point,
|
||||
nersc_csum,scidac_csuma,scidac_csumb);
|
||||
} else {
|
||||
GaugeSimpleUnmunger<fobj3D,sobj> munge;
|
||||
BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
|
||||
nersc_csum,scidac_csuma,scidac_csumb);
|
||||
}
|
||||
header.checksum = nersc_csum;
|
||||
if ( grid->IsBoss() ) {
|
||||
writeHeader(header,file);
|
||||
@ -287,8 +305,7 @@ public:
|
||||
header.plaquette=0.0;
|
||||
MachineCharacteristics(header);
|
||||
|
||||
uint64_t offset;
|
||||
|
||||
uint64_t offset;
|
||||
#ifdef RNG_RANLUX
|
||||
header.floating_point = std::string("UINT64");
|
||||
header.data_type = std::string("RANLUX48");
|
||||
@ -328,7 +345,7 @@ public:
|
||||
|
||||
GridBase *grid = parallel.Grid();
|
||||
|
||||
uint64_t offset = readHeader(file,grid,header);
|
||||
uint64_t offset = readHeader(file,grid,header);
|
||||
|
||||
FieldMetaData clone(header);
|
||||
|
||||
|
@ -68,9 +68,16 @@ public:
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Support for MADWF tricks
|
||||
///////////////////////////////////////////////////////////////
|
||||
RealD Mass(void) { return mass; };
|
||||
RealD Mass(void) { return (mass_plus + mass_minus) / 2.0; };
|
||||
RealD MassPlus(void) { return mass_plus; };
|
||||
RealD MassMinus(void) { return mass_minus; };
|
||||
void SetMass(RealD _mass) {
|
||||
mass=_mass;
|
||||
mass_plus=mass_minus=_mass;
|
||||
SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs
|
||||
} ;
|
||||
void SetMass(RealD _mass_plus, RealD _mass_minus) {
|
||||
mass_plus=_mass_plus;
|
||||
mass_minus=_mass_minus;
|
||||
SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs
|
||||
} ;
|
||||
void P(const FermionField &psi, FermionField &chi);
|
||||
@ -108,7 +115,7 @@ public:
|
||||
void MeooeDag5D (const FermionField &in, FermionField &out);
|
||||
|
||||
// protected:
|
||||
RealD mass;
|
||||
RealD mass_plus, mass_minus;
|
||||
|
||||
// Save arguments to SetCoefficientsInternal
|
||||
Vector<Coeff_t> _gamma;
|
||||
|
@ -47,7 +47,7 @@ CayleyFermion5D<Impl>::CayleyFermion5D(GaugeField &_Umu,
|
||||
FiveDimRedBlackGrid,
|
||||
FourDimGrid,
|
||||
FourDimRedBlackGrid,_M5,p),
|
||||
mass(_mass)
|
||||
mass_plus(_mass), mass_minus(_mass)
|
||||
{
|
||||
}
|
||||
|
||||
@ -209,8 +209,8 @@ void CayleyFermion5D<Impl>::M5D (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
Vector<Coeff_t> diag (Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] =mass;
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass_minus;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] =mass_plus;
|
||||
M5D(psi,chi,chi,lower,diag,upper);
|
||||
}
|
||||
template<class Impl>
|
||||
@ -220,8 +220,8 @@ void CayleyFermion5D<Impl>::Meooe5D (const FermionField &psi, FermionField &D
|
||||
Vector<Coeff_t> diag = bs;
|
||||
Vector<Coeff_t> upper= cs;
|
||||
Vector<Coeff_t> lower= cs;
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
upper[Ls-1]=-mass_minus*upper[Ls-1];
|
||||
lower[0] =-mass_plus*lower[0];
|
||||
M5D(psi,psi,Din,lower,diag,upper);
|
||||
}
|
||||
// FIXME Redunant with the above routine; check this and eliminate
|
||||
@ -235,8 +235,8 @@ template<class Impl> void CayleyFermion5D<Impl>::Meo5D (const FermionField &
|
||||
upper[i]=-ceo[i];
|
||||
lower[i]=-ceo[i];
|
||||
}
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
upper[Ls-1]=-mass_minus*upper[Ls-1];
|
||||
lower[0] =-mass_plus*lower[0];
|
||||
M5D(psi,psi,chi,lower,diag,upper);
|
||||
}
|
||||
template<class Impl>
|
||||
@ -250,8 +250,8 @@ void CayleyFermion5D<Impl>::Mooee (const FermionField &psi, FermionField &
|
||||
upper[i]=-cee[i];
|
||||
lower[i]=-cee[i];
|
||||
}
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
upper[Ls-1]=-mass_minus*upper[Ls-1];
|
||||
lower[0] =-mass_plus*lower[0];
|
||||
M5D(psi,psi,chi,lower,diag,upper);
|
||||
}
|
||||
template<class Impl>
|
||||
@ -266,9 +266,9 @@ void CayleyFermion5D<Impl>::MooeeDag (const FermionField &psi, FermionField &
|
||||
// Assemble the 5d matrix
|
||||
if ( s==0 ) {
|
||||
upper[s] = -cee[s+1] ;
|
||||
lower[s] = mass*cee[Ls-1];
|
||||
lower[s] = mass_minus*cee[Ls-1];
|
||||
} else if ( s==(Ls-1)) {
|
||||
upper[s] = mass*cee[0];
|
||||
upper[s] = mass_plus*cee[0];
|
||||
lower[s] = -cee[s-1];
|
||||
} else {
|
||||
upper[s]=-cee[s+1];
|
||||
@ -291,8 +291,8 @@ void CayleyFermion5D<Impl>::M5Ddag (const FermionField &psi, FermionField &chi)
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0);
|
||||
Vector<Coeff_t> lower(Ls,-1.0);
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
upper[Ls-1]=-mass_plus*upper[Ls-1];
|
||||
lower[0] =-mass_minus*lower[0];
|
||||
M5Ddag(psi,chi,chi,lower,diag,upper);
|
||||
}
|
||||
|
||||
@ -307,9 +307,9 @@ void CayleyFermion5D<Impl>::MeooeDag5D (const FermionField &psi, FermionField
|
||||
for (int s=0;s<Ls;s++){
|
||||
if ( s== 0 ) {
|
||||
upper[s] = cs[s+1];
|
||||
lower[s] =-mass*cs[Ls-1];
|
||||
lower[s] =-mass_minus*cs[Ls-1];
|
||||
} else if ( s==(Ls-1) ) {
|
||||
upper[s] =-mass*cs[0];
|
||||
upper[s] =-mass_plus*cs[0];
|
||||
lower[s] = cs[s-1];
|
||||
} else {
|
||||
upper[s] = cs[s+1];
|
||||
@ -552,7 +552,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
|
||||
|
||||
lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
|
||||
|
||||
leem[i]=mass*cee[Ls-1]/bee[0];
|
||||
leem[i]=mass_minus*cee[Ls-1]/bee[0];
|
||||
for(int j=0;j<i;j++) {
|
||||
assert(bee[j+1]!=Coeff_t(0.0));
|
||||
leem[i]*= aee[j]/bee[j+1];
|
||||
@ -560,7 +560,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
|
||||
|
||||
uee[i] =-aee[i]/bee[i]; // up-diag entry on the ith row
|
||||
|
||||
ueem[i]=mass;
|
||||
ueem[i]=mass_plus;
|
||||
for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
|
||||
ueem[i]*= aee[0]/bee[0];
|
||||
|
||||
@ -573,7 +573,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
|
||||
}
|
||||
|
||||
{
|
||||
Coeff_t delta_d=mass*cee[Ls-1];
|
||||
Coeff_t delta_d=mass_minus*cee[Ls-1];
|
||||
for(int j=0;j<Ls-1;j++) {
|
||||
assert(bee[j] != Coeff_t(0.0));
|
||||
delta_d *= cee[j]/bee[j];
|
||||
@ -642,6 +642,10 @@ void CayleyFermion5D<Impl>::ContractConservedCurrent( PropagatorField &q_in_1,
|
||||
Current curr_type,
|
||||
unsigned int mu)
|
||||
{
|
||||
|
||||
assert(mass_plus == mass_minus);
|
||||
RealD mass = mass_plus;
|
||||
|
||||
#if (!defined(GRID_HIP))
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
@ -777,6 +781,8 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
|
||||
assert(mu>=0);
|
||||
assert(mu<Nd);
|
||||
|
||||
assert(mass_plus == mass_minus);
|
||||
RealD mass = mass_plus;
|
||||
|
||||
#if 0
|
||||
int tshift = (mu == Nd-1) ? 1 : 0;
|
||||
|
@ -65,8 +65,11 @@ CompactWilsonCloverFermion<Impl, CloverHelpers>::CompactWilsonCloverFermion(Gaug
|
||||
csw_r /= clover_anisotropy.xi_0;
|
||||
|
||||
ImportGauge(_Umu);
|
||||
if (open_boundaries)
|
||||
if (open_boundaries) {
|
||||
this->BoundaryMaskEven.Checkerboard() = Even;
|
||||
this->BoundaryMaskOdd.Checkerboard() = Odd;
|
||||
CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
|
@ -4,12 +4,13 @@ Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/WilsonFermion.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2022
|
||||
|
||||
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Fabian Joswig <fabian.joswig@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
|
||||
@ -599,11 +600,47 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
|
||||
Current curr_type,
|
||||
unsigned int mu)
|
||||
{
|
||||
if(curr_type != Current::Vector)
|
||||
{
|
||||
std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
conformable(_grid, q_in_1.Grid());
|
||||
conformable(_grid, q_in_2.Grid());
|
||||
conformable(_grid, q_out.Grid());
|
||||
assert(0);
|
||||
auto UGrid= this->GaugeGrid();
|
||||
|
||||
PropagatorField tmp_shifted(UGrid);
|
||||
PropagatorField g5Lg5(UGrid);
|
||||
PropagatorField R(UGrid);
|
||||
PropagatorField gmuR(UGrid);
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT,
|
||||
};
|
||||
Gamma gmu=Gamma(Gmu[mu]);
|
||||
|
||||
g5Lg5=g5*q_in_1*g5;
|
||||
tmp_shifted=Cshift(q_in_2,mu,1);
|
||||
Impl::multLinkField(R,this->Umu,tmp_shifted,mu);
|
||||
gmuR=gmu*R;
|
||||
|
||||
q_out=adj(g5Lg5)*R;
|
||||
q_out-=adj(g5Lg5)*gmuR;
|
||||
|
||||
tmp_shifted=Cshift(q_in_1,mu,1);
|
||||
Impl::multLinkField(g5Lg5,this->Umu,tmp_shifted,mu);
|
||||
g5Lg5=g5*g5Lg5*g5;
|
||||
R=q_in_2;
|
||||
gmuR=gmu*R;
|
||||
|
||||
q_out-=adj(g5Lg5)*R;
|
||||
q_out-=adj(g5Lg5)*gmuR;
|
||||
}
|
||||
|
||||
|
||||
@ -617,9 +654,51 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
|
||||
unsigned int tmax,
|
||||
ComplexField &lattice_cmplx)
|
||||
{
|
||||
if(curr_type != Current::Vector)
|
||||
{
|
||||
std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
int tshift = (mu == Nd-1) ? 1 : 0;
|
||||
unsigned int LLt = GridDefaultLatt()[Tp];
|
||||
conformable(_grid, q_in.Grid());
|
||||
conformable(_grid, q_out.Grid());
|
||||
assert(0);
|
||||
auto UGrid= this->GaugeGrid();
|
||||
|
||||
PropagatorField tmp(UGrid);
|
||||
PropagatorField Utmp(UGrid);
|
||||
PropagatorField L(UGrid);
|
||||
PropagatorField zz (UGrid);
|
||||
zz=Zero();
|
||||
LatticeInteger lcoor(UGrid); LatticeCoordinate(lcoor,Nd-1);
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT,
|
||||
};
|
||||
Gamma gmu=Gamma(Gmu[mu]);
|
||||
|
||||
tmp = Cshift(q_in,mu,1);
|
||||
Impl::multLinkField(Utmp,this->Umu,tmp,mu);
|
||||
tmp = ( Utmp*lattice_cmplx - gmu*Utmp*lattice_cmplx ); // Forward hop
|
||||
tmp = where((lcoor>=tmin),tmp,zz); // Mask the time
|
||||
q_out = where((lcoor<=tmax),tmp,zz); // Position of current complicated
|
||||
|
||||
tmp = q_in *lattice_cmplx;
|
||||
tmp = Cshift(tmp,mu,-1);
|
||||
Impl::multLinkField(Utmp,this->Umu,tmp,mu+Nd); // Adjoint link
|
||||
tmp = -( Utmp + gmu*Utmp );
|
||||
// Mask the time
|
||||
if (tmax == LLt - 1 && tshift == 1){ // quick fix to include timeslice 0 if tmax + tshift is over the last timeslice
|
||||
unsigned int t0 = 0;
|
||||
tmp = where(((lcoor==t0) || (lcoor>=tmin+tshift)),tmp,zz);
|
||||
} else {
|
||||
tmp = where((lcoor>=tmin+tshift),tmp,zz);
|
||||
}
|
||||
q_out+= where((lcoor<=tmax+tshift),tmp,zz); // Position of current complicated
|
||||
}
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
@ -55,12 +55,12 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
|
||||
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);
|
||||
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) {
|
||||
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) {
|
||||
|
||||
GridBase *grid = Umu.Grid();
|
||||
|
||||
@ -122,6 +122,8 @@ public:
|
||||
|
||||
}
|
||||
}
|
||||
std::cout << GridLogError << "Gauge fixing did not converge in " << maxiter << " iterations." << std::endl;
|
||||
if (err_on_no_converge) assert(0);
|
||||
};
|
||||
static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
|
||||
GridBase *grid = U[0].Grid();
|
||||
|
@ -125,7 +125,6 @@ public:
|
||||
return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
|
||||
}
|
||||
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// average over all x,y,z the temporal loop
|
||||
//////////////////////////////////////////////////
|
||||
@ -165,7 +164,7 @@ public:
|
||||
|
||||
double vol = Umu.Grid()->gSites();
|
||||
|
||||
return p.real() / vol / 4.0 / 3.0;
|
||||
return p.real() / vol / (4.0 * Nc ) ;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
|
@ -81,8 +81,8 @@ int main (int argc, char ** argv)
|
||||
Vector<Coeff_t> diag = Dw.bs;
|
||||
Vector<Coeff_t> upper= Dw.cs;
|
||||
Vector<Coeff_t> lower= Dw.cs;
|
||||
upper[Ls-1]=-Dw.mass*upper[Ls-1];
|
||||
lower[0] =-Dw.mass*lower[0];
|
||||
upper[Ls-1]=-Dw.mass_minus*upper[Ls-1];
|
||||
lower[0] =-Dw.mass_plus*lower[0];
|
||||
|
||||
LatticeFermion r_eo(FGrid);
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
|
@ -159,7 +159,7 @@ case ${ac_ZMOBIUS} in
|
||||
esac
|
||||
############### Nc
|
||||
AC_ARG_ENABLE([Nc],
|
||||
[AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])],
|
||||
[AC_HELP_STRING([--enable-Nc=2|3|4|5], [enable number of colours])],
|
||||
[ac_Nc=${enable_Nc}], [ac_Nc=3])
|
||||
|
||||
case ${ac_Nc} in
|
||||
|
@ -147,7 +147,7 @@ int main (int argc, char ** argv)
|
||||
Complex p = TensorRemove(Tp);
|
||||
std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
|
||||
|
||||
Complex LinkTraceScale(1.0/vol/4.0/3.0);
|
||||
Complex LinkTraceScale(1.0/vol/4.0/(Real)Nc);
|
||||
TComplex Tl = sum(LinkTrace);
|
||||
Complex l = TensorRemove(Tl);
|
||||
std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
|
||||
@ -157,8 +157,10 @@ int main (int argc, char ** argv)
|
||||
Complex ll= TensorRemove(TcP);
|
||||
std::cout<<GridLogMessage << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
|
||||
|
||||
std::string clone2x3("./ckpoint_clone2x3.4000");
|
||||
std::string clone3x3("./ckpoint_clone3x3.4000");
|
||||
const string stNc = to_string( Nc ) ;
|
||||
const string stNcM1 = to_string( Nc-1 ) ;
|
||||
std::string clone2x3("./ckpoint_clone"+stNcM1+"x"+stNc+".4000");
|
||||
std::string clone3x3("./ckpoint_clone"+stNc+"x"+stNc+".4000");
|
||||
|
||||
NerscIO::writeConfiguration(Umu,clone3x3,0,precision32);
|
||||
NerscIO::writeConfiguration(Umu,clone2x3,1,precision32);
|
||||
|
@ -9,6 +9,7 @@ Copyright (C) 2015
|
||||
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||
Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
|
||||
|
||||
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
|
||||
@ -42,14 +43,14 @@ directory
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
;
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
std::vector<int> latt({4, 4, 4, 8});
|
||||
GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
|
||||
latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
|
||||
GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
|
||||
|
||||
@ -60,15 +61,19 @@ int main(int argc, char** argv) {
|
||||
<< std::endl;
|
||||
SU2::printGenerators();
|
||||
std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl;
|
||||
|
||||
// guard as this code fails to compile for Nc != 3
|
||||
#if (Nc == 3)
|
||||
|
||||
SU2Adjoint::printGenerators();
|
||||
SU2::testGenerators();
|
||||
SU2Adjoint::testGenerators();
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
SU3::printGenerators();
|
||||
std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl;
|
||||
SU3Adjoint::printGenerators();
|
||||
@ -111,12 +116,10 @@ int main(int argc, char** argv) {
|
||||
|
||||
// AdjointRepresentation has the predefined number of colours Nc
|
||||
// Representations<FundamentalRepresentation, AdjointRepresentation, TwoIndexSymmetricRepresentation> RepresentationTypes(grid);
|
||||
|
||||
|
||||
LatticeGaugeField U(grid), V(grid);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V);
|
||||
|
||||
|
||||
// Adjoint representation
|
||||
// Test group structure
|
||||
// (U_f * V_f)_r = U_r * V_r
|
||||
@ -127,17 +130,17 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeMatrix Vmu = peekLorentz(V,mu);
|
||||
pokeLorentz(UV,Umu*Vmu, mu);
|
||||
}
|
||||
|
||||
|
||||
AdjRep.update_representation(UV);
|
||||
typename AdjointRep<Nc>::LatticeField UVr = AdjRep.U; // (U_f * V_f)_r
|
||||
|
||||
|
||||
|
||||
|
||||
AdjRep.update_representation(U);
|
||||
typename AdjointRep<Nc>::LatticeField Ur = AdjRep.U; // U_r
|
||||
|
||||
|
||||
AdjRep.update_representation(V);
|
||||
typename AdjointRep<Nc>::LatticeField Vr = AdjRep.U; // V_r
|
||||
|
||||
|
||||
typename AdjointRep<Nc>::LatticeField UrVr(grid);
|
||||
UrVr = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
@ -145,10 +148,10 @@ int main(int argc, char** argv) {
|
||||
typename AdjointRep<Nc>::LatticeMatrix Vrmu = peekLorentz(Vr,mu);
|
||||
pokeLorentz(UrVr,Urmu*Vrmu, mu);
|
||||
}
|
||||
|
||||
|
||||
typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
|
||||
std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Adjoint representation) : " << norm2(Diff_check) << std::endl;
|
||||
|
||||
|
||||
// Check correspondence of algebra and group transformations
|
||||
// Create a random vector
|
||||
SU3::LatticeAlgebraVector h_adj(grid);
|
||||
@ -156,32 +159,31 @@ int main(int argc, char** argv) {
|
||||
random(gridRNG,h_adj);
|
||||
h_adj = real(h_adj);
|
||||
SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
|
||||
|
||||
|
||||
// Re-extract h_adj
|
||||
SU3::LatticeAlgebraVector h_adj2(grid);
|
||||
SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
|
||||
SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2;
|
||||
std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
|
||||
|
||||
|
||||
// Exponentiate
|
||||
typename AdjointRep<Nc>::LatticeMatrix Uadj(grid);
|
||||
Uadj = expMat(Ar, 1.0, 16);
|
||||
|
||||
|
||||
typename AdjointRep<Nc>::LatticeMatrix uno(grid);
|
||||
uno = 1.0;
|
||||
// Check matrix Uadj, must be real orthogonal
|
||||
typename AdjointRep<Nc>::LatticeMatrix Ucheck = Uadj - conjugate(Uadj);
|
||||
std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck)
|
||||
<< std::endl;
|
||||
|
||||
<< std::endl;
|
||||
|
||||
Ucheck = Uadj * adj(Uadj) - uno;
|
||||
std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck)
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
Ucheck = adj(Uadj) * Uadj - uno;
|
||||
std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
<< std::endl;
|
||||
|
||||
// Construct the fundamental matrix in the group
|
||||
SU3::LatticeMatrix Af(grid);
|
||||
SU3::FundamentalLieAlgebraMatrix(h_adj,Af);
|
||||
@ -193,72 +195,65 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeMatrix UnitCheck(grid);
|
||||
UnitCheck = Ufund * adj(Ufund) - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
UnitCheck = adj(Ufund) * Ufund - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck)
|
||||
<< std::endl;
|
||||
|
||||
<< std::endl;
|
||||
|
||||
// Tranform to the adjoint representation
|
||||
U = Zero(); // fill this with only one direction
|
||||
pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields
|
||||
|
||||
|
||||
AdjRep.update_representation(U);
|
||||
Ur = AdjRep.U; // U_r
|
||||
typename AdjointRep<Nc>::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj
|
||||
|
||||
|
||||
typename AdjointRep<Nc>::LatticeMatrix Diff_check_mat = Ur0 - Uadj;
|
||||
std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl;
|
||||
|
||||
|
||||
|
||||
|
||||
// TwoIndexRep tests
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "* eS^{ij} base for SU(2)" << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Dimension of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
|
||||
SU2TwoIndexSymm::printBase();
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
|
||||
SU2TwoIndexSymm::printGenerators();
|
||||
std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl;
|
||||
std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl;
|
||||
SU2TwoIndexSymm::testGenerators();
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "* eAS^{ij} base for SU(2)" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
|
||||
SU2TwoIndexAntiSymm::printBase();
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
|
||||
SU2TwoIndexAntiSymm::printGenerators();
|
||||
std::cout << GridLogMessage << "Test of Two Index anti-Symmetric Generators: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
|
||||
SU2TwoIndexAntiSymm::testGenerators();
|
||||
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
// Projectors
|
||||
SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid);
|
||||
random(gridRNG,Gauss2);
|
||||
@ -276,13 +271,13 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeAlgebraVector diff2 = ha - hb;
|
||||
std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
// Projectors
|
||||
SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid);
|
||||
random(gridRNG,Gauss2a);
|
||||
@ -300,11 +295,11 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeAlgebraVector diff2a = ha - hb;
|
||||
std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "Two index Symmetric: Checking Group Structure"
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
// Testing HMC representation classes
|
||||
TwoIndexRep< Nc, Symmetric > TIndexRep(grid);
|
||||
|
||||
@ -313,7 +308,7 @@ int main(int argc, char** argv) {
|
||||
LatticeGaugeField U2(grid), V2(grid);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
|
||||
|
||||
|
||||
LatticeGaugeField UV2(grid);
|
||||
UV2 = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
@ -321,16 +316,16 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
|
||||
pokeLorentz(UV2,Umu2*Vmu2, mu);
|
||||
}
|
||||
|
||||
|
||||
TIndexRep.update_representation(UV2);
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeField UVr2 = TIndexRep.U; // (U_f * V_f)_r
|
||||
|
||||
|
||||
TIndexRep.update_representation(U2);
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2 = TIndexRep.U; // U_r
|
||||
|
||||
|
||||
TIndexRep.update_representation(V2);
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeField Vr2 = TIndexRep.U; // V_r
|
||||
|
||||
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2Vr2(grid);
|
||||
Ur2Vr2 = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
@ -338,11 +333,11 @@ int main(int argc, char** argv) {
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Vrmu2 = peekLorentz(Vr2,mu);
|
||||
pokeLorentz(Ur2Vr2,Urmu2*Vrmu2, mu);
|
||||
}
|
||||
|
||||
|
||||
typename TwoIndexRep< Nc, Symmetric >::LatticeField Diff_check2 = UVr2 - Ur2Vr2;
|
||||
std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index Symmetric): " << norm2(Diff_check2) << std::endl;
|
||||
|
||||
|
||||
|
||||
|
||||
// Check correspondence of algebra and group transformations
|
||||
// Create a random vector
|
||||
SU3::LatticeAlgebraVector h_sym(grid);
|
||||
@ -350,34 +345,31 @@ int main(int argc, char** argv) {
|
||||
random(gridRNG,h_sym);
|
||||
h_sym = real(h_sym);
|
||||
SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
|
||||
|
||||
|
||||
// Re-extract h_sym
|
||||
SU3::LatticeAlgebraVector h_sym2(grid);
|
||||
SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
|
||||
SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
|
||||
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
|
||||
|
||||
|
||||
|
||||
// Exponentiate
|
||||
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix U2iS(grid);
|
||||
U2iS = expMat(Ar_sym, 1.0, 16);
|
||||
|
||||
|
||||
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix uno2iS(grid);
|
||||
uno2iS = 1.0;
|
||||
// Check matrix U2iS, must be real orthogonal
|
||||
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ucheck2iS = U2iS - conjugate(U2iS);
|
||||
std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iS)
|
||||
<< std::endl;
|
||||
|
||||
<< std::endl;
|
||||
|
||||
Ucheck2iS = U2iS * adj(U2iS) - uno2iS;
|
||||
std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iS)
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
Ucheck2iS = adj(U2iS) * U2iS - uno2iS;
|
||||
std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iS)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
|
||||
<< std::endl;
|
||||
|
||||
// Construct the fundamental matrix in the group
|
||||
SU3::LatticeMatrix Af_sym(grid);
|
||||
SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
|
||||
@ -386,147 +378,137 @@ int main(int argc, char** argv) {
|
||||
SU3::LatticeMatrix UnitCheck2(grid);
|
||||
UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
|
||||
<< std::endl;
|
||||
<< std::endl;
|
||||
UnitCheck2 = adj(Ufund2) * Ufund2 - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
<< std::endl;
|
||||
|
||||
|
||||
// Tranform to the 2Index Sym representation
|
||||
U = Zero(); // fill this with only one direction
|
||||
pokeLorentz(U,Ufund2,0); // the representation transf acts on full gauge fields
|
||||
|
||||
|
||||
TIndexRep.update_representation(U);
|
||||
Ur2 = TIndexRep.U; // U_r
|
||||
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ur02 = peekLorentz(Ur2,0); // this should be the same as U2iS
|
||||
|
||||
|
||||
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Diff_check_mat2 = Ur02 - U2iS;
|
||||
std::cout << GridLogMessage << "Projections structure check group difference (Two Index Symmetric): " << norm2(Diff_check_mat2) << std::endl;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
if (TwoIndexRep<Nc, AntiSymmetric >::Dimension != 1){
|
||||
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
std::cout << GridLogMessage << "*********************************************"
|
||||
<< std::endl;
|
||||
|
||||
|
||||
std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure"
|
||||
<< std::endl;
|
||||
// Testing HMC representation classes
|
||||
TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid);
|
||||
std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure"
|
||||
<< std::endl;
|
||||
// Testing HMC representation classes
|
||||
TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid);
|
||||
|
||||
|
||||
// Test group structure
|
||||
// (U_f * V_f)_r = U_r * V_r
|
||||
LatticeGaugeField U2A(grid), V2A(grid);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
|
||||
// Test group structure
|
||||
// (U_f * V_f)_r = U_r * V_r
|
||||
LatticeGaugeField U2A(grid), V2A(grid);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
|
||||
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
|
||||
|
||||
LatticeGaugeField UV2A(grid);
|
||||
UV2A = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
|
||||
SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
|
||||
pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
|
||||
LatticeGaugeField UV2A(grid);
|
||||
UV2A = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
|
||||
SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
|
||||
pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
|
||||
}
|
||||
|
||||
TIndexRep.update_representation(UV2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r
|
||||
|
||||
TIndexRep.update_representation(U2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r
|
||||
|
||||
TIndexRep.update_representation(V2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid);
|
||||
Ur2Vr2A = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu);
|
||||
pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu);
|
||||
}
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A;
|
||||
std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index anti-Symmetric): " << norm2(Diff_check2A) << std::endl;
|
||||
|
||||
|
||||
// Check correspondence of algebra and group transformations
|
||||
// Create a random vector
|
||||
SU3::LatticeAlgebraVector h_Asym(grid);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
|
||||
random(gridRNG,h_Asym);
|
||||
h_Asym = real(h_Asym);
|
||||
SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
|
||||
|
||||
// Re-extract h_sym
|
||||
SU3::LatticeAlgebraVector h_Asym2(grid);
|
||||
SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
|
||||
SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
|
||||
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
|
||||
|
||||
|
||||
// Exponentiate
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid);
|
||||
U2iAS = expMat(Ar_Asym, 1.0, 16);
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid);
|
||||
uno2iAS = 1.0;
|
||||
// Check matrix U2iS, must be real orthogonal
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS);
|
||||
std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
|
||||
Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS;
|
||||
std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS;
|
||||
std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
|
||||
// Construct the fundamental matrix in the group
|
||||
SU3::LatticeMatrix Af_Asym(grid);
|
||||
SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
|
||||
SU3::LatticeMatrix Ufund2A(grid);
|
||||
Ufund2A = expMat(Af_Asym, 1.0, 16);
|
||||
SU3::LatticeMatrix UnitCheck2A(grid);
|
||||
UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
|
||||
<< std::endl;
|
||||
UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
// Tranform to the 2Index Sym representation
|
||||
U = Zero(); // fill this with only one direction
|
||||
pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields
|
||||
|
||||
TIndexRepA.update_representation(U);
|
||||
Ur2A = TIndexRepA.U; // U_r
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS;
|
||||
std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl;
|
||||
|
||||
} else {
|
||||
std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests "
|
||||
"because representation is trivial (dim = 1)"
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
TIndexRep.update_representation(UV2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r
|
||||
|
||||
TIndexRep.update_representation(U2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r
|
||||
|
||||
TIndexRep.update_representation(V2A);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid);
|
||||
Ur2Vr2A = Zero();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu);
|
||||
pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu);
|
||||
}
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A;
|
||||
std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index anti-Symmetric): " << norm2(Diff_check2A) << std::endl;
|
||||
|
||||
|
||||
// Check correspondence of algebra and group transformations
|
||||
// Create a random vector
|
||||
SU3::LatticeAlgebraVector h_Asym(grid);
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
|
||||
random(gridRNG,h_Asym);
|
||||
h_Asym = real(h_Asym);
|
||||
SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
|
||||
|
||||
// Re-extract h_sym
|
||||
SU3::LatticeAlgebraVector h_Asym2(grid);
|
||||
SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
|
||||
SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
|
||||
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
|
||||
|
||||
|
||||
// Exponentiate
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid);
|
||||
U2iAS = expMat(Ar_Asym, 1.0, 16);
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid);
|
||||
uno2iAS = 1.0;
|
||||
// Check matrix U2iS, must be real orthogonal
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS);
|
||||
std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
|
||||
Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS;
|
||||
std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS;
|
||||
std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
|
||||
// Construct the fundamental matrix in the group
|
||||
SU3::LatticeMatrix Af_Asym(grid);
|
||||
SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
|
||||
SU3::LatticeMatrix Ufund2A(grid);
|
||||
Ufund2A = expMat(Af_Asym, 1.0, 16);
|
||||
SU3::LatticeMatrix UnitCheck2A(grid);
|
||||
UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
|
||||
<< std::endl;
|
||||
UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f;
|
||||
std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A)
|
||||
<< std::endl;
|
||||
|
||||
|
||||
// Tranform to the 2Index Sym representation
|
||||
U = Zero(); // fill this with only one direction
|
||||
pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields
|
||||
|
||||
TIndexRepA.update_representation(U);
|
||||
Ur2A = TIndexRepA.U; // U_r
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS
|
||||
|
||||
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS;
|
||||
std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl;
|
||||
|
||||
} else {
|
||||
std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests "
|
||||
"because representation is trivial (dim = 1)"
|
||||
<< std::endl;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
@ -122,14 +122,15 @@ int main (int argc, char ** argv)
|
||||
std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl;
|
||||
tmp = U*adj(U) - ident;
|
||||
std::cout << "Unitarity check before projection " << norm2(tmp)<<std::endl;
|
||||
|
||||
#if (Nc == 3)
|
||||
ProjectSU3(U);
|
||||
detU= Determinant(U) ;
|
||||
detU= detU -1.0;
|
||||
std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl;
|
||||
tmp = U*adj(U) - ident;
|
||||
std::cout << "Unitarity check after projection " << norm2(tmp)<<std::endl;
|
||||
|
||||
#endif
|
||||
|
||||
ProjectSUn(UU);
|
||||
detUU= Determinant(UU);
|
||||
detUU= detUU -1.0;
|
||||
|
253
tests/core/Test_wilson_conserved_current.cc
Normal file
253
tests/core/Test_wilson_conserved_current.cc
Normal file
@ -0,0 +1,253 @@
|
||||
/*************************************************************************************
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_cayley_cg.cc
|
||||
|
||||
Copyright (C) 2022
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Fabian Joswig <fabian.joswig@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 */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
|
||||
template<class What>
|
||||
void TestConserved(What & Dw,
|
||||
LatticeGaugeField &Umu,
|
||||
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
|
||||
GridParallelRNG *RNG4);
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT,
|
||||
Gamma::Algebra::Gamma5
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
int threads = GridThread::GetThreads();
|
||||
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
GridParallelRNG RNG4(UGrid);
|
||||
std::vector<int> seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
if( argc > 1 && argv[1][0] != '-' )
|
||||
{
|
||||
std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
|
||||
FieldMetaData header;
|
||||
NerscIO::readConfiguration(Umu, header, argv[1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
|
||||
SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
}
|
||||
|
||||
typename WilsonCloverFermionR::ImplParams params;
|
||||
WilsonAnisotropyCoefficients anis;
|
||||
RealD mass = 0.1;
|
||||
RealD csw_r = 1.0;
|
||||
RealD csw_t = 1.0;
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonFermionR Dw(Umu,*UGrid,*UrbGrid,mass,params);
|
||||
TestConserved<WilsonFermionR>(Dw,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonCloverFermionR Dwc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonCloverFermionR>(Dwc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"CompactWilsonCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
CompactWilsonCloverFermionR Dwcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonCloverFermionR>(Dwcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonExpCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonExpCloverFermionR Dewc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonExpCloverFermionR>(Dewc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"CompactWilsonExpCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
CompactWilsonExpCloverFermionR Dewcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonExpCloverFermionR>(Dewcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
|
||||
|
||||
template<class Action>
|
||||
void TestConserved(Action & Dw,
|
||||
LatticeGaugeField &Umu,
|
||||
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
|
||||
GridParallelRNG *RNG4)
|
||||
{
|
||||
LatticePropagator phys_src(UGrid);
|
||||
LatticePropagator seqsrc(UGrid);
|
||||
LatticePropagator prop4(UGrid);
|
||||
LatticePropagator Vector_mu(UGrid);
|
||||
LatticeComplex SV (UGrid);
|
||||
LatticeComplex VV (UGrid);
|
||||
LatticePropagator seqprop(UGrid);
|
||||
|
||||
SpinColourMatrix kronecker; kronecker=1.0;
|
||||
Coordinate coor({0,0,0,0});
|
||||
phys_src=Zero();
|
||||
pokeSite(kronecker,phys_src,coor);
|
||||
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-16,100000);
|
||||
SchurRedBlackDiagTwoSolve<LatticeFermion> schur(CG);
|
||||
ZeroGuesser<LatticeFermion> zpg;
|
||||
for(int s=0;s<Nd;s++){
|
||||
for(int c=0;c<Nc;c++){
|
||||
LatticeFermion src4 (UGrid);
|
||||
PropToFerm<Action>(src4,phys_src,s,c);
|
||||
|
||||
LatticeFermion result4(UGrid); result4=Zero();
|
||||
schur(Dw,src4,result4,zpg);
|
||||
std::cout<<GridLogMessage<<"spin "<<s<<" color "<<c<<" norm2(sourc4d) "<<norm2(src4)
|
||||
<<" norm2(result4d) "<<norm2(result4)<<std::endl;
|
||||
FermToProp<Action>(prop4,result4,s,c);
|
||||
}
|
||||
}
|
||||
|
||||
auto curr = Current::Vector;
|
||||
const int mu_J=0;
|
||||
const int t_J=0;
|
||||
|
||||
LatticeComplex ph (UGrid); ph=1.0;
|
||||
|
||||
Dw.SeqConservedCurrent(prop4,
|
||||
seqsrc,
|
||||
phys_src,
|
||||
curr,
|
||||
mu_J,
|
||||
t_J,
|
||||
t_J,// whole lattice
|
||||
ph);
|
||||
|
||||
for(int s=0;s<Nd;s++){
|
||||
for(int c=0;c<Nc;c++){
|
||||
|
||||
LatticeFermion src4 (UGrid);
|
||||
PropToFerm<Action>(src4,seqsrc,s,c);
|
||||
|
||||
LatticeFermion result4(UGrid); result4=Zero();
|
||||
schur(Dw,src4,result4,zpg);
|
||||
|
||||
FermToProp<Action>(seqprop,result4,s,c);
|
||||
}
|
||||
}
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
Gamma gT(Gamma::Algebra::GammaT);
|
||||
|
||||
std::vector<TComplex> sumSV;
|
||||
std::vector<TComplex> sumVV;
|
||||
|
||||
Dw.ContractConservedCurrent(prop4,prop4,Vector_mu,phys_src,Current::Vector,Tdir);
|
||||
|
||||
SV = trace(Vector_mu); // Scalar-Vector conserved current
|
||||
VV = trace(gT*Vector_mu); // (local) Vector-Vector conserved current
|
||||
|
||||
// Spatial sum
|
||||
sliceSum(SV,sumSV,Tdir);
|
||||
sliceSum(VV,sumVV,Tdir);
|
||||
|
||||
const int Nt{static_cast<int>(sumSV.size())};
|
||||
|
||||
std::cout<<GridLogMessage<<"Vector Ward identity by timeslice (~ 0)"<<std::endl;
|
||||
for(int t=0;t<Nt;t++){
|
||||
std::cout<<GridLogMessage <<" t "<<t<<" SV "<<real(TensorRemove(sumSV[t]))<<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
|
||||
assert(abs(real(TensorRemove(sumSV[t]))) < 1e-10);
|
||||
assert(abs(real(TensorRemove(sumVV[t]))) < 1e-2);
|
||||
}
|
||||
|
||||
///////////////////////////////
|
||||
// 3pt vs 2pt check
|
||||
///////////////////////////////
|
||||
{
|
||||
Gamma::Algebra gA = Gamma::Algebra::Identity;
|
||||
Gamma g(gA);
|
||||
|
||||
LatticePropagator cur(UGrid);
|
||||
LatticePropagator tmp(UGrid);
|
||||
LatticeComplex c(UGrid);
|
||||
SpinColourMatrix qSite;
|
||||
peekSite(qSite, seqprop, coor);
|
||||
|
||||
Complex test_S, test_V, check_S, check_V;
|
||||
|
||||
std::vector<TComplex> check_buf;
|
||||
|
||||
test_S = trace(qSite*g);
|
||||
test_V = trace(qSite*g*Gamma::gmu[mu_J]);
|
||||
|
||||
Dw.ContractConservedCurrent(prop4,prop4,cur,phys_src,curr,mu_J);
|
||||
|
||||
c = trace(cur*g);
|
||||
sliceSum(c, check_buf, Tp);
|
||||
check_S = TensorRemove(check_buf[t_J]);
|
||||
|
||||
auto gmu=Gamma::gmu[mu_J];
|
||||
c = trace(cur*g*gmu);
|
||||
sliceSum(c, check_buf, Tp);
|
||||
check_V = TensorRemove(check_buf[t_J]);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage << std::setprecision(14)<<"Test S = " << abs(test_S) << std::endl;
|
||||
std::cout<<GridLogMessage << "Test V = " << abs(test_V) << std::endl;
|
||||
std::cout<<GridLogMessage << "Check S = " << abs(check_S) << std::endl;
|
||||
std::cout<<GridLogMessage << "Check V = " << abs(check_V) << std::endl;
|
||||
|
||||
// Check difference = 0
|
||||
check_S = check_S - test_S;
|
||||
check_V = check_V - test_V;
|
||||
|
||||
std::cout<<GridLogMessage << "Consistency check for sequential conserved " <<std::endl;
|
||||
std::cout<<GridLogMessage << "Diff S = " << abs(check_S) << std::endl;
|
||||
assert(abs(check_S) < 1e-8);
|
||||
std::cout<<GridLogMessage << "Diff V = " << abs(check_V) << std::endl;
|
||||
assert(abs(check_V) < 1e-8);
|
||||
}
|
||||
|
||||
}
|
@ -132,8 +132,8 @@ int main(int argc, char **argv) {
|
||||
|
||||
// Checkpointer definition
|
||||
CheckpointerParameters CPparams(Reader);
|
||||
//TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
|
||||
TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar);
|
||||
TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
|
||||
//TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar); this breaks for compilation without lime
|
||||
|
||||
RNGModuleParameters RNGpar(Reader);
|
||||
TheHMC.Resources.SetRNGSeeds(RNGpar);
|
||||
|
@ -74,10 +74,10 @@ int main(int argc, char **argv) {
|
||||
|
||||
// Checkpointer definition
|
||||
CheckpointerParameters CPparams(Reader);
|
||||
//TheHMC.Resources.LoadNerscCheckpointer(CPparams);
|
||||
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
|
||||
|
||||
// Store metadata in the Scidac checkpointer
|
||||
TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar);
|
||||
// Store metadata in the Scidac checkpointer - obviously breaks without LIME
|
||||
//TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar);
|
||||
|
||||
RNGModuleParameters RNGpar(Reader);
|
||||
TheHMC.Resources.SetRNGSeeds(RNGpar);
|
||||
|
@ -37,6 +37,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
#ifdef HAVE_LIME
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis>
|
||||
class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos<Fobj,CComplex,nbasis>
|
||||
{
|
||||
@ -249,3 +251,11 @@ int main (int argc, char ** argv) {
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
#else
|
||||
|
||||
int main( void )
|
||||
{
|
||||
return 0 ;
|
||||
}
|
||||
|
||||
#endif // HAVE_LIME_H
|
||||
|
@ -36,7 +36,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
#ifdef HAVE_LIME
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis>
|
||||
class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos<Fobj,CComplex,nbasis>
|
||||
@ -250,3 +251,11 @@ int main (int argc, char ** argv) {
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
#else
|
||||
|
||||
int main( void )
|
||||
{
|
||||
return 0 ;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
@ -93,8 +93,16 @@ int main(int argc, char** argv) {
|
||||
// Setup of Dirac Matrix and Operator //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu);
|
||||
|
||||
LatticeGaugeField Umu(Grid_f);
|
||||
#if (Nc==2)
|
||||
SU2::HotConfiguration(pRNG_f, Umu);
|
||||
#elif (defined Nc==3)
|
||||
SU3::HotConfiguration(pRNG_f, Umu);
|
||||
#elif (defined Nc==4)
|
||||
SU4::HotConfiguration(pRNG_f, Umu);
|
||||
#elif (defined Nc==5)
|
||||
SU5::HotConfiguration(pRNG_f, Umu);
|
||||
#endif
|
||||
RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15;
|
||||
|
||||
RealD mass = -0.30;
|
||||
|
@ -34,6 +34,7 @@ using namespace Grid;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
#ifdef HAVE_LIME
|
||||
typedef typename DomainWallFermionR::FermionField FermionField;
|
||||
typedef typename DomainWallFermionR::ComplexField ComplexField;
|
||||
typename DomainWallFermionR::ImplParams params;
|
||||
@ -237,4 +238,5 @@ int main (int argc, char ** argv)
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
#endif // HAVE_LIME
|
||||
}
|
||||
|
Loading…
x
Reference in New Issue
Block a user