portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-03 21:44:33 +00:00
Code Releases Activity

Clean up

Browse Source
This commit is contained in:
Peter Boyle
2023-05-19 21:20:45 -04:00
parent 9955bf9daf
commit 29a4bfe5e5
1 changed files with 118 additions and 260 deletions
Download Patch File Download Diff File Expand all files Collapse all files

370
Grid/qcd/smearing/GaugeConfigurationMasked.h
View File

@@ -6,19 +6,12 @@
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
template<class T> void Dump(Lattice<T> & lat, std::string s,Coordinate site = Coordinate({0,0,0,0}))
{
typename T::scalar_object tmp;
peekSite(tmp,lat,site);
std::cout << " GRID "<<s<<" "<<tmp<<std::endl;
}
/*! /*!
@brief Smeared configuration masked container @brief Smeared configuration masked container
Modified for a multi-subset smearing (aka Luscher Flowed HMC) Modified for a multi-subset smearing (aka Luscher Flowed HMC)
*/ */
template <class Gimpl> template <class Gimpl>
class SmearedConfigurationMasked class SmearedConfigurationMasked : public SmearedConfiguration<Gimpl>
{ {
public: public:
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
@@ -29,33 +22,28 @@ private:
std::vector<GaugeField> SmearedSet; std::vector<GaugeField> SmearedSet;
std::vector<LatticeLorentzComplex> masks; std::vector<LatticeLorentzComplex> masks;
///////////////////////////////////////////////////
// Additions
// Could just set masks to 1 in case of oldsmearing
// and share the code. Pass in a mask creation object with alternate constructor -- SmearMaskMaker
///////////////////////////////////////////////////
typedef typename SU3Adjoint::AMatrix AdjMatrix; typedef typename SU3Adjoint::AMatrix AdjMatrix;
typedef typename SU3Adjoint::LatticeAdjMatrix AdjMatrixField; typedef typename SU3Adjoint::LatticeAdjMatrix AdjMatrixField;
typedef typename SU3Adjoint::LatticeAdjVector AdjVectorField; typedef typename SU3Adjoint::LatticeAdjVector AdjVectorField;
/*
set_zero(f_det_basis); // Adjoint vector to GaugeField force
for (int a = 0; a < 8; ++a) { void InsertForce(GaugeField &Fdet,AdjVectorField &Fdet_nu,int nu)
const ColorMatrix uc = uc_pre * ts[a] * uc_post; {
for (int c = 0; c < 8; ++c) { Complex ci(0,1);
const ColorMatrix d_n = make_tr_less_anti_herm_matrix(ts[c] * uc); GaugeLinkField Fdet_pol(Fdet.Grid());
const array<double, 8> d_n_b = Fdet_pol=Zero();
basis_projection_anti_hermitian_matrix(d_n); for(int e=0;e<8;e++){
for (int b = 0; b < 8; ++b) { ColourMatrix te;
f_det_basis[a] += n_e_mp_inv_j_x_mat(c, b) * d_n_b[b]; SU3::generatorQlat(e, te);
auto tmp=peekColour(Fdet_nu,e);
Fdet_pol=Fdet_pol + ci*tmp*te; // but norm of te is different.. why?
} }
pokeLorentz(Fdet, Fdet_pol, nu);
} }
} void Compute_MpInvJx_dNxxdSy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
*/
void AdjointDeriv2(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
{ {
GaugeLinkField UtaU(PlaqL.Grid()); GaugeLinkField UtaU(PlaqL.Grid());
GaugeLinkField D(PlaqL.Grid()); GaugeLinkField D(PlaqL.Grid());
GaugeLinkField aa(PlaqL.Grid());
AdjMatrixField Dbc(PlaqL.Grid()); AdjMatrixField Dbc(PlaqL.Grid());
LatticeComplex tmp(PlaqL.Grid()); LatticeComplex tmp(PlaqL.Grid());
const int Ngen = SU3Adjoint::Dimension; const int Ngen = SU3Adjoint::Dimension;
@@ -69,24 +57,18 @@ private:
for(int c=0;c<Ngen;c++) { for(int c=0;c<Ngen;c++) {
SU3::generatorQlat(c, tc); SU3::generatorQlat(c, tc);
D = Ta( (2.0)*ci*tc *UtaU); D = Ta( (2.0)*ci*tc *UtaU);
// Dump(D," (-4.0)*Ta( tc *UtaU)");
for(int b=0;b<Ngen;b++){ for(int b=0;b<Ngen;b++){
SU3::generatorQlat(b, tb); SU3::generatorQlat(b, tb);
tmp =-trace(ci*tb*D); tmp =-trace(ci*tb*D);
PokeIndex<ColourIndex>(Dbc,tmp,b,c); // Adjoint rep PokeIndex<ColourIndex>(Dbc,tmp,b,c); // Adjoint rep
// Dump(tmp," -trace(ci*tb*D)");
} }
} }
// Dump(Dbc," Dbc ");
// Dump(MpInvJx," MpInvJx ");
tmp = trace(MpInvJx * Dbc); tmp = trace(MpInvJx * Dbc);
// Dump(tmp," trace(MpInvJx * Dbc) ");
PokeIndex<ColourIndex>(Fdet2,tmp,a); PokeIndex<ColourIndex>(Fdet2,tmp,a);
} }
Dump(Fdet2," Fdet2 ");
} }
void AdjointDeriv(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR,AdjMatrixField &NxAd) void ComputeNxy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR,AdjMatrixField &NxAd)
{ {
GaugeLinkField Nx(PlaqL.Grid()); GaugeLinkField Nx(PlaqL.Grid());
const int Ngen = SU3Adjoint::Dimension; const int Ngen = SU3Adjoint::Dimension;
@@ -95,12 +77,11 @@ private:
ColourMatrix tc; ColourMatrix tc;
for(int b=0;b<Ngen;b++) { for(int b=0;b<Ngen;b++) {
SU3::generatorQlat(b, tb); SU3::generatorQlat(b, tb);
// Qlat Tb = 2i Tb^Grid
Nx = (2.0)*Ta( adj(PlaqL)*ci*tb * PlaqR ); Nx = (2.0)*Ta( adj(PlaqL)*ci*tb * PlaqR );
for(int c=0;c<Ngen;c++) { for(int c=0;c<Ngen;c++) {
SU3::generatorQlat(c, tc); SU3::generatorQlat(c, tc);
auto tmp =closure( -trace(ci*tc*Nx)); // Luchang's norm: (2Tc) (2Td) N^db = -2 delta cd N^db // - was important auto tmp =closure( -trace(ci*tc*Nx));
PokeIndex<ColourIndex>(NxAd,tmp,c,b); // Adjoint rep PokeIndex<ColourIndex>(NxAd,tmp,c,b);
} }
} }
} }
@@ -116,6 +97,7 @@ private:
} }
} }
public: public:
void logDetJacobianForceLevel(const GaugeField &U, GaugeField &force ,int smr) void logDetJacobianForceLevel(const GaugeField &U, GaugeField &force ,int smr)
{ {
GridBase* grid = U.Grid(); GridBase* grid = U.Grid();
@@ -128,7 +110,7 @@ public:
std::vector<GaugeLinkField> Umu(Nd,grid); std::vector<GaugeLinkField> Umu(Nd,grid);
GaugeLinkField Cmu(grid); // U and staple; C contains factor of epsilon GaugeLinkField Cmu(grid); // U and staple; C contains factor of epsilon
GaugeLinkField Zx(grid); // U times Staple, contains factor of epsilon GaugeLinkField Zx(grid); // U times Staple, contains factor of epsilon
GaugeLinkField Nx(grid); // Nxx fundamental space GaugeLinkField Nxx(grid); // Nxx fundamental space
GaugeLinkField Utmp(grid); GaugeLinkField Utmp(grid);
GaugeLinkField PlaqL(grid); GaugeLinkField PlaqL(grid);
GaugeLinkField PlaqR(grid); GaugeLinkField PlaqR(grid);
@@ -136,11 +118,12 @@ public:
AdjMatrix TRb; AdjMatrix TRb;
ColourMatrix Ident; ColourMatrix Ident;
LatticeComplex cplx(grid); LatticeComplex cplx(grid);
AdjVectorField AlgV(grid);
AdjVectorField AlgVtmp(grid); AdjVectorField dJdXe_nMpInv(grid);
AdjVectorField dJdXe_nMpInv_y(grid);
AdjMatrixField MpAd(grid); // Mprime luchang's notes AdjMatrixField MpAd(grid); // Mprime luchang's notes
AdjMatrixField MpAdInv(grid); // Mprime inverse AdjMatrixField MpAdInv(grid); // Mprime inverse
AdjMatrixField NxAd(grid); // Nxx in adjoint space AdjMatrixField NxxAd(grid); // Nxx in adjoint space
AdjMatrixField JxAd(grid); AdjMatrixField JxAd(grid);
AdjMatrixField ZxAd(grid); AdjMatrixField ZxAd(grid);
AdjMatrixField mZxAd(grid); AdjMatrixField mZxAd(grid);
@@ -153,42 +136,49 @@ public:
} }
int mu= (smr/2) %Nd; int mu= (smr/2) %Nd;
////////////////////////////////////////////////////////////////////////////////
// Mask the gauge field
////////////////////////////////////////////////////////////////////////////////
auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask
// Mask the gauge field
Umsk = U; Umsk = U;
ApplyMask(Umsk,smr); ApplyMask(Umsk,smr);
Utmp = peekLorentz(Umsk,mu); Utmp = peekLorentz(Umsk,mu);
////////////////////////////////////////////////////////////////////////////////
// Retrieve the eps/rho parameter(s) -- could allow all different but not so far
////////////////////////////////////////////////////////////////////////////////
double rho=StoutSmearing->SmearRho[1];
int idx=0;
for(int mu=0;mu<4;mu++){
for(int nu=0;nu<4;nu++){
if ( mu!=nu) assert(StoutSmearing->SmearRho[idx]==rho);
else assert(StoutSmearing->SmearRho[idx]==0.0);
idx++;
}}
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// Assemble the N matrix // Assemble the N matrix
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// Computes ALL the staples -- could compute one only here // Computes ALL the staples -- could compute one only and do it here
StoutSmearing->BaseSmear(C, U); StoutSmearing->BaseSmear(C, U);
double rho=0.1;
Cmu = peekLorentz(C, mu); Cmu = peekLorentz(C, mu);
Dump(Cmu,std::string(" Cmu "));
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// Assemble Luscher exp diff map J matrix // Assemble Luscher exp diff map J matrix
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// Ta so Z lives in Lie algabra // Ta so Z lives in Lie algabra
Zx = Ta(Cmu * adj(Umu[mu])); Zx = Ta(Cmu * adj(Umu[mu]));
// Dump(Zx,std::string("Zx"));
// Move Z to the Adjoint Rep == make_adjoint_representation // Move Z to the Adjoint Rep == make_adjoint_representation
ZxAd = Zero(); ZxAd = Zero();
for(int b=0;b<8;b++) { for(int b=0;b<8;b++) {
// Adj group sets traceless antihermitian T's -- Guido, really???? // Adj group sets traceless antihermitian T's -- Guido, really????
// Is the mapping of these the same? Same structure constants
// Might never have been checked.
SU3::generatorQlat(b, tb); // Fund group sets traceless hermitian T's SU3::generatorQlat(b, tb); // Fund group sets traceless hermitian T's
SU3Adjoint::generatorQlat(b,TRb); SU3Adjoint::generatorQlat(b,TRb);
TRb=-TRb; TRb=-TRb;
cplx = 2.0*trace(ci*tb*Zx); // my convention 1/2 delta ba cplx = 2.0*trace(ci*tb*Zx); // my convention 1/2 delta ba
ZxAd = ZxAd + cplx * TRb; // is this right? YES - Guido used Anti herm Ta's and with bloody wrong sign. ZxAd = ZxAd + cplx * TRb; // is this right? YES - Guido used Anti herm Ta's and with bloody wrong sign.
} }
Dump(ZxAd,std::string("ZxAd"));
////////////////////////////////////// //////////////////////////////////////
// J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)! // J(x) = 1 + Sum_k=1..N (-Zac)^k/(k+1)!
@@ -202,17 +192,9 @@ public:
kpfac = kpfac /(k+1); kpfac = kpfac /(k+1);
JxAd = JxAd + X * kpfac; JxAd = JxAd + X * kpfac;
} }
Dump(JxAd,std::string("JxAd"));
////////////////////////////////////// //////////////////////////////////////
// dJ(x)/dxe = d/dxe Sum_k=1..N X^k/(k+1)! // dJ(x)/dxe
// = 1/2! te
// + 1/3! (te x + x te ) )
// + 1/4! (te x x + x te x + x x te )
// + 1/5! (te x x x + x te x x + x x te x + x x x te )
// Iterate:
// teX_n = teX_{n-1} x
// S_n = x S_{n-1} + teX_{n}
////////////////////////////////////// //////////////////////////////////////
std::vector<AdjMatrixField> dJdX; dJdX.resize(8,grid); std::vector<AdjMatrixField> dJdX; dJdX.resize(8,grid);
AdjMatrixField tbXn(grid); AdjMatrixField tbXn(grid);
@@ -225,7 +207,7 @@ public:
for(int b=0;b<8;b++){ for(int b=0;b<8;b++){
aunit = ComplexD(1.0); aunit = ComplexD(1.0);
SU3Adjoint::generatorQlat(b, TRb); //dt2 SU3Adjoint::generatorQlat(b, TRb); //dt2
Dump(ZxAd,std::string("ZxAd"));
X = (-1.0)*ZxAd; X = (-1.0)*ZxAd;
t2 = X; t2 = X;
dt2 = TRb; dt2 = TRb;
@@ -235,212 +217,131 @@ public:
t2 = X * t3; t2 = X * t3;
dt2 = TRb * t3 + X * dt3; dt2 = TRb * t3 + X * dt3;
} }
// dt3 = .5 * dt2; dJdX[b] = -dt2;
dJdX[b] = -dt2; // sign and 2x ?
Dump(dJdX[b],std::string("dJdX"));
} }
/*
X = (-1.0)*ZxAd; //x=t2
// n=1 point
tbXn=TRb;
sumXtbX= TRb;
RealD kpfac = 1.0/2.0;
dJdX[b] = sumXtbX *kpfac;
for(int k=2;k<12;k++){
kpfac = kpfac /(k+1);
tbXn = tbXn*X;
sumXtbX = X*sumXtbX + tbXn;
dJdX[b] = dJdX[b] + sumXtbX *kpfac;
}
*/
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Mask Umu for this link // Mask Umu for this link
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Nx = (2.0)*Ta( adj(PlaqL)*ci*tb * PlaqR );
PlaqL = Ident; PlaqL = Ident;
PlaqR = Utmp*adj(Cmu); PlaqR = Utmp*adj(Cmu);
AdjointDeriv(PlaqL,PlaqR,NxAd); ComputeNxy(PlaqL,PlaqR,NxxAd);
Dump(NxAd,std::string("NxAd"));
//////////////////////////// ////////////////////////////
// Mab // Mab
//////////////////////////// ////////////////////////////
// Mab = Complex(1.0,0.0);
// Mab = Mab - Jac * Ncb;
MpAd = Complex(1.0,0.0); MpAd = Complex(1.0,0.0);
MpAd = MpAd - JxAd * NxAd; MpAd = MpAd - JxAd * NxxAd;
Dump(MpAd,std::string("MpAd"));
///////////////////////// /////////////////////////
// invert the 8x8 // invert the 8x8
///////////////////////// /////////////////////////
MpAdInv = Inverse(MpAd); MpAdInv = Inverse(MpAd);
Dump(MpAdInv,std::string("MpAdInv"));
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Alternate way of creating // Nxx Mp^-1
// May need to keep the +nu and -nu plaq fields distinct to stop collision
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
AdjVectorField FdetV(grid);
AdjVectorField Fdet1_nu(grid);
AdjVectorField Fdet2_nu(grid);
AdjVectorField Fdet2_mu(grid);
AdjVectorField Fdet1_mu(grid);
AdjMatrixField nMpInv(grid); AdjMatrixField nMpInv(grid);
nMpInv= NxAd *MpAdInv; nMpInv= NxxAd *MpAdInv;
Dump(nMpInv," nMpInv ");
AdjMatrixField MpInvJx(grid); AdjMatrixField MpInvJx(grid);
AdjMatrixField MpInvJx_nu(grid); AdjMatrixField MpInvJx_nu(grid);
MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor
Dump(MpInvJx," MpInvJx ");
AdjVectorField FdetV(grid); Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
AdjVectorField FdetV2(grid); Fdet2_mu=FdetV;
AdjVectorField FdetV2_mu(grid); Fdet1_mu=Zero();
// First shot at the Fdet2
AdjointDeriv2(PlaqL,PlaqR,MpInvJx,FdetV);
FdetV2_mu=FdetV;
Dump(FdetV,std::string(" FdetV2xx_mu "));
for(int e =0 ; e<8 ; e++){ for(int e =0 ; e<8 ; e++){
LatticeComplexD tr(grid); LatticeComplexD tr(grid);
ColourMatrix te; ColourMatrix te;
SU3::generatorQlat(e, te); SU3::generatorQlat(e, te);
tr = trace(dJdX[e] * nMpInv); tr = trace(dJdX[e] * nMpInv);
pokeColour(AlgV,tr,e); pokeColour(dJdXe_nMpInv,tr,e);
// std::cout << " ***** tr " <<e<<std::endl;
// Dump(tr,std::string("tr"));
} }
///////////////////////////////
// Mask it off // Mask it off
///////////////////////////////
auto tmp=PeekIndex<LorentzIndex>(masks[smr],mu); auto tmp=PeekIndex<LorentzIndex>(masks[smr],mu);
AlgV = AlgV*tmp; dJdXe_nMpInv = dJdXe_nMpInv*tmp;
Dump(AlgV,std::string("AlgV"));
// AlgV needs to multiply: // dJdXe_nMpInv needs to multiply:
// NxAd (site local) (1) // Nxx_mu (site local) (1)
// NfmuPlus (site local) one site forward in each nu direction (3) // Nxy_mu one site forward in each nu direction (3)
// NfmuMinus (site local) one site backward in each nu direction(3) // Nxy_mu one site backward in each nu direction (3)
// Nfnu (site local) 0,0 ; +mu,0; 0,-nu; +mu-nu [ 3x4 = 12] // Nxy_nu 0,0 ; +mu,0; 0,-nu; +mu-nu [ 3x4 = 12]
// 19 terms. // 19 terms.
AdjVectorField AlgVmu_p(grid); AlgVmu_p=Zero(); AdjMatrixField Nxy(grid);
AdjVectorField AlgVmu_m(grid); AlgVmu_m=Zero();
AdjVectorField AlgVnu(grid); AlgVnu=Zero();
GaugeLinkField FmuPlus(grid);
GaugeLinkField FmuMinus(grid);
GaugeLinkField Fnu(grid);
GaugeLinkField Fnumu(grid);
std::vector<AdjMatrixField> Nfperp(Nd,grid); // Why needed vs nu?
AdjMatrixField NfmuPlus(grid);
AdjMatrixField NfmuMinus(grid);
for(int d=0;d<Nd;d++){
Nfperp[d] = Zero();
}
NfmuPlus=Zero();
NfmuMinus=Zero();
//////////////////////////////////////////////////////
// six force inserts x 3 directions from OTHER links y!=x || mu!=nu
//
// To avoid collison of x-y pairs need to keep +ve nu and -ve nu separate, at least for mu force
// FmuPlus, FmuMinus and Fnu/Fperp
//////////////////////////////////////////////////////
GaugeField Fdet1(grid); GaugeField Fdet1(grid);
GaugeField Fdet2(grid); GaugeField Fdet2(grid);
GaugeLinkField Fdet1_pol(grid); GaugeLinkField Fdet_pol(grid); // one polarisation
GaugeLinkField Fdet2_pol(grid);
AdjVectorField FdetV_acc(grid);
AdjVectorField FdetV2_acc(grid);
AdjVectorField FdetV_mu(grid); FdetV_mu=Zero();
for(int nu=0;nu<Nd;nu++){ for(int nu=0;nu<Nd;nu++){
if (nu!=mu) { if (nu!=mu) {
LatticeComplexD tr(grid);
int Lnu = grid->GlobalDimensions()[nu];
Coordinate coormu({0,0,0,0}); coormu[mu]=1;
Coordinate coornu({0,0,0,0}); coornu[nu]=1;
Coordinate coornnu({0,0,0,0}); coornnu[nu]=Lnu-1;
Coordinate coormunu({0,0,0,0}); coormunu[mu]=1; coormunu[nu]=1;
Coordinate coormunnu({0,0,0,0}); coormunnu[mu]=1; coormunnu[nu]=Lnu-1;
///////////////// +ve nu ///////////////// ///////////////// +ve nu /////////////////
// __ // __
// | | // | |
// x== // nu polarisation -- clockwise // x== // nu polarisation -- clockwise
std::cout << " ********************* "<<std::endl;
std::cout << " nu+ = "<<nu<<std::endl;
std::cout << " ********************* "<<std::endl;
PlaqL=Ident; PlaqL=Ident;
PlaqR=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu, PlaqR=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu,
Gimpl::CovShiftForward(Umu[mu], mu, Gimpl::CovShiftForward(Umu[mu], mu,
Gimpl::CovShiftBackward(Umu[nu], nu, Gimpl::CovShiftBackward(Umu[nu], nu,
Gimpl::CovShiftIdentityBackward(Utmp, mu)))); Gimpl::CovShiftIdentityBackward(Utmp, mu))));
// Nx = -2.0*ci*trace(tc*Ta( adj(PlaqL)*ci*tb * PlaqR ))); dJdXe_nMpInv_y = dJdXe_nMpInv;
AdjointDeriv(PlaqL,PlaqR,Nfperp[nu]); ComputeNxy(PlaqL,PlaqR,Nxy);
Fdet1_nu = transpose(Nxy)*dJdXe_nMpInv_y;
AlgVnu = AlgV;
PlaqR=(-1.0)*PlaqR; PlaqR=(-1.0)*PlaqR;
FdetV = transpose(Nfperp[nu])*AlgVnu; Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
Dump(FdetV,std::string("FdetVxy_nu ; y=x ")); // OK Fdet2_nu = FdetV;
FdetV_acc = FdetV;
AdjointDeriv2(PlaqL,PlaqR,MpInvJx,FdetV2);
FdetV2_acc = FdetV2;
Dump(FdetV2,std::string("FdetV2xy_nu ; y=x "));
// x== // x==
// | | // | |
// .__| // nu polarisation -- anticlockwise // .__| // nu polarisation -- anticlockwise
// Nx = (2.0)*Ta( adj(PlaqL)*ci*tb * PlaqR );
PlaqR=(rho)*Gimpl::CovShiftForward(Umu[nu], nu, PlaqR=(rho)*Gimpl::CovShiftForward(Umu[nu], nu,
Gimpl::CovShiftBackward(Umu[mu], mu, Gimpl::CovShiftBackward(Umu[mu], mu,
Gimpl::CovShiftIdentityBackward(Umu[nu], nu))); Gimpl::CovShiftIdentityBackward(Umu[nu], nu)));
PlaqL=Gimpl::CovShiftIdentityBackward(Utmp, mu); PlaqL=Gimpl::CovShiftIdentityBackward(Utmp, mu);
AdjointDeriv(PlaqL, PlaqR,Nfperp[nu]); dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
ComputeNxy(PlaqL, PlaqR,Nxy);
Fdet1_nu = Fdet1_nu+transpose(Nxy)*dJdXe_nMpInv_y;
AlgVnu = Cshift(AlgV,mu,-1);
FdetV = transpose(Nfperp[nu])*AlgVnu;
Dump(FdetV,std::string("FdetVxy_nu ; +mu "),coormu); // OK
FdetV_acc = FdetV_acc + FdetV;
MpInvJx_nu = Cshift(MpInvJx,mu,-1); MpInvJx_nu = Cshift(MpInvJx,mu,-1);
AdjointDeriv2(PlaqL,PlaqR,MpInvJx_nu,FdetV2); Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
FdetV2_acc = FdetV2_acc+FdetV2; Fdet2_nu = Fdet2_nu+FdetV;
Dump(FdetV2,std::string("FdetV2xy_nu ; +mu "),coormu);
///////////////// -ve nu ///////////////// ///////////////// -ve nu /////////////////
// __ // __
// | | // | |
// x== // nu polarisation -- clockwise // x== // nu polarisation -- clockwise
std::cout << " ********************* "<<std::endl;
std::cout << " nu- = "<<nu<<std::endl;
std::cout << " ********************* "<<std::endl;
PlaqL=(rho)* Gimpl::CovShiftForward(Umu[mu], mu, PlaqL=(rho)* Gimpl::CovShiftForward(Umu[mu], mu,
Gimpl::CovShiftForward(Umu[nu], nu, Gimpl::CovShiftForward(Umu[nu], nu,
Gimpl::CovShiftIdentityBackward(Utmp, mu))); Gimpl::CovShiftIdentityBackward(Utmp, mu)));
PlaqR = Gimpl::CovShiftIdentityForward(Umu[nu], nu); PlaqR = Gimpl::CovShiftIdentityForward(Umu[nu], nu);
AdjointDeriv(PlaqL,PlaqR,Nfperp[nu]);
AlgVnu = Cshift(AlgV,nu,1); dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
ComputeNxy(PlaqL,PlaqR,Nxy);
FdetV = transpose(Nfperp[nu])*AlgVnu; Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
Dump(FdetV,std::string("FdetVxy_nu ; -nu"),coornnu);
FdetV_acc = FdetV_acc + FdetV;
MpInvJx_nu = Cshift(MpInvJx,nu,1); MpInvJx_nu = Cshift(MpInvJx,nu,1);
AdjointDeriv2(PlaqL,PlaqR,MpInvJx_nu,FdetV2); Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
FdetV2_acc = FdetV2_acc+FdetV2; Fdet2_nu = Fdet2_nu+FdetV;
Dump(FdetV2,std::string("FdetV2xy_nu ; -nu "),coornnu);
// x== // x==
// | | // | |
@@ -448,43 +349,26 @@ public:
PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu, PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[nu], nu,
Gimpl::CovShiftIdentityBackward(Utmp, mu)); Gimpl::CovShiftIdentityBackward(Utmp, mu));
PlaqR=Gimpl::CovShiftBackward(Umu[mu], mu, PlaqR=Gimpl::CovShiftBackward(Umu[mu], mu,
Gimpl::CovShiftIdentityForward(Umu[nu], nu)); Gimpl::CovShiftIdentityForward(Umu[nu], nu));
AdjointDeriv(PlaqL,PlaqR,Nfperp[nu]); dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,mu,-1);
dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv_y,nu,1);
AlgVnu = Cshift(AlgV,mu,-1); ComputeNxy(PlaqL,PlaqR,Nxy);
AlgVnu = Cshift(AlgVnu,nu,1); Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
FdetV = transpose(Nfperp[nu])*AlgVnu;
Dump(FdetV,std::string("FdetVxy_nu; -nu +mu"),coormunnu);
FdetV_acc = FdetV_acc + FdetV;
MpInvJx_nu = Cshift(MpInvJx,mu,-1); MpInvJx_nu = Cshift(MpInvJx,mu,-1);
MpInvJx_nu = Cshift(MpInvJx_nu,nu,1); MpInvJx_nu = Cshift(MpInvJx_nu,nu,1);
AdjointDeriv2(PlaqL,PlaqR,MpInvJx_nu,FdetV2); Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
FdetV2_acc = FdetV2_acc+FdetV2; Fdet2_nu = Fdet2_nu+FdetV;
Dump(FdetV2,std::string("FdetV2xy_nu ; -nu +mu "),coormunu);
/////////////////////////////////////////////////////////////////////
// Set up the determinant force contribution in 3x3 algebra basis // Set up the determinant force contribution in 3x3 algebra basis
Fdet1_pol=Zero(); /////////////////////////////////////////////////////////////////////
for(int e=0;e<8;e++){ InsertForce(Fdet1,Fdet1_nu,nu);
ColourMatrix te; InsertForce(Fdet2,Fdet2_nu,nu);
SU3::generatorQlat(e, te);
auto tmp=peekColour(FdetV_acc,e);
Fdet1_pol=Fdet1_pol + ci*tmp*te; // but norm of te is different.. why?
}
pokeLorentz(Fdet1, Fdet1_pol, nu);
Fdet2_pol=Zero();
for(int e=0;e<8;e++){
ColourMatrix te;
SU3::generatorQlat(e, te);
auto tmp=peekColour(FdetV2_acc,e);
Fdet2_pol=Fdet2_pol + ci*tmp*te; // but norm of te is different.. why?
}
pokeLorentz(Fdet2, Fdet2_pol, nu);
////////////////////////////////////////////////// //////////////////////////////////////////////////
// Parallel direction terms // Parallel direction terms
@@ -496,21 +380,18 @@ public:
PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu, PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu,
Gimpl::CovShiftBackward(Umu[nu], nu, Gimpl::CovShiftBackward(Umu[nu], nu,
Gimpl::CovShiftIdentityBackward(Utmp, mu))); Gimpl::CovShiftIdentityBackward(Utmp, mu)));
PlaqR=Gimpl::CovShiftIdentityBackward(Umu[nu], nu); PlaqR=Gimpl::CovShiftIdentityBackward(Umu[nu], nu);
AdjointDeriv(PlaqL,PlaqR,NfmuPlus); dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,-1);
AlgVnu = Cshift(AlgV,nu,-1); ComputeNxy(PlaqL,PlaqR,Nxy);
Fdet1_mu = Fdet1_mu + transpose(Nxy)*dJdXe_nMpInv_y;
FdetV = transpose(NfmuPlus)*AlgVnu;
Dump(FdetV,std::string("FdetV mu ; +nu "),coornu);
FdetV_mu = FdetV_mu + FdetV;
MpInvJx_nu = Cshift(MpInvJx,nu,-1); MpInvJx_nu = Cshift(MpInvJx,nu,-1);
AdjointDeriv2(PlaqL,PlaqR,MpInvJx_nu,FdetV2);
FdetV2_mu = FdetV2_mu+FdetV2; Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
Dump(FdetV2,std::string("FdetV2xy_mu ; +nu "),coornu); Fdet2_mu = Fdet2_mu+FdetV;
// __ // __
// " | // " |
@@ -519,50 +400,28 @@ public:
PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu, PlaqL=(-rho)*Gimpl::CovShiftForward(Umu[mu], mu,
Gimpl::CovShiftForward(Umu[nu], nu, Gimpl::CovShiftForward(Umu[nu], nu,
Gimpl::CovShiftIdentityBackward(Utmp, mu))); Gimpl::CovShiftIdentityBackward(Utmp, mu)));
PlaqR=Gimpl::CovShiftIdentityForward(Umu[nu], nu); PlaqR=Gimpl::CovShiftIdentityForward(Umu[nu], nu);
AdjointDeriv(PlaqL,PlaqR,NfmuMinus); dJdXe_nMpInv_y = Cshift(dJdXe_nMpInv,nu,1);
AlgVnu = Cshift(AlgV,nu,1); ComputeNxy(PlaqL,PlaqR,Nxy);
Fdet1_mu = Fdet1_mu + transpose(Nxy)*dJdXe_nMpInv_y;
FdetV = transpose(NfmuMinus)*AlgVnu;
Dump(FdetV,std::string("FdetV_xy mu ; -nu "),coornnu);
FdetV_mu = FdetV_mu + FdetV;
MpInvJx_nu = Cshift(MpInvJx,nu,1); MpInvJx_nu = Cshift(MpInvJx,nu,1);
AdjointDeriv2(PlaqL,PlaqR,MpInvJx_nu,FdetV2);
FdetV2_mu = FdetV2_mu+FdetV2; Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
Dump(FdetV2,std::string("FdetV2xy_mu ; -nu "),coornnu); Fdet2_mu = Fdet2_mu+FdetV;
} }
} }
FdetV = transpose(NxAd)*AlgV; Fdet1_mu = Fdet1_mu + transpose(NxxAd)*dJdXe_nMpInv;
Dump(FdetV,std::string(" FdetVxx_mu "));
FdetV_mu = FdetV_mu + FdetV;
Fdet1_pol=Zero();
for(int e=0;e<8;e++){
ColourMatrix te;
SU3::generatorQlat(e, te);
auto tmp=peekColour(FdetV_mu,e);
Fdet1_pol=Fdet1_pol + ci*tmp*te; // but norm of te is different.. why?
}
pokeLorentz(Fdet1, Fdet1_pol, mu);
Fdet2_pol=Zero(); InsertForce(Fdet1,Fdet1_mu,mu);
for(int e=0;e<8;e++){ InsertForce(Fdet2,Fdet2_mu,mu);
ColourMatrix te;
SU3::generatorQlat(e, te);
auto tmp=peekColour(FdetV2_mu,e);
Fdet2_pol=Fdet2_pol + ci*tmp*te; // but norm of te is different.. why?
}
pokeLorentz(Fdet2, Fdet2_pol, mu);
force = Fdet1 + Fdet2; force = Fdet1 + Fdet2;
// force = Fdet1;
// force = Fdet2;
} }
RealD logDetJacobianLevel(const GaugeField &U,int smr) RealD logDetJacobianLevel(const GaugeField &U,int smr)
@@ -776,7 +635,6 @@ private:
return SigmaK; return SigmaK;
} }
//////////////////////////////////////// ////////////////////////////////////////
// INHERIT THESE // INHERIT THESE
//////////////////////////////////////// ////////////////////////////////////////