diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index c27254aa..8f80f08f 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,25 +1,26 @@ -#include -#include -#include -#include -#include -#include +#include +#include #include +#include +#include #include +#include #include #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include #include #include #include +#include +#include +#include +#include +#include +#include #include -#include -#include -#include -#include -#include diff --git a/extras/Hadrons/Modules/MSource/Laplacian.hpp b/extras/Hadrons/Modules/MSource/Laplacian.hpp new file mode 100644 index 00000000..d5f866dd --- /dev/null +++ b/extras/Hadrons/Modules/MSource/Laplacian.hpp @@ -0,0 +1,153 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/Laplacian.hpp + +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 */ + +#ifndef Hadrons_MSource_Laplacian_hpp_ +#define Hadrons_MSource_Laplacian_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Laplacian smearing source + ----------------------------- + + * options: + - source: name of source object to be smeared (string) + - N: number of steps (integer) + - alpha: smearing parameter (real) + + */ + +/****************************************************************************** + * Z2 stochastic source * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class LaplacianPar : Serializable +{ + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianPar, + std::string, source, + std::string, gauge, + unsigned int, N, + double, alpha); +}; + +template +class TLaplacian : public Module +{ + public: + FERM_TYPE_ALIASES(FImpl, ); + + public: + // constructor + TLaplacian(const std::string name); + // destructor + virtual ~TLaplacian(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Lap, TLaplacian, MSource); + +/****************************************************************************** + * TLaplacian template implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TLaplacian::TLaplacian(const std::string name) + : Module(name) +{ +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TLaplacian::getInput(void) +{ + std::vector in = {par().source, par().gauge}; + + return in; +} + +template +std::vector TLaplacian::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TLaplacian::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TLaplacian::execute(void) +{ + + FermionField source(env().getGrid()), tmp(env().getGrid()); + PropagatorField &SmrSrc = *env().template createLattice(getName()); + PropagatorField &fullSrc = *env().template getObject(par().source); + auto &U = *env().template getObject(par().gauge); + Laplacian LaplaceOperator(env().getGrid()); + LaplaceOperator.ImportGauge(U); + double prefactor = par().alpha / (double)(par().N); + + for (unsigned int s = 0; s < Ns; ++s) + { + for (unsigned int c = 0; c < Nc; ++c) + { + PropToFerm(source, fullSrc, s, c); + for (int smr = 0; smr < par().N; ++smr) + { + LaplaceOperator.M(source, tmp); + source += prefactor * tmp; + } + FermToProp(SmrSrc, source, s, c); + } + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSource_Z2_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 669b08ba..b3614499 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,38 +1,39 @@ modules_cc =\ - Modules/MContraction/WeakHamiltonianEye.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \ - Modules/MGauge/Load.cc \ + Modules/MContraction/WeakHamiltonianEye.cc \ + Modules/MScalar/FreeProp.cc \ + Modules/MScalar/ChargedProp.cc \ + Modules/MGauge/Unit.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/StochEm.cc \ - Modules/MGauge/Unit.cc \ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/FreeProp.cc + Modules/MGauge/Load.cc modules_hpp =\ - Modules/MAction/DWF.hpp \ - Modules/MAction/Wilson.hpp \ - Modules/MContraction/Baryon.hpp \ - Modules/MContraction/DiscLoop.hpp \ - Modules/MContraction/Gamma3pt.hpp \ - Modules/MContraction/Meson.hpp \ + Modules/MLoop/NoiseLoop.hpp \ + Modules/MFermion/GaugeProp.hpp \ Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MContraction/Baryon.hpp \ Modules/MContraction/WeakHamiltonianNonEye.hpp \ Modules/MContraction/WeakNeutral4ptDisc.hpp \ - Modules/MFermion/GaugeProp.hpp \ - Modules/MGauge/Load.hpp \ - Modules/MGauge/Random.hpp \ - Modules/MGauge/StochEm.hpp \ - Modules/MGauge/Unit.hpp \ - Modules/MLoop/NoiseLoop.hpp \ + Modules/MContraction/Gamma3pt.hpp \ + Modules/MSource/Z2.hpp \ + Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Point.hpp \ + Modules/MSource/Wall.hpp \ + Modules/MSource/Laplacian.hpp \ + Modules/MSolver/RBPrecCG.hpp \ Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ Modules/MScalar/Scalar.hpp \ - Modules/MSink/Point.hpp \ - Modules/MSolver/RBPrecCG.hpp \ - Modules/MSource/Point.hpp \ - Modules/MSource/SeqGamma.hpp \ - Modules/MSource/Wall.hpp \ - Modules/MSource/Z2.hpp + Modules/MAction/DWF.hpp \ + Modules/MAction/Wilson.hpp \ + Modules/MGauge/StochEm.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MGauge/Load.hpp \ + Modules/MSink/Point.hpp diff --git a/lib/qcd/action/Action.h b/lib/qcd/action/Action.h index 7272c90d..1297939d 100644 --- a/lib/qcd/action/Action.h +++ b/lib/qcd/action/Action.h @@ -47,4 +47,8 @@ Author: paboyle //////////////////////////////////////// #include +//////////////////////////////////////////////////////////////////////// +// Laplacian on fermion fields +//////////////////////////////////////////////////////////////////////// +#include #endif diff --git a/lib/qcd/action/ActionCore.h b/lib/qcd/action/ActionCore.h index 7a5caf15..8c7e8bc7 100644 --- a/lib/qcd/action/ActionCore.h +++ b/lib/qcd/action/ActionCore.h @@ -53,7 +53,7 @@ directory // Utility functions //////////////////////////////////////////// #include -#include +#include diff --git a/lib/qcd/action/fermion/DomainWallEOFAFermion.cc b/lib/qcd/action/fermion/DomainWallEOFAFermion.cc index dd8a500d..5c90e6b5 100644 --- a/lib/qcd/action/fermion/DomainWallEOFAFermion.cc +++ b/lib/qcd/action/fermion/DomainWallEOFAFermion.cc @@ -60,11 +60,11 @@ namespace QCD { Approx::zolotarev_free(zdata); } - /*************************************************************** - /* Additional EOFA operators only called outside the inverter. - /* Since speed is not essential, simple axpby-style - /* implementations should be fine. - /***************************************************************/ + /* + Additional EOFA operators only called outside the inverter. + Since speed is not essential, simple axpby-style + implementations should be fine. + */ template void DomainWallEOFAFermion::Omega(const FermionField& psi, FermionField& Din, int sign, int dag) { @@ -115,9 +115,9 @@ namespace QCD { return(norm2(chi)); } - /******************************************************************** - /* Performance critical fermion operators called inside the inverter - /********************************************************************/ + + // Performance critical fermion operators called inside the inverter + template void DomainWallEOFAFermion::M5D(const FermionField& psi, FermionField& chi) diff --git a/lib/qcd/action/fermion/MobiusEOFAFermion.cc b/lib/qcd/action/fermion/MobiusEOFAFermion.cc index 085fa988..8c1fd7a8 100644 --- a/lib/qcd/action/fermion/MobiusEOFAFermion.cc +++ b/lib/qcd/action/fermion/MobiusEOFAFermion.cc @@ -77,11 +77,11 @@ namespace QCD { } } - /*************************************************************** - /* Additional EOFA operators only called outside the inverter. - /* Since speed is not essential, simple axpby-style - /* implementations should be fine. - /***************************************************************/ + /* + Additional EOFA operators only called outside the inverter. + Since speed is not essential, simple axpby-style + implementations should be fine. + */ template void MobiusEOFAFermion::Omega(const FermionField& psi, FermionField& Din, int sign, int dag) { @@ -193,9 +193,9 @@ namespace QCD { return(norm2(chi)); } - /******************************************************************** - /* Performance critical fermion operators called inside the inverter - /********************************************************************/ + + // Performance critical fermion operators called inside the inverter + template void MobiusEOFAFermion::M5D(const FermionField& psi, FermionField& chi) diff --git a/lib/qcd/utils/CovariantAdjointLaplacian.h b/lib/qcd/utils/CovariantAdjointLaplacian.h new file mode 100644 index 00000000..1d3a2053 --- /dev/null +++ b/lib/qcd/utils/CovariantAdjointLaplacian.h @@ -0,0 +1,209 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/scalar/CovariantAdjointLaplacian.h + +Copyright (C) 2016 + +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 */ + +#ifndef COVARIANT_ADJOINT_LAPLACIAN_H +#define COVARIANT_ADJOINT_LAPLACIAN_H + +namespace Grid +{ +namespace QCD +{ + +struct LaplacianParams : Serializable +{ + GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams, + RealD, lo, + RealD, hi, + int, MaxIter, + RealD, tolerance, + int, degree, + int, precision); + + // constructor + LaplacianParams(RealD lo = 0.0, + RealD hi = 1.0, + int maxit = 1000, + RealD tol = 1.0e-8, + int degree = 10, + int precision = 64) + : lo(lo), + hi(hi), + MaxIter(maxit), + tolerance(tol), + degree(degree), + precision(precision){}; +}; + +//////////////////////////////////////////////////////////// +// Laplacian operator L on adjoint fields +// +// phi: adjoint field +// L: D_mu^dag D_mu +// +// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag + +// U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu) +// -2phi(x)] +// +// Operator designed to be encapsulated by +// an HermitianLinearOperator<.. , ..> +//////////////////////////////////////////////////////////// + +template +class LaplacianAdjointField : public Metric +{ + OperatorFunction &Solver; + LaplacianParams param; + MultiShiftFunction PowerHalf; + MultiShiftFunction PowerInvHalf; + +public: + INHERIT_GIMPL_TYPES(Impl); + + LaplacianAdjointField(GridBase *grid, OperatorFunction &S, LaplacianParams &p, const RealD k = 1.0) + : U(Nd, grid), Solver(S), param(p), kappa(k) + { + AlgRemez remez(param.lo, param.hi, param.precision); + std::cout << GridLogMessage << "Generating degree " << param.degree << " for x^(1/2)" << std::endl; + remez.generateApprox(param.degree, 1, 2); + PowerHalf.Init(remez, param.tolerance, false); + PowerInvHalf.Init(remez, param.tolerance, true); + }; + + void Mdir(const GaugeField &, GaugeField &, int, int) { assert(0); } + void Mdiag(const GaugeField &, GaugeField &) { assert(0); } + + void ImportGauge(const GaugeField &_U) + { + for (int mu = 0; mu < Nd; mu++) + { + U[mu] = PeekIndex(_U, mu); + } + } + + void M(const GaugeField &in, GaugeField &out) + { + // in is an antihermitian matrix + // test + //GaugeField herm = in + adj(in); + //std::cout << "AHermiticity: " << norm2(herm) << std::endl; + + GaugeLinkField tmp(in._grid); + GaugeLinkField tmp2(in._grid); + GaugeLinkField sum(in._grid); + + for (int nu = 0; nu < Nd; nu++) + { + sum = zero; + GaugeLinkField in_nu = PeekIndex(in, nu); + GaugeLinkField out_nu(out._grid); + for (int mu = 0; mu < Nd; mu++) + { + tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]); + tmp2 = adj(U[mu]) * in_nu * U[mu]; + sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu; + } + out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum; + PokeIndex(out, out_nu, nu); + } + } + + void MDeriv(const GaugeField &in, GaugeField &der) + { + // in is anti-hermitian + RealD factor = -kappa / (double(4 * Nd)); + + for (int mu = 0; mu < Nd; mu++) + { + GaugeLinkField der_mu(der._grid); + der_mu = zero; + for (int nu = 0; nu < Nd; nu++) + { + GaugeLinkField in_nu = PeekIndex(in, nu); + der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu; + } + // the minus sign comes by using the in_nu instead of the + // adjoint in the last multiplication + PokeIndex(der, -2.0 * factor * der_mu, mu); + } + } + + // separating this temporarily + void MDeriv(const GaugeField &left, const GaugeField &right, + GaugeField &der) + { + // in is anti-hermitian + RealD factor = -kappa / (double(4 * Nd)); + + for (int mu = 0; mu < Nd; mu++) + { + GaugeLinkField der_mu(der._grid); + der_mu = zero; + for (int nu = 0; nu < Nd; nu++) + { + GaugeLinkField left_nu = PeekIndex(left, nu); + GaugeLinkField right_nu = PeekIndex(right, nu); + der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu; + der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu; + } + PokeIndex(der, -factor * der_mu, mu); + } + } + + void Minv(const GaugeField &in, GaugeField &inverted) + { + HermitianLinearOperator, GaugeField> HermOp(*this); + Solver(HermOp, in, inverted); + } + + void MSquareRoot(GaugeField &P) + { + GaugeField Gp(P._grid); + HermitianLinearOperator, GaugeField> HermOp(*this); + ConjugateGradientMultiShift msCG(param.MaxIter, PowerHalf); + msCG(HermOp, P, Gp); + P = Gp; + } + + void MInvSquareRoot(GaugeField &P) + { + GaugeField Gp(P._grid); + HermitianLinearOperator, GaugeField> HermOp(*this); + ConjugateGradientMultiShift msCG(param.MaxIter, PowerInvHalf); + msCG(HermOp, P, Gp); + P = Gp; + } + +private: + RealD kappa; + std::vector U; +}; + +} // QCD +} // Grid +#endif diff --git a/lib/qcd/utils/CovariantLaplacian.h b/lib/qcd/utils/CovariantLaplacian.h index 0c99b03e..188ab2d5 100644 --- a/lib/qcd/utils/CovariantLaplacian.h +++ b/lib/qcd/utils/CovariantLaplacian.h @@ -30,168 +30,57 @@ directory #ifndef COVARIANT_LAPLACIAN_H #define COVARIANT_LAPLACIAN_H -namespace Grid { -namespace QCD { - -struct LaplacianParams : Serializable { - GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams, - RealD, lo, - RealD, hi, - int, MaxIter, - RealD, tolerance, - int, degree, - int, precision); - - // constructor - LaplacianParams(RealD lo = 0.0, - RealD hi = 1.0, - int maxit = 1000, - RealD tol = 1.0e-8, - int degree = 10, - int precision = 64) - : lo(lo), - hi(hi), - MaxIter(maxit), - tolerance(tol), - degree(degree), - precision(precision){}; -}; - - +namespace Grid +{ +namespace QCD +{ //////////////////////////////////////////////////////////// -// Laplacian operator L on adjoint fields +// Laplacian operator L on fermion fields // -// phi: adjoint field -// L: D_mu^dag D_mu +// phi: fermion field // -// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag + -// U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu) -// -2phi(x)] +// L phi(x) = Sum_mu [ U_mu(x) phi(x+mu) + U_mu(x-mu) phi(x-mu) - 2phi(x)] // // Operator designed to be encapsulated by // an HermitianLinearOperator<.. , ..> //////////////////////////////////////////////////////////// +// has to inherit from a fermion implementation template -class LaplacianAdjointField: public Metric { - OperatorFunction &Solver; - LaplacianParams param; - MultiShiftFunction PowerHalf; - MultiShiftFunction PowerInvHalf; +class Laplacian +{ +public: + INHERIT_IMPL_TYPES(Impl); - public: - INHERIT_GIMPL_TYPES(Impl); + // add a bool to smear only in the spatial directions + Laplacian(GridBase *grid, bool spatial = false) + : U(Nd, grid), spatial_laplacian(spatial){}; - LaplacianAdjointField(GridBase* grid, OperatorFunction& S, LaplacianParams& p, const RealD k = 1.0) - : U(Nd, grid), Solver(S), param(p), kappa(k){ - AlgRemez remez(param.lo,param.hi,param.precision); - std::cout<(_U, mu); - } } - void M(const GaugeField& in, GaugeField& out) { - // in is an antihermitian matrix - // test - //GaugeField herm = in + adj(in); - //std::cout << "AHermiticity: " << norm2(herm) << std::endl; + void M(const FermionField &in, FermionField &out) + { + int dims = spatial_laplacian ? (Nd - 1) : Nd; - GaugeLinkField tmp(in._grid); - GaugeLinkField tmp2(in._grid); - GaugeLinkField sum(in._grid); - - for (int nu = 0; nu < Nd; nu++) { - sum = zero; - GaugeLinkField in_nu = PeekIndex(in, nu); - GaugeLinkField out_nu(out._grid); - for (int mu = 0; mu < Nd; mu++) { - tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]); - tmp2 = adj(U[mu]) * in_nu * U[mu]; - sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu; - } - out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum; - PokeIndex(out, out_nu, nu); - } + out = -2.0 * dims * in; + // eventually speed up with the stencil operator, if necessary + for (int mu = 0; mu < dims; mu++) + out += Impl::CovShiftForward(U[mu], mu, in) + Impl::CovShiftBackward(U[mu], mu, in); } - void MDeriv(const GaugeField& in, GaugeField& der) { - // in is anti-hermitian - RealD factor = -kappa / (double(4 * Nd)); - - for (int mu = 0; mu < Nd; mu++){ - GaugeLinkField der_mu(der._grid); - der_mu = zero; - for (int nu = 0; nu < Nd; nu++){ - GaugeLinkField in_nu = PeekIndex(in, nu); - der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu; - } - // the minus sign comes by using the in_nu instead of the - // adjoint in the last multiplication - PokeIndex(der, -2.0 * factor * der_mu, mu); - } - } - - // separating this temporarily - void MDeriv(const GaugeField& left, const GaugeField& right, - GaugeField& der) { - // in is anti-hermitian - RealD factor = -kappa / (double(4 * Nd)); - - for (int mu = 0; mu < Nd; mu++) { - GaugeLinkField der_mu(der._grid); - der_mu = zero; - for (int nu = 0; nu < Nd; nu++) { - GaugeLinkField left_nu = PeekIndex(left, nu); - GaugeLinkField right_nu = PeekIndex(right, nu); - der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu; - der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu; - } - PokeIndex(der, -factor * der_mu, mu); - } - } - - void Minv(const GaugeField& in, GaugeField& inverted){ - HermitianLinearOperator,GaugeField> HermOp(*this); - Solver(HermOp, in, inverted); - } - - void MSquareRoot(GaugeField& P){ - GaugeField Gp(P._grid); - HermitianLinearOperator,GaugeField> HermOp(*this); - ConjugateGradientMultiShift msCG(param.MaxIter,PowerHalf); - msCG(HermOp,P,Gp); - P = Gp; - } - - void MInvSquareRoot(GaugeField& P){ - GaugeField Gp(P._grid); - HermitianLinearOperator,GaugeField> HermOp(*this); - ConjugateGradientMultiShift msCG(param.MaxIter,PowerInvHalf); - msCG(HermOp,P,Gp); - P = Gp; - } - - - - private: - RealD kappa; +private: + bool spatial_laplacian; std::vector U; -}; - -} -} +}; // Laplacian +} // QCD +} // Grid #endif diff --git a/tests/core/Test_laplacian.cc b/tests/core/Test_laplacian.cc new file mode 100644 index 00000000..67c1d1d4 --- /dev/null +++ b/tests/core/Test_laplacian.cc @@ -0,0 +1,105 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_laplacian.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 + +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(&Grid); + + int threads = GridThread::GetThreads(); + std::cout<({45,12,81,9})); + + + std::vector point({0,0,0,0}); + + LatticeFermion src (&Grid); //random(pRNG,src); + SpinColourVectorD Sp; + for (unsigned int s = 0; s < Ns; ++s) + for (unsigned int c = 0; c < Nc; ++c) + Sp()(s)(c) = 1; + + src = zero; + pokeSite(Sp,src,point); + + LatticeFermion result(&Grid); result=zero; + LatticeFermion tmp(&Grid); tmp=zero; + + // Gauge configuration + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + std::cout< LaplaceOperator(src._grid); + LaplaceOperator.ImportGauge(Umu); + LaplaceOperator.M(src, result); + + std::cout << "Source vector" << std::endl; + std::cout << src << std::endl; + + std::cout << "Result vector" << std::endl; + std::cout << result << std::endl; + + std::cout< + + 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. + *******************************************************************************/ + +#include + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.01, .04, .2 , .25 , .3 }; + unsigned int nt = GridDefaultLatt()[Tp]; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + globalPar.genetic.maxGen = 1000; + globalPar.genetic.maxCstGen = 200; + globalPar.genetic.popSize = 20; + globalPar.genetic.mutationRate = .1; + application.setPar(globalPar); + + + // gauge field + application.createModule("gauge"); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + } + for (unsigned int t = 0; t < nt; t += 1) + { + std::string srcName; + std::string lapName; + std::vector qName; + std::vector> seqName; + + // Z2 source + MSource::Z2::Par z2Par; + z2Par.tA = t; + z2Par.tB = t; + srcName = "z2_" + std::to_string(t); + application.createModule(srcName, z2Par); + + // Example of smearing of the source + MSource::Laplacian::Par LapPar; + LapPar.N = 10; + LapPar.alpha = 0.1; + LapPar.source = srcName; + LapPar.gauge = "gauge"; + lapName = "z2smr_" + std::to_string(t); + application.createModule(lapName, LapPar); + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // sequential sources + MSource::SeqGamma::Par seqPar; + qName.push_back("QZ2_" + flavour[i] + "_" + std::to_string(t)); + seqPar.q = qName[i]; + seqPar.tA = (t + nt/4) % nt; + seqPar.tB = (t + nt/4) % nt; + seqPar.mom = "1. 0. 0. 0."; + seqName.push_back(std::vector(Nd)); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + seqPar.gamma = 0x1 << mu; + seqName[i][mu] = "G" + std::to_string(seqPar.gamma) + + "_" + std::to_string(seqPar.tA) + "-" + + qName[i]; + application.createModule(seqName[i][mu], seqPar); + } + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = srcName; + application.createModule(qName[i], quarkPar); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + quarkPar.source = seqName[i][mu]; + seqName[i][mu] = "Q_" + flavour[i] + "-" + seqName[i][mu]; + application.createModule(seqName[i][mu], quarkPar); + } + } + + // contractions + MContraction::Meson::Par mesPar; + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = qName[i]; + mesPar.q2 = qName[j]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("meson_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j], + mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = 0; j < flavour.size(); ++j) + for (unsigned int mu = 0; mu < Nd; ++mu) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_" + + std::to_string(mu); + mesPar.q1 = qName[i]; + mesPar.q2 = seqName[j][mu]; + mesPar.gammas = "all"; + mesPar.sink = "sink"; + application.createModule("3pt_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j] + + "_" + + std::to_string(mu), + mesPar); + } + } + + // execution + application.saveParameterFile("meson3pt.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}