From 1c5b7a6be5d2d65d268cb0af048b963d078109ec Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Thu, 3 Nov 2016 16:26:56 +0000 Subject: [PATCH] Staggered phases first cut, c1, c2, u0 --- benchmarks/Benchmark_staggered.cc | 133 ++++++++++++++++++ lib/qcd/action/Actions.h | 3 + lib/qcd/action/fermion/FermionOperatorImpl.h | 48 ++++++- .../fermion/ImprovedStaggeredFermion.cc | 51 ++++++- .../action/fermion/ImprovedStaggeredFermion.h | 4 + 5 files changed, 227 insertions(+), 12 deletions(-) create mode 100644 benchmarks/Benchmark_staggered.cc diff --git a/benchmarks/Benchmark_staggered.cc b/benchmarks/Benchmark_staggered.cc new file mode 100644 index 00000000..c647cda0 --- /dev/null +++ b/benchmarks/Benchmark_staggered.cc @@ -0,0 +1,133 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./benchmarks/Benchmark_wilson.cc + + Copyright (C) 2015 + +Author: Peter Boyle +Author: paboyle + + 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 + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + + int threads = GridThread::GetThreads(); + std::cout< seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedRandomDevice(); + + typedef typename ImprovedStaggeredFermionR::FermionField FermionField; + typename ImprovedStaggeredFermionR::ImplParams params; + + FermionField src (&Grid); random(pRNG,src); + FermionField result(&Grid); result=zero; + FermionField ref(&Grid); ref=zero; + FermionField tmp(&Grid); tmp=zero; + FermionField err(&Grid); tmp=zero; + LatticeGaugeField Umu(&Grid); random(pRNG,Umu); + std::vector U(4,&Grid); + + double volume=1; + for(int mu=0;mu(Umu,U[nn],nn); + } +#endif + + for(int mu=0;mu(Umu,mu); + } + ref = zero; + /* + { // Naive wilson implementation + ref = zero; + for(int mu=0;mu GparityMobiusFermionR; typedef MobiusFermion GparityMobiusFermionF; typedef MobiusFermion GparityMobiusFermionD; +typedef ImprovedStaggeredFermion ImprovedStaggeredFermionR; +typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; +typedef ImprovedStaggeredFermion ImprovedStaggeredFermionD; }} /////////////////////////////////////////////////////////////////////////////// diff --git a/lib/qcd/action/fermion/FermionOperatorImpl.h b/lib/qcd/action/fermion/FermionOperatorImpl.h index 0b162f77..98d2f859 100644 --- a/lib/qcd/action/fermion/FermionOperatorImpl.h +++ b/lib/qcd/action/fermion/FermionOperatorImpl.h @@ -520,14 +520,17 @@ PARALLEL_FOR_LOOP INHERIT_GIMPL_TYPES(Gimpl); + template using iImplScalar = iScalar > >; template using iImplSpinor = iScalar > >; template using iImplHalfSpinor = iVector >, Ngp>; template using iImplDoubledGaugeField = iVector >, Nds>; + typedef iImplScalar SiteComplex; typedef iImplSpinor SiteSpinor; typedef iImplHalfSpinor SiteHalfSpinor; typedef iImplDoubledGaugeField SiteDoubledGaugeField; + typedef Lattice ComplexField; typedef Lattice FermionField; typedef Lattice DoubledGaugeField; @@ -564,15 +567,46 @@ PARALLEL_FOR_LOOP conformable(Uds._grid, GaugeGrid); conformable(Umu._grid, GaugeGrid); GaugeLinkField U(GaugeGrid); + GaugeLinkField UU(GaugeGrid); + GaugeLinkField UUU(GaugeGrid); + GaugeLinkField Udag(GaugeGrid); + GaugeLinkField UUUdag(GaugeGrid); for (int mu = 0; mu < Nd; mu++) { - U = PeekIndex(Umu, mu); + + // Staggered Phase. + ComplexField coor(GaugeGrid); + ComplexField phases(GaugeGrid); + ComplexField x(GaugeGrid); LatticeCoordinate(x,0); + ComplexField y(GaugeGrid); LatticeCoordinate(y,1); + ComplexField z(GaugeGrid); LatticeCoordinate(z,2); + ComplexField t(GaugeGrid); LatticeCoordinate(t,3); + + SiteComplex zz(0.0,0.0); + SiteComplex one(1.0,0.0); + + phases = one; + if ( mu == 1 ) phases = where( mod(x ,2)== zz, phases,-phases); + if ( mu == 2 ) phases = where( mod(x+y ,2)== zz, phases,-phases); + if ( mu == 3 ) phases = where( mod(x+y+z,2)== zz, phases,-phases); + + U = PeekIndex(Umu, mu); + UU = Gimpl::CovShiftForward(U,mu,U); + UUU= Gimpl::CovShiftForward(U,mu,UU); + + U = U *phases; + UUU = UUU *phases; + PokeIndex(Uds, U, mu); - PokeIndex(UUUds, U, mu); - std::cout << GridLogMessage << " NOT created the treble links for staggered yet" <(Uds, U, mu + 4); - PokeIndex(UUUds, U, mu+4); + PokeIndex(UUUds, UUU, mu); + + std::cout << GridLogMessage << " Created the treble links for staggered Naik term" <(Uds, Udag, mu + 4); + PokeIndex(UUUds, UUUdag, mu+4); + } } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 73cc272a..d0927222 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -43,6 +43,7 @@ ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, template ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, + RealD _c1, RealD _c2,RealD _u0, const ImplParams &p) : Kernels(p), _grid(&Fgrid), @@ -51,6 +52,9 @@ ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Umu, GridC StencilEven(&Hgrid, npoint, Even, directions, displacements), // source is Even StencilOdd(&Hgrid, npoint, Odd, directions, displacements), // source is Odd mass(_mass), + c1(_c1), + c2(_c2), + u0(_u0), Lebesgue(_grid), LebesgueEvenOdd(_cbgrid), Umu(&Fgrid), @@ -63,15 +67,52 @@ ImprovedStaggeredFermion::ImprovedStaggeredFermion(GaugeField &_Umu, GridC ImportGauge(_Umu); } + //////////////////////////////////////////////////////////// + // Momentum space propagator should be + // https://arxiv.org/pdf/hep-lat/9712010.pdf + // + // mom space action. + // gamma_mu i ( c1 sin pmu + c2 sin 3 pmu ) + m + // + // must track through staggered flavour/spin reduction in literature to + // turn to free propagator for the one component chi field, a la page 4/5 + // of above link to implmement fourier based solver. + //////////////////////////////////////////////////////////// template void ImprovedStaggeredFermion::ImportGauge(const GaugeField &_Umu) { - GaugeField HUmu(_Umu._grid); - HUmu = _Umu * (-0.5); - Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, HUmu); + + GaugeLinkField U(GaugeGrid); + + //////////////////////////////////////////////////////// + // Double Store should take two fields for Naik and one hop separately. + //////////////////////////////////////////////////////// + Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, _Umu); + + + //////////////////////////////////////////////////////// + // Apply scale factors to get the right fermion Kinetic term + // + // 0.5 ( U p(x+mu) - Udag(x-mu) p(x-mu) ) + //////////////////////////////////////////////////////// + for (int mu = 0; mu < Nd; mu++) { + + U = PeekIndex(Umu, mu); + PokeIndex(Umu, U*( 0.5*c1/u0), mu ); + + U = PeekIndex(Umu, mu+4); + PokeIndex(Umu, U*(-0.5*c1/u0), mu+4); + + U = PeekIndex(UUUmu, mu); + PokeIndex(UUUmu, U*( 0.5*c2/u0/u0/u0), mu ); + + U = PeekIndex(UUUmu, mu+4); + PokeIndex(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4); + } + pickCheckerboard(Even, UmuEven, Umu); - pickCheckerboard(Odd, UmuOdd, Umu); + pickCheckerboard(Odd, UmuOdd , Umu); pickCheckerboard(Even, UUUmuEven, UUUmu); - pickCheckerboard(Odd, UUUmuOdd, UUUmu); + pickCheckerboard(Odd, UUUmuOdd, UUUmu); } ///////////////////////////// diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index bf05240e..72f01bf3 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -108,6 +108,7 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS // Constructor ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass, + RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0, const ImplParams &p = ImplParams()); // DoubleStore impl dependent @@ -122,6 +123,9 @@ class ImprovedStaggeredFermion : public StaggeredKernels, public ImprovedS // any other parameters of action ??? RealD mass; + RealD u0; + RealD c1; + RealD c2; GridBase *_grid; GridBase *_cbgrid;