From 82aecbf4cfdb76a44c90855292e5c08628e48b62 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Tue, 5 Apr 2022 15:26:39 +0100 Subject: [PATCH] Test_wilson_conserved_current added --- tests/core/Test_wilson_conserved_current.cc | 290 ++++++++++++++++++++ 1 file changed, 290 insertions(+) create mode 100644 tests/core/Test_wilson_conserved_current.cc diff --git a/tests/core/Test_wilson_conserved_current.cc b/tests/core/Test_wilson_conserved_current.cc new file mode 100644 index 00000000..3eba4fdc --- /dev/null +++ b/tests/core/Test_wilson_conserved_current.cc @@ -0,0 +1,290 @@ +/************************************************************************************* +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_cayley_cg.cc + +Copyright (C) 2022 + +Author: Peter Boyle +Author: Fabian Joswig + +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 +#include + +using namespace std; +using namespace Grid; + + +template +void TestConserved(What & Ddwf, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5, + What *Ddwfrev=nullptr); + + 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< omegas; + std::vector < ComplexD > omegasrev(Ls); + + 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); + + + GridCartesian * UGridF = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd,vComplexF::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian * UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF); + GridCartesian * FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF); + GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF); + + + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + std::vector seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + if( argc > 1 && argv[1][0] != '-' ) + { + std::cout<::HotConfiguration(RNG4,Umu); + } + + RealD mass=0.3; + RealD M5 =1.0; + std::cout<(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + std::cout<(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} + + + +template +void TestConserved(Action & Ddwf, + LatticeGaugeField &Umu, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5, + Action * Ddwfrev) +{ + LatticePropagator phys_src(UGrid); + LatticePropagator seqsrc(FGrid); + LatticePropagator prop5(FGrid); + LatticePropagator prop5rev(FGrid); + LatticePropagator prop4(UGrid); + LatticePropagator Axial_mu(UGrid); + LatticePropagator Vector_mu(UGrid); + LatticeComplex PA (UGrid); + LatticeComplex SV (UGrid); + LatticeComplex VV (UGrid); + LatticeComplex PJ5q(UGrid); + LatticeComplex PP (UGrid); + LatticePropagator seqprop(UGrid); + + SpinColourMatrix kronecker; kronecker=1.0; + Coordinate coor({0,0,0,0}); + phys_src=Zero(); + pokeSite(kronecker,phys_src,coor); + + ConjugateGradient CG(1.0e-16,100000); + SchurRedBlackDiagTwoSolve schur(CG); + ZeroGuesser zpg; + for(int s=0;s(src4,phys_src,s,c); + + LatticeFermion src5 (FGrid); + Ddwf.ImportPhysicalFermionSource(src4,src5); + + LatticeFermion result5(FGrid); result5=Zero(); + schur(Ddwf,src5,result5,zpg); + std::cout<(prop5,result5,s,c); + + LatticeFermion result4(UGrid); + Ddwf.ExportPhysicalFermionSolution(result5,result4); + FermToProp(prop4,result4,s,c); + + if( Ddwfrev ) { + Ddwfrev->ImportPhysicalFermionSource(src4,src5); + result5 = Zero(); + schur(*Ddwfrev,src5,result5,zpg); + } + FermToProp(prop5rev,result5,s,c); + } + } + + auto curr = Current::Vector; + const int mu_J=0; + const int t_J=0; + + LatticeComplex ph (UGrid); ph=1.0; + + Ddwf.SeqConservedCurrent(prop5, + seqsrc, + phys_src, + curr, + mu_J, + t_J, + t_J,// whole lattice + ph); + + for(int s=0;s(src5,seqsrc,s,c); + + LatticeFermion result5(FGrid); result5=Zero(); + schur(Ddwf,src5,result5,zpg); + + LatticeFermion result4(UGrid); + Ddwf.ExportPhysicalFermionSolution(result5,result4); + FermToProp(seqprop,result4,s,c); + } + } + + Gamma g5(Gamma::Algebra::Gamma5); + Gamma gT(Gamma::Algebra::GammaT); + + std::vector sumPA; + std::vector sumSV; + std::vector sumVV; + std::vector sumPP; + std::vector sumPJ5q; + + Ddwf.ContractConservedCurrent(prop5rev,prop5,Axial_mu,phys_src,Current::Axial,Tdir); + Ddwf.ContractConservedCurrent(prop5rev,prop5,Vector_mu,phys_src,Current::Vector,Tdir); + Ddwf.ContractJ5q(prop5,PJ5q); + + PA = trace(g5*Axial_mu); // Pseudoscalar-Axial conserved current + SV = trace(Vector_mu); // Scalar-Vector conserved current + VV = trace(gT*Vector_mu); // (local) Vector-Vector conserved current + PP = trace(adj(prop4)*prop4); // Pseudoscalar density + + // Spatial sum + sliceSum(PA,sumPA,Tdir); + sliceSum(SV,sumSV,Tdir); + sliceSum(VV,sumVV,Tdir); + sliceSum(PP,sumPP,Tdir); + sliceSum(PJ5q,sumPJ5q,Tdir); + + const int Nt{static_cast(sumPA.size())}; + std::cout< check_buf; + + test_S = trace(qSite*g); + test_V = trace(qSite*g*Gamma::gmu[mu_J]); + + Ddwf.ContractConservedCurrent(prop5rev,prop5,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<