mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Mobius force improvement
This commit is contained in:
paboyle
parent
90f6bc16bb
commit
78ef10e60f
@ -320,7 +320,7 @@ void CayleyFermion5D<Impl>::MDeriv (GaugeField &mat,const FermionField &U,const
|
|||||||
this->DhopDeriv(mat,U,Din,dag);
|
this->DhopDeriv(mat,U,Din,dag);
|
||||||
} else {
|
} else {
|
||||||
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
||||||
MeooeDag5D(U,Din);
|
Meooe5D(U,Din);
|
||||||
this->DhopDeriv(mat,Din,V,dag);
|
this->DhopDeriv(mat,Din,V,dag);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
@ -335,8 +335,8 @@ void CayleyFermion5D<Impl>::MoeDeriv(GaugeField &mat,const FermionField &U,const
|
|||||||
this->DhopDerivOE(mat,U,Din,dag);
|
this->DhopDerivOE(mat,U,Din,dag);
|
||||||
} else {
|
} else {
|
||||||
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
||||||
MeooeDag5D(U,Din);
|
Meooe5D(U,Din);
|
||||||
this->DhopDerivOE(mat,Din,V,dag);
|
this->DhopDerivOE(mat,Din,V,dag);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
template<class Impl>
|
template<class Impl>
|
||||||
@ -350,7 +350,7 @@ void CayleyFermion5D<Impl>::MeoDeriv(GaugeField &mat,const FermionField &U,const
|
|||||||
this->DhopDerivEO(mat,U,Din,dag);
|
this->DhopDerivEO(mat,U,Din,dag);
|
||||||
} else {
|
} else {
|
||||||
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
// U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
|
||||||
MeooeDag5D(U,Din);
|
Meooe5D(U,Din);
|
||||||
this->DhopDerivEO(mat,Din,V,dag);
|
this->DhopDerivEO(mat,Din,V,dag);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|||||||
287
tests/core/Test_mobius_even_odd.cc
Normal file
287
tests/core/Test_mobius_even_odd.cc
Normal file
@ -0,0 +1,287 @@
|
|||||||
|
/*************************************************************************************
|
||||||
|
|
||||||
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
|
Source file: ./tests/Test_dwf_even_odd.cc
|
||||||
|
|
||||||
|
Copyright (C) 2015
|
||||||
|
|
||||||
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||||
|
|
||||||
|
This program is free software; you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation; either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
This program is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License along
|
||||||
|
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||||
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||||
|
|
||||||
|
See the full license in the file "LICENSE" in the top level distribution directory
|
||||||
|
*************************************************************************************/
|
||||||
|
/* END LEGAL */
|
||||||
|
#include <Grid/Grid.h>
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
using namespace Grid;
|
||||||
|
using namespace Grid::QCD;
|
||||||
|
|
||||||
|
template<class d>
|
||||||
|
struct scal {
|
||||||
|
d internal;
|
||||||
|
};
|
||||||
|
|
||||||
|
Gamma::Algebra Gmu [] = {
|
||||||
|
Gamma::Algebra::GammaX,
|
||||||
|
Gamma::Algebra::GammaY,
|
||||||
|
Gamma::Algebra::GammaZ,
|
||||||
|
Gamma::Algebra::GammaT
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
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;
|
||||||
|
|
||||||
|
|
||||||
|
const int Ls=10;
|
||||||
|
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||||
|
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||||
|
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||||
|
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||||
|
|
||||||
|
std::vector<int> seeds4({1,2,3,4});
|
||||||
|
std::vector<int> seeds5({5,6,7,8});
|
||||||
|
|
||||||
|
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||||
|
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||||
|
|
||||||
|
LatticeFermion src (FGrid); random(RNG5,src);
|
||||||
|
LatticeFermion phi (FGrid); random(RNG5,phi);
|
||||||
|
LatticeFermion chi (FGrid); random(RNG5,chi);
|
||||||
|
LatticeFermion result(FGrid); result=zero;
|
||||||
|
LatticeFermion ref(FGrid); ref=zero;
|
||||||
|
LatticeFermion tmp(FGrid); tmp=zero;
|
||||||
|
LatticeFermion err(FGrid); tmp=zero;
|
||||||
|
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
|
||||||
|
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||||
|
|
||||||
|
// Only one non-zero (y)
|
||||||
|
Umu=zero;
|
||||||
|
for(int nn=0;nn<Nd;nn++){
|
||||||
|
random(RNG4,U[nn]);
|
||||||
|
if ( nn>0 )
|
||||||
|
U[nn]=zero;
|
||||||
|
PokeIndex<LorentzIndex>(Umu,U[nn],nn);
|
||||||
|
}
|
||||||
|
|
||||||
|
RealD mass=0.1;
|
||||||
|
RealD M5 =1.8;
|
||||||
|
std::vector < std::complex<double> > omegas;
|
||||||
|
#if 0
|
||||||
|
for(int i=0;i<Ls;i++){
|
||||||
|
double imag = 0.;
|
||||||
|
if (i==0) imag=1.;
|
||||||
|
if (i==Ls-1) imag=-1.;
|
||||||
|
std::complex<double> temp (0.25+0.01*i, imag*0.01);
|
||||||
|
omegas.push_back(temp);
|
||||||
|
}
|
||||||
|
#else
|
||||||
|
omegas.push_back( std::complex<double>(1.45806438985048,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(1.18231318389348,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.830951166685955,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.542352409156791,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.341985020453729,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.21137902619029,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.126074299502912,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
|
||||||
|
omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
|
||||||
|
omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
|
||||||
|
#endif
|
||||||
|
|
||||||
|
MobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5);
|
||||||
|
// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||||
|
|
||||||
|
LatticeFermion src_e (FrbGrid);
|
||||||
|
LatticeFermion src_o (FrbGrid);
|
||||||
|
LatticeFermion r_e (FrbGrid);
|
||||||
|
LatticeFermion r_o (FrbGrid);
|
||||||
|
LatticeFermion r_eo (FGrid);
|
||||||
|
LatticeFermion r_eeoo(FGrid);
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||||
|
|
||||||
|
pickCheckerboard(Even,src_e,src);
|
||||||
|
pickCheckerboard(Odd,src_o,src);
|
||||||
|
|
||||||
|
Ddwf.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
|
||||||
|
Ddwf.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
|
||||||
|
setCheckerboard(r_eo,r_o);
|
||||||
|
setCheckerboard(r_eo,r_e);
|
||||||
|
|
||||||
|
Ddwf.Mooee(src_e,r_e); std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
|
||||||
|
Ddwf.Mooee(src_o,r_o); std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
|
||||||
|
setCheckerboard(r_eeoo,r_e);
|
||||||
|
setCheckerboard(r_eeoo,r_o);
|
||||||
|
|
||||||
|
r_eo=r_eo+r_eeoo;
|
||||||
|
Ddwf.M(src,ref);
|
||||||
|
|
||||||
|
// std::cout<<GridLogMessage << r_eo<<std::endl;
|
||||||
|
// std::cout<<GridLogMessage << ref <<std::endl;
|
||||||
|
|
||||||
|
err= ref - r_eo;
|
||||||
|
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
|
||||||
|
|
||||||
|
LatticeComplex cerr(FGrid);
|
||||||
|
cerr = localInnerProduct(err,err);
|
||||||
|
// std::cout<<GridLogMessage << cerr<<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Test MooeeDagger is the dagger of Mooee by requiring "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
|
||||||
|
LatticeFermion chi_e (FrbGrid);
|
||||||
|
LatticeFermion chi_o (FrbGrid);
|
||||||
|
|
||||||
|
LatticeFermion dchi_e (FrbGrid);
|
||||||
|
LatticeFermion dchi_o (FrbGrid);
|
||||||
|
|
||||||
|
LatticeFermion phi_e (FrbGrid);
|
||||||
|
LatticeFermion phi_o (FrbGrid);
|
||||||
|
|
||||||
|
LatticeFermion dphi_e (FrbGrid);
|
||||||
|
LatticeFermion dphi_o (FrbGrid);
|
||||||
|
|
||||||
|
pickCheckerboard(Even,chi_e,chi);
|
||||||
|
pickCheckerboard(Odd ,chi_o,chi);
|
||||||
|
pickCheckerboard(Even,phi_e,phi);
|
||||||
|
pickCheckerboard(Odd ,phi_o,phi);
|
||||||
|
|
||||||
|
Ddwf.Mooee(chi_e,dchi_o);
|
||||||
|
Ddwf.Mooee(chi_o,dchi_e);
|
||||||
|
Ddwf.MooeeDag(phi_e,dphi_o);
|
||||||
|
Ddwf.MooeeDag(phi_o,dphi_e);
|
||||||
|
|
||||||
|
ComplexD pDce = innerProduct(phi_e,dchi_e);
|
||||||
|
ComplexD pDco = innerProduct(phi_o,dchi_o);
|
||||||
|
ComplexD cDpe = innerProduct(chi_e,dphi_e);
|
||||||
|
ComplexD cDpo = innerProduct(chi_o,dphi_o);
|
||||||
|
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
|
pickCheckerboard(Even,chi_e,chi);
|
||||||
|
pickCheckerboard(Odd ,chi_o,chi);
|
||||||
|
pickCheckerboard(Even,phi_e,phi);
|
||||||
|
pickCheckerboard(Odd ,phi_o,phi);
|
||||||
|
|
||||||
|
Ddwf.Meooe(chi_e,dchi_o);
|
||||||
|
Ddwf.Meooe(chi_o,dchi_e);
|
||||||
|
Ddwf.MeooeDag(phi_e,dphi_o);
|
||||||
|
Ddwf.MeooeDag(phi_o,dphi_e);
|
||||||
|
|
||||||
|
pDce = innerProduct(phi_e,dchi_e);
|
||||||
|
pDco = innerProduct(phi_o,dchi_o);
|
||||||
|
cDpe = innerProduct(chi_e,dphi_e);
|
||||||
|
cDpo = innerProduct(chi_o,dphi_o);
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
|
||||||
|
pickCheckerboard(Even,chi_e,chi);
|
||||||
|
pickCheckerboard(Odd ,chi_o,chi);
|
||||||
|
|
||||||
|
Ddwf.Mooee(chi_e,src_e);
|
||||||
|
Ddwf.MooeeInv(src_e,phi_e);
|
||||||
|
|
||||||
|
Ddwf.Mooee(chi_o,src_o);
|
||||||
|
Ddwf.MooeeInv(src_o,phi_o);
|
||||||
|
|
||||||
|
setCheckerboard(phi,phi_e);
|
||||||
|
setCheckerboard(phi,phi_o);
|
||||||
|
|
||||||
|
err = phi-chi;
|
||||||
|
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
|
||||||
|
pickCheckerboard(Even,chi_e,chi);
|
||||||
|
pickCheckerboard(Odd ,chi_o,chi);
|
||||||
|
|
||||||
|
Ddwf.MooeeDag(chi_e,src_e);
|
||||||
|
Ddwf.MooeeInvDag(src_e,phi_e);
|
||||||
|
|
||||||
|
Ddwf.MooeeDag(chi_o,src_o);
|
||||||
|
Ddwf.MooeeInvDag(src_o,phi_o);
|
||||||
|
|
||||||
|
setCheckerboard(phi,phi_e);
|
||||||
|
setCheckerboard(phi,phi_o);
|
||||||
|
|
||||||
|
err = phi-chi;
|
||||||
|
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
|
||||||
|
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
|
||||||
|
|
||||||
|
random(RNG5,phi);
|
||||||
|
random(RNG5,chi);
|
||||||
|
pickCheckerboard(Even,chi_e,chi);
|
||||||
|
pickCheckerboard(Odd ,chi_o,chi);
|
||||||
|
pickCheckerboard(Even,phi_e,phi);
|
||||||
|
pickCheckerboard(Odd ,phi_o,phi);
|
||||||
|
RealD t1,t2;
|
||||||
|
|
||||||
|
|
||||||
|
SchurDiagMooeeOperator<MobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||||
|
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
|
||||||
|
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
|
||||||
|
|
||||||
|
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
|
||||||
|
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
|
||||||
|
|
||||||
|
pDce = innerProduct(phi_e,dchi_e);
|
||||||
|
pDco = innerProduct(phi_o,dchi_o);
|
||||||
|
cDpe = innerProduct(chi_e,dphi_e);
|
||||||
|
cDpo = innerProduct(chi_o,dphi_o);
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
|
||||||
|
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
|
||||||
|
|
||||||
|
Grid_finalize();
|
||||||
|
}
|
||||||
153
tests/forces/Test_mobius_force.cc
Normal file
153
tests/forces/Test_mobius_force.cc
Normal file
@ -0,0 +1,153 @@
|
|||||||
|
/*************************************************************************************
|
||||||
|
|
||||||
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
|
Source file: ./tests/Test_dwf_force.cc
|
||||||
|
|
||||||
|
Copyright (C) 2015
|
||||||
|
|
||||||
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
|
||||||
|
This program is free software; you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation; either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
This program is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License along
|
||||||
|
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||||
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||||
|
|
||||||
|
See the full license in the file "LICENSE" in the top level distribution directory
|
||||||
|
*************************************************************************************/
|
||||||
|
/* END LEGAL */
|
||||||
|
#include <Grid/Grid.h>
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
using namespace Grid;
|
||||||
|
using namespace Grid::QCD;
|
||||||
|
|
||||||
|
int main (int argc, char ** argv)
|
||||||
|
{
|
||||||
|
Grid_init(&argc,&argv);
|
||||||
|
|
||||||
|
std::vector<int> latt_size = GridDefaultLatt();
|
||||||
|
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||||
|
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||||
|
|
||||||
|
const int Ls=8;
|
||||||
|
|
||||||
|
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||||
|
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||||
|
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||||
|
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||||
|
|
||||||
|
std::vector<int> seeds4({1,2,3,4});
|
||||||
|
std::vector<int> seeds5({5,6,7,8});
|
||||||
|
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||||
|
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||||
|
|
||||||
|
int threads = GridThread::GetThreads();
|
||||||
|
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
||||||
|
|
||||||
|
LatticeFermion phi (FGrid); gaussian(RNG5,phi);
|
||||||
|
LatticeFermion Mphi (FGrid);
|
||||||
|
LatticeFermion MphiPrime (FGrid);
|
||||||
|
|
||||||
|
LatticeGaugeField U(UGrid);
|
||||||
|
|
||||||
|
SU3::HotConfiguration(RNG4,U);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Unmodified matrix element
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD mass=0.01;
|
||||||
|
RealD M5=1.8;
|
||||||
|
RealD b=0.5;
|
||||||
|
RealD c=0.5;
|
||||||
|
MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
|
||||||
|
Ddwf.M (phi,Mphi);
|
||||||
|
|
||||||
|
ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p
|
||||||
|
|
||||||
|
// get the deriv of phidag MdagM phi with respect to "U"
|
||||||
|
LatticeGaugeField UdSdU(UGrid);
|
||||||
|
LatticeGaugeField tmp(UGrid);
|
||||||
|
|
||||||
|
|
||||||
|
Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
|
||||||
|
Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
|
||||||
|
|
||||||
|
LatticeFermion Ftmp (FGrid);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Modify the gauge field a little
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD dt = 0.0001;
|
||||||
|
|
||||||
|
LatticeColourMatrix mommu(UGrid);
|
||||||
|
LatticeColourMatrix forcemu(UGrid);
|
||||||
|
LatticeGaugeField mom(UGrid);
|
||||||
|
LatticeGaugeField Uprime(UGrid);
|
||||||
|
|
||||||
|
for(int mu=0;mu<Nd;mu++){
|
||||||
|
|
||||||
|
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
|
||||||
|
|
||||||
|
PokeIndex<LorentzIndex>(mom,mommu,mu);
|
||||||
|
|
||||||
|
// fourth order exponential approx
|
||||||
|
parallel_for(auto i=mom.begin();i<mom.end();i++){
|
||||||
|
Uprime[i](mu) =
|
||||||
|
U[i](mu)
|
||||||
|
+ mom[i](mu)*U[i](mu)*dt
|
||||||
|
+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
|
||||||
|
;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
Ddwf.ImportGauge(Uprime);
|
||||||
|
Ddwf.M (phi,MphiPrime);
|
||||||
|
|
||||||
|
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);
|
||||||
|
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
// Use derivative to estimate dS
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
|
||||||
|
|
||||||
|
LatticeComplex dS(UGrid); dS = zero;
|
||||||
|
for(int mu=0;mu<Nd;mu++){
|
||||||
|
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||||
|
mommu=Ta(mommu)*2.0;
|
||||||
|
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int mu=0;mu<Nd;mu++){
|
||||||
|
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
|
||||||
|
mommu = PeekIndex<LorentzIndex>(mom,mu);
|
||||||
|
|
||||||
|
// Update PF action density
|
||||||
|
dS = dS+trace(mommu*forcemu)*dt;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
Complex dSpred = sum(dS);
|
||||||
|
|
||||||
|
std::cout << GridLogMessage << " S "<<S<<std::endl;
|
||||||
|
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
|
||||||
|
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
|
||||||
|
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
|
||||||
|
|
||||||
|
std::cout<< GridLogMessage << "Done" <<std::endl;
|
||||||
|
Grid_finalize();
|
||||||
|
}
|
||||||
@ -317,6 +317,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
Grid::QCD::LatticeComplex rect(UGrid);
|
Grid::QCD::LatticeComplex rect(UGrid);
|
||||||
Grid::QCD::TComplex trect;
|
Grid::QCD::TComplex trect;
|
||||||
Grid::QCD::Complex crect;
|
Grid::QCD::Complex crect;
|
||||||
|
Grid::RealD rrect;
|
||||||
Grid::RealD vol = UGrid->gSites();
|
Grid::RealD vol = UGrid->gSites();
|
||||||
for(int mu=0;mu<Grid::QCD::Nd;mu++){
|
for(int mu=0;mu<Grid::QCD::Nd;mu++){
|
||||||
for(int nu=0;nu<Grid::QCD::Nd;nu++){
|
for(int nu=0;nu<Grid::QCD::Nd;nu++){
|
||||||
@ -325,7 +326,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
|
Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
|
||||||
trect = Grid::sum(rect);
|
trect = Grid::sum(rect);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
|
||||||
|
|
||||||
Grid::GridStopWatch Peter;
|
Grid::GridStopWatch Peter;
|
||||||
Grid::GridStopWatch Azusa;
|
Grid::GridStopWatch Azusa;
|
||||||
@ -355,7 +357,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
|
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect=real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
// __
|
// __
|
||||||
@ -370,7 +373,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
|
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect=real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
// __
|
// __
|
||||||
// |__ __ |
|
// |__ __ |
|
||||||
@ -384,7 +388,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
|
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
// __ ___
|
// __ ___
|
||||||
@ -399,7 +404,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
TrStap = Grid::trace (U[mu]*Stap);
|
TrStap = Grid::trace (U[mu]*Stap);
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
// --
|
// --
|
||||||
@ -423,7 +429,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
SumTrStap += TrStap;
|
SumTrStap += TrStap;
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -441,11 +448,13 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
TrStap = Grid::trace (U[mu]*Stap);
|
TrStap = Grid::trace (U[mu]*Stap);
|
||||||
trect = Grid::sum(TrStap);
|
trect = Grid::sum(TrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
trect = Grid::sum(SumTrStap);
|
trect = Grid::sum(SumTrStap);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
|
||||||
}
|
}
|
||||||
Peter.Stop();
|
Peter.Stop();
|
||||||
Azusa.Start();
|
Azusa.Start();
|
||||||
@ -489,7 +498,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
// __ __
|
// __ __
|
||||||
// |__ |
|
// |__ |
|
||||||
@ -501,7 +511,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
// __
|
// __
|
||||||
// |__ __ |
|
// |__ __ |
|
||||||
@ -513,7 +524,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
// __
|
// __
|
||||||
@ -526,7 +538,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
|
|
||||||
// 1(mu) x 2 (nu) ** this was ok
|
// 1(mu) x 2 (nu) ** this was ok
|
||||||
@ -542,7 +555,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
|
|
||||||
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
rrect = real(crect);
|
||||||
|
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
// 1(mu) x 2 (nu) ** this was ok
|
// 1(mu) x 2 (nu) ** this was ok
|
||||||
//
|
//
|
||||||
@ -570,8 +584,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
|||||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||||
trect = Grid::sum(RectPlaq_d);
|
trect = Grid::sum(RectPlaq_d);
|
||||||
crect = Grid::TensorRemove(trect);
|
crect = Grid::TensorRemove(trect);
|
||||||
|
rrect = real(crect);
|
||||||
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
|
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
|
||||||
|
|
||||||
}
|
}
|
||||||
Azusa.Stop();
|
Azusa.Stop();
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user