diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index f1b8820e..ee97e96d 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -183,11 +183,13 @@ namespace Grid { virtual RealD Mpc (const Field &in, Field &out) =0; virtual RealD MpcDag (const Field &in, Field &out) =0; virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { - Field tmp(in._grid); + Field tmp(in._grid); + tmp.checkerboard = in.checkerboard; ni=Mpc(in,tmp); no=MpcDag(tmp,out); } virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + out.checkerboard = in.checkerboard; MpcDagMpc(in,out,n1,n2); } virtual void HermOp(const Field &in, Field &out){ @@ -215,13 +217,15 @@ namespace Grid { public: SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){}; virtual RealD Mpc (const Field &in, Field &out) { - Field tmp(in._grid); -// std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; + Field tmp(in._grid); + tmp.checkerboard = !in.checkerboard; + //std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl; _Mat.Meooe(in,tmp); _Mat.MooeeInv(tmp,out); _Mat.Meooe(out,tmp); + //std::cout << "cb in " << in.checkerboard << " cb out " << out.checkerboard << std::endl; _Mat.Mooee(in,out); return axpy_norm(out,-1.0,tmp,out); } diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index fff970a2..3ec90e06 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -84,14 +84,14 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); // Compute the Clover Operator acting on Colour and Spin - CloverTerm = fillCloverYZ(Bx); - CloverTerm += fillCloverXZ(By); - CloverTerm += fillCloverXY(Bz); - CloverTerm += fillCloverXT(Ex); - CloverTerm += fillCloverYT(Ey); - CloverTerm += fillCloverZT(Ez); - CloverTerm *= (0.5) * csw; - CloverTerm += (4.0 + this->mass); + // multiply here by the clover coefficients for the anisotropy + CloverTerm = fillCloverYZ(Bx) * csw_r; + CloverTerm += fillCloverXZ(By) * csw_r; + CloverTerm += fillCloverXY(Bz) * csw_r; + CloverTerm += fillCloverXT(Ex) * csw_t; + CloverTerm += fillCloverYT(Ey) * csw_t; + CloverTerm += fillCloverZT(Ez) * csw_t; + CloverTerm += diag_mass; int lvol = _Umu._grid->lSites(); int DimRep = Impl::Dimension; @@ -145,7 +145,6 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) template void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) { - conformable(in, out); this->MooeeInternal(in, out, DaggerNo, InverseNo); } @@ -158,14 +157,12 @@ void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &o template void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) { - conformable(in,out); this->MooeeInternal(in, out, DaggerNo, InverseYes); } template void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) { - conformable(in,out); this->MooeeInternal(in, out, DaggerYes, InverseYes); } @@ -228,88 +225,7 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie template void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) { - - GridBase *grid = mat._grid; - - //GaugeLinkField Lambdaodd(grid), Lambdaeven(grid), tmp(grid); - //Lambdaodd = zero; //Yodd*dag(Xodd)+Xodd*dag(Yodd); // I have to peek spin and decide the color structure - //Lambdaeven = zero; //Teven*dag(Xeven)+Xeven*dag(Yeven) + 2*(Dee^-1) - - GaugeLinkField Lambda(grid), tmp(grid); - Lambda = zero; - - conformable(mat._grid, X._grid); - conformable(Y._grid, X._grid); - - std::vector C1p(Nd, grid), C2p(Nd, grid), C3p(Nd, grid), C4p(Nd, grid); - std::vector C1m(Nd, grid), C2m(Nd, grid), C3m(Nd, grid), C4m(Nd, grid); - std::vector U(Nd, mat._grid); - - for (int mu = 0; mu < Nd; mu++) - { - U[mu] = PeekIndex(mat, mu); - C1p[mu] = zero; - C2p[mu] = zero; - C3p[mu] = zero; - C4p[mu] = zero; - C1m[mu] = zero; - C2m[mu] = zero; - C3m[mu] = zero; - C4m[mu] = zero; - } - - /* - PARALLEL_FOR_LOOP - for (int i = 0; i < CloverTerm._grid->oSites(); i++) - { - T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()()); - T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()()); - T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); - T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); - } -*/ - - for (int i = 0; i < 4; i++) - { //spin - for (int j = 0; j < 4; j++) - { //spin - - for (int mu = 0; mu < 4; mu++) - { //color - for (int nu = 0; nu < 4; nu++) - { //color - - // insertion in upper staple - tmp = Lambda * U[nu]; - C1p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - - tmp = Lambda * U[mu]; - C2p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - - tmp = Impl::CovShiftIdentityForward(Lambda, nu) * U[nu]; - C3p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); - - tmp = Lambda; - C4p[mu] += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * tmp; - - // insertion in lower staple - tmp = Lambda * U[nu]; - C1m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); - - tmp = Lambda * U[mu]; - C2m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, U[nu])), mu); - - tmp = Lambda * U[nu]; - C3m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); - - tmp = Lambda; - C4m[mu] += Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * tmp; - } - } - } - } - - //Still implementing. Have to be tested, and understood how to project EO + assert(0); } // Derivative parts diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index cd13b225..268564c0 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -6,8 +6,8 @@ Copyright (C) 2017 - Author: paboyle Author: Guido Cossu + Author: David Preti <> 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 @@ -37,6 +37,22 @@ namespace Grid namespace QCD { +/////////////////////////////////////////////////////////////////// +// Wilson Clover +// +// Operator ( with anisotropy coefficients): +// +// Q = 1 + (Nd-1)/xi_0 + m +// + W_t + (nu/xi_0) * W_s +// - 1/2*[ csw_t * sum_s (sigma_ts F_ts) + (csw_s/xi_0) * sum_ss (sigma_ss F_ss) ] +// +// s spatial, t temporal directions. +// where W_t and W_s are the temporal and spatial components of the +// Wilson Dirac operator +// +// csw_r = csw_t to recover the isotropic version +////////////////////////////////////////////////////////////////// + template class WilsonCloverFermion : public WilsonFermion { @@ -55,28 +71,43 @@ public: // Constructors WilsonCloverFermion(GaugeField &_Umu, GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, - RealD _mass, - RealD _csw, - const ImplParams &p = ImplParams()) : WilsonFermion(_Umu, - Fgrid, - Hgrid, - _mass, p), - CloverTerm(&Fgrid), - CloverTermInv(&Fgrid), - CloverTermEven(&Hgrid), - CloverTermOdd(&Hgrid), - CloverTermInvEven(&Hgrid), - CloverTermInvOdd(&Hgrid), - CloverTermDagEven(&Hgrid), - CloverTermDagOdd(&Hgrid), - CloverTermInvDagEven(&Hgrid), - CloverTermInvDagOdd(&Hgrid) + const RealD _mass, + const RealD _csw_r = 0.0, + const RealD _csw_t = 0.0, + const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(), + const ImplParams &impl_p = ImplParams()) : WilsonFermion(_Umu, + Fgrid, + Hgrid, + _mass, impl_p, clover_anisotropy), + CloverTerm(&Fgrid), + CloverTermInv(&Fgrid), + CloverTermEven(&Hgrid), + CloverTermOdd(&Hgrid), + CloverTermInvEven(&Hgrid), + CloverTermInvOdd(&Hgrid), + CloverTermDagEven(&Hgrid), + CloverTermDagOdd(&Hgrid), + CloverTermInvDagEven(&Hgrid), + CloverTermInvDagOdd(&Hgrid) { - csw = _csw; assert(Nd == 4); // require 4 dimensions - if (csw == 0) - std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl; + if (clover_anisotropy.isAnisotropic) + { + csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0; + diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); + } + else + { + csw_r = _csw_r * 0.5; + diag_mass = 4.0 + _mass; + } + csw_t = _csw_t * 0.5; + + if (csw_r == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl; + if (csw_t == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl; ImportGauge(_Umu); } @@ -105,15 +136,15 @@ public: GaugeField clover_force(force._grid); PropagatorField Lambda(force._grid); - // Here we are hitting some performance issues: + // Guido: Here we are hitting some performance issues: // need to extract the components of the DoubledGaugeField // for each call // Possible solution // Create a vector object to store them? (cons: wasting space) std::vector U(Nd, this->Umu._grid); - + Impl::extractLinkField(U, this->Umu); - + force = zero; // Derivative of the Wilson hopping term this->DhopDeriv(force, X, Y, dag); @@ -121,10 +152,9 @@ public: /////////////////////////////////////////////////////////// // Clover term derivative /////////////////////////////////////////////////////////// - Impl::outerProductImpl(Lambda, X, Y); + Impl::outerProductImpl(Lambda, X, Y); //std::cout << "Lambda:" << Lambda << std::endl; - Gamma::Algebra sigma[] = { Gamma::Algebra::SigmaXY, Gamma::Algebra::SigmaXZ, @@ -148,25 +178,34 @@ public: */ int count = 0; - clover_force = zero; + clover_force = zero; for (int mu = 0; mu < 4; mu++) { force_mu = zero; for (int nu = 0; nu < 4; nu++) { - if (mu == nu) continue; + if (mu == nu) + continue; + + RealD factor; + if (nu == 4 || mu == 4) + { + factor = 2.0 * csw_t; + } + else + { + factor = 2.0 * csw_r; + } PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked - Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok - force_mu -= Cmunu(U, lambda, mu, nu); // checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= factor*Cmunu(U, lambda, mu, nu); // checked count++; } pokeLorentz(clover_force, U[mu] * force_mu, mu); } - clover_force *= csw; + //clover_force *= csw; force += clover_force; - - } // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis @@ -176,15 +215,15 @@ public: GaugeLinkField out(lambda._grid), tmp(lambda._grid); // insertion in upper staple // please check redundancy of shift operations - + // C1+ tmp = lambda * U[nu]; out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - + // C2+ tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu); out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - + // C3+ tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu); out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); @@ -213,16 +252,17 @@ public: private: // here fixing the 4 dimensions, make it more general? - RealD csw; // Clover coefficient + RealD csw_r; // Clover coefficient - spatial + RealD csw_t; // Clover coefficient - temporal + RealD diag_mass; // Mass term CloverFieldType CloverTerm, CloverTermInv; // Clover term CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO - // eventually these two can be compressed into 6x6 blocks instead of the 12x12 + // eventually these can be compressed into 6x6 blocks instead of the 12x12 // using the DeGrand-Rossi basis for the gamma matrices - CloverFieldType fillCloverYZ(const GaugeLinkField &F) { CloverFieldType T(F._grid); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 7b2b5206..ca5eba8b 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -45,10 +45,11 @@ class WilsonFermionStatic { }; struct WilsonAnisotropyCoefficients{ - bool isAnisotropic; - int t_direction; - double xi_0; - double nu; + GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients, + bool, isAnisotropic, + int, t_direction, + double, xi_0, + double, nu); WilsonAnisotropyCoefficients(): isAnisotropic(false), diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc index 9a55f6b2..9281e298 100644 --- a/tests/core/Test_wilson_clover.cc +++ b/tests/core/Test_wilson_clover.cc @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -24,302 +24,334 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ - /* END LEGAL */ +/* END LEGAL */ #include using namespace std; using namespace Grid; using namespace Grid::QCD; -int main (int argc, char ** argv) +int main(int argc, char **argv) { - Grid_init(&argc,&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); + 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(&Grid); int threads = GridThread::GetThreads(); - std::cout< seeds({1,2,3,4}); - GridParallelRNG pRNG(&Grid); + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); - typedef typename WilsonCloverFermionR::FermionField FermionField; - typename WilsonCloverFermionR::ImplParams params; + typedef typename WilsonCloverFermionR::FermionField FermionField; + typename WilsonCloverFermionR::ImplParams params; + WilsonAnisotropyCoefficients anis; - FermionField src (&Grid); random(pRNG,src); - FermionField result(&Grid); result=zero; - FermionField result2(&Grid); result2=zero; - FermionField ref(&Grid); ref=zero; - FermionField tmp(&Grid); tmp=zero; - FermionField err(&Grid); err=zero; - FermionField err2(&Grid); err2=zero; - FermionField phi (&Grid); random(pRNG,phi); - FermionField chi (&Grid); random(pRNG,chi); - LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); - std::vector U(4,&Grid); + FermionField src(&Grid); + random(pRNG, src); + FermionField result(&Grid); + result = zero; + FermionField result2(&Grid); + result2 = zero; + FermionField ref(&Grid); + ref = zero; + FermionField tmp(&Grid); + tmp = zero; + FermionField err(&Grid); + err = zero; + FermionField err2(&Grid); + err2 = zero; + FermionField phi(&Grid); + random(pRNG, phi); + FermionField chi(&Grid); + random(pRNG, chi); + LatticeGaugeField Umu(&Grid); + SU3::HotConfiguration(pRNG, Umu); + std::vector U(4, &Grid); - - double volume=1; - for(int mu=0;mu * = < chi | Deo^dag| phi> "< * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; - FermionField dchi_e (&RBGrid); - FermionField dchi_o (&RBGrid); + FermionField chi_e(&RBGrid); + FermionField chi_o(&RBGrid); - FermionField phi_e (&RBGrid); - FermionField phi_o (&RBGrid); + FermionField dchi_e(&RBGrid); + FermionField dchi_o(&RBGrid); - FermionField dphi_e (&RBGrid); - FermionField dphi_o (&RBGrid); + FermionField phi_e(&RBGrid); + FermionField phi_o(&RBGrid); - pickCheckerboard(Even,chi_e,chi); - pickCheckerboard(Odd ,chi_o,chi); - pickCheckerboard(Even,phi_e,phi); - pickCheckerboard(Odd ,phi_o,phi); + FermionField dphi_e(&RBGrid); + FermionField dphi_o(&RBGrid); - Dwc.Meooe(chi_e,dchi_o); - Dwc.Meooe(chi_o,dchi_e); - Dwc.MeooeDag(phi_e,dphi_o); - Dwc.MeooeDag(phi_o,dphi_e); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); - 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); + Dwc.Meooe(chi_e, dchi_o); + Dwc.Meooe(chi_o, dchi_e); + Dwc.MeooeDag(phi_e, dphi_o); + Dwc.MeooeDag(phi_o, dphi_e); - std::cout< seeds2({5,6,7,8}); - GridParallelRNG pRNG2(&Grid); pRNG2.SeedFixedIntegers(seeds2); + std::vector seeds2({5, 6, 7, 8}); + GridParallelRNG pRNG2(&Grid); + pRNG2.SeedFixedIntegers(seeds2); LatticeColourMatrix Omega(&Grid); LatticeColourMatrix ShiftedOmega(&Grid); - LatticeGaugeField U_prime(&Grid); U_prime=zero; - LatticeColourMatrix U_prime_mu(&Grid); U_prime_mu=zero; + LatticeGaugeField U_prime(&Grid); + U_prime = zero; + LatticeColourMatrix U_prime_mu(&Grid); + U_prime_mu = zero; SU::LieRandomize(pRNG2, Omega, 1.0); - for (int mu=0;mu + +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 + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + // Possibile to create the module by hand + // hardcoding parameters or using a Reader + + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6 ; + WilsonGaugeActionR Waction(beta); + + // temporarily need a gauge field + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + + ConjugateGradient CG(1.0e-8, 2000); + + TwoFlavourEvenOddPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main + + + + + + + diff --git a/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc new file mode 100644 index 00000000..322bb304 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc @@ -0,0 +1,126 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2017 + +Author: Guido Cossu + +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 + +int main(int argc, char **argv) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + typedef WilsonImplR FermionImplPolicy; + typedef WilsonCloverFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_lat"; + CPparams.rng_prefix = "ckpoint_rng"; + CPparams.saveInterval = 5; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + // need wrappers of the fermionic classes + // that have a complex construction + // standard + RealD beta = 5.6; + WilsonGaugeActionR Waction(beta); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + Real mass = 0.01; + Real csw = 1.0; + + FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw); + ConjugateGradient CG(1.0e-8, 5000); + + TwoFlavourPseudoFermionAction Nf2(FermOp, CG, CG); + + // Set smearing (true/false), default: false + Nf2.is_smeared = false; + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Nf2); + + ActionLevel Level2(4); + Level2.push_back(&Waction); + + TheHMC.TheAction.push_back(Level1); + TheHMC.TheAction.push_back(Level2); + ///////////////////////////////////////////////////////////// + + /* + double rho = 0.1; // smearing parameter + int Nsmear = 2; // number of smearing levels + Smear_Stout Stout(rho); + SmearedConfiguration SmearingPolicy( + UGrid, Nsmear, Stout); + */ + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + TheHMC.Run(); // no smearing + // TheHMC.Run(SmearingPolicy); // for smearing + + Grid_finalize(); + +} // main diff --git a/tests/qdpxx/Test_qdpxx_wilson.cc b/tests/qdpxx/Test_qdpxx_wilson.cc index a084ebad..29e9c9ce 100644 --- a/tests/qdpxx/Test_qdpxx_wilson.cc +++ b/tests/qdpxx/Test_qdpxx_wilson.cc @@ -274,10 +274,10 @@ public: Chroma::CloverFermActParams p; p.Mass = _mq; p.clovCoeffR = QDP::Real(1.0); - p.clovCoeffT = QDP::Real(1.0); + p.clovCoeffT = QDP::Real(2.0); p.u0 = QDP::Real(1.0); AnisoParam_t _apar; - _apar.anisoP = false; + _apar.anisoP = true; _apar.t_dir = 3; // in 4d _apar.xi_0 = 2.0; _apar.nu = 1.0; @@ -422,10 +422,15 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD if (action == WilsonClover) { - Grid::RealD _csw = 1.0; + Grid::RealD _csw_r = 1.0; + Grid::RealD _csw_t = 2.0; WilsonAnisotropyCoefficients anis; - WilsonImplParams implParam; - Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw); + anis.isAnisotropic = true; + anis.t_direction = 3; + anis.xi_0 = 2.0; + anis.nu = 1.0; + WilsonImplParams CloverImplParam; + Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam); Wf.ImportGauge(Umu); std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;