portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Hadrons: templatisation of the fermion implementation

Browse Source
This commit is contained in:
Antonin Portelli 2016-12-05 16:47:29 +09:00
parent 7ae734103e
commit d1f7c6b94e
18 changed files with 483 additions and 638 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
extras/Hadrons/Global.hpp
View File

@ -56,10 +56,19 @@ using Grid::operator<<;
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
// type aliases // type aliases
typedef FermionOperator<FIMPL> FMat; //typedef FermionOperator<FIMPL> FMat;
typedef FIMPL::FermionField FermionField; //typedef FIMPL::FermionField FermionField;
typedef FIMPL::PropagatorField PropagatorField; //typedef FIMPL::PropagatorField PropagatorField;
typedef std::function<void(FermionField &, const FermionField &)> SolverFn; //typedef std::function<void(FermionField &, const FermionField &)> SolverFn;
#define TYPE_ALIASES(FImpl, suffix)\
typedef FermionOperator<FImpl> FMat##suffix; \
typedef typename FImpl::FermionField FermionField##suffix; \
typedef typename FImpl::PropagatorField PropagatorField##suffix; \
typedef typename FImpl::SitePropagator SitePropagator##suffix; \
typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\
typedef std::function<void(FermionField##suffix &, \
const FermionField##suffix &)> SolverFn##suffix;
// logger // logger
class HadronsLogger: public Logger class HadronsLogger: public Logger

2
extras/Hadrons/Modules.hpp
View File

@ -4,7 +4,7 @@
#include <Grid/Hadrons/Modules/MGauge/Load.hpp> #include <Grid/Hadrons/Modules/MGauge/Load.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp> #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp> #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MQuark.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp> #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp> #include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp> #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/Quark.hpp>

82
extras/Hadrons/Modules/MAction/DWF.cc
View File

@ -1,82 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/DWF.cc
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MAction/DWF.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MAction;
/******************************************************************************
* DWF implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
DWF::DWF(const std::string name)
: Module<DWFPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> DWF::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
std::vector<std::string> DWF::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void DWF::setup(void)
{
unsigned int size;
size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
env().registerObject(getName(), size, par().Ls);
}
// execution ///////////////////////////////////////////////////////////////////
void DWF::execute(void)
{
LOG(Message) << "Setting up domain wall fermion matrix with m= "
<< par().mass << ", M5= " << par().M5 << " and Ls= "
<< par().Ls << " using gauge field '" << par().gauge << "'"
<< std::endl;
env().createGrid(par().Ls);
auto &U = *env().getObject<LatticeGaugeField>(par().gauge);
auto &g4 = *env().getGrid();
auto &grb4 = *env().getRbGrid();
auto &g5 = *env().getGrid(par().Ls);
auto &grb5 = *env().getRbGrid(par().Ls);
FMat *fMatPt = new DomainWallFermion<FIMPL>(U, g5, grb5, g4, grb4,
par().mass, par().M5);
env().setObject(getName(), fMatPt);
}

66
extras/Hadrons/Modules/MAction/DWF.hpp
View File

@ -49,13 +49,16 @@ public:
double , M5); double , M5);
}; };
class DWF: public Module<DWFPar> template <typename FImpl>
class TDWF: public Module<DWFPar>
{ {
public:
TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
DWF(const std::string name); TDWF(const std::string name);
// destructor // destructor
virtual ~DWF(void) = default; virtual ~TDWF(void) = default;
// dependency relation // dependency relation
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -65,6 +68,63 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* DWF template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TDWF<FImpl>::TDWF(const std::string name)
: Module<DWFPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TDWF<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TDWF<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDWF<FImpl>::setup(void)
{
unsigned int size;
size = 3*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
env().registerObject(getName(), size, par().Ls);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TDWF<FImpl>::execute(void)
{
LOG(Message) << "Setting up domain wall fermion matrix with m= "
<< par().mass << ", M5= " << par().M5 << " and Ls= "
<< par().Ls << " using gauge field '" << par().gauge << "'"
<< std::endl;
env().createGrid(par().Ls);
auto &U = *env().template getObject<LatticeGaugeField>(par().gauge);
auto &g4 = *env().getGrid();
auto &grb4 = *env().getRbGrid();
auto &g5 = *env().getGrid(par().Ls);
auto &grb5 = *env().getRbGrid(par().Ls);
FMat *fMatPt = new DomainWallFermion<FImpl>(U, g5, grb5, g4, grb4,
par().mass, par().M5);
env().setObject(getName(), fMatPt);
}
typedef TDWF<FIMPL> DWF;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(DWF, MAction); MODULE_REGISTER_NS(DWF, MAction);

76
extras/Hadrons/Modules/MAction/Wilson.cc
View File

@ -1,76 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Wilson.cc
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MAction/Wilson.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MAction;
/******************************************************************************
* Wilson implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
Wilson::Wilson(const std::string name)
: Module<WilsonPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> Wilson::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
std::vector<std::string> Wilson::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void Wilson::setup(void)
{
unsigned int size;
size = 3*env().lattice4dSize<WilsonFermionR::DoubledGaugeField>();
env().registerObject(getName(), size);
}
// execution ///////////////////////////////////////////////////////////////////
void Wilson::execute()
{
LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
<< " using gauge field '" << par().gauge << "'" << std::endl;
auto &U = *env().getObject<LatticeGaugeField>(par().gauge);
auto &grid = *env().getGrid();
auto &gridRb = *env().getRbGrid();
FMat *fMatPt = new WilsonFermion<FIMPL>(U, grid, gridRb, par().mass);
env().setObject(getName(), fMatPt);
}

64
extras/Hadrons/Modules/MAction/Wilson.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Wilson.hpp Source file: programs/Hadrons/TWilson.hpp
Copyright (C) 2016 Copyright (C) 2016
@ -35,7 +35,7 @@ directory.
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
/****************************************************************************** /******************************************************************************
* Wilson quark action * * TWilson quark action *
******************************************************************************/ ******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction) BEGIN_MODULE_NAMESPACE(MAction)
@ -47,13 +47,16 @@ public:
double , mass); double , mass);
}; };
class Wilson: public Module<WilsonPar> template <typename FImpl>
class TWilson: public Module<WilsonPar>
{ {
public:
TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
Wilson(const std::string name); TWilson(const std::string name);
// destructor // destructor
virtual ~Wilson(void) = default; virtual ~TWilson(void) = default;
// dependencies/products // dependencies/products
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -63,6 +66,57 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* TWilson template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TWilson<FImpl>::TWilson(const std::string name)
: Module<WilsonPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TWilson<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TWilson<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWilson<FImpl>::setup(void)
{
unsigned int size;
size = 3*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
env().registerObject(getName(), size);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWilson<FImpl>::execute()
{
LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
<< " using gauge field '" << par().gauge << "'" << std::endl;
auto &U = *env().template getObject<LatticeGaugeField>(par().gauge);
auto &grid = *env().getGrid();
auto &gridRb = *env().getRbGrid();
FMat *fMatPt = new WilsonFermion<FImpl>(U, grid, gridRb, par().mass);
env().setObject(getName(), fMatPt);
}
typedef TWilson<FIMPL> Wilson;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(Wilson, MAction); MODULE_REGISTER_NS(Wilson, MAction);

94
extras/Hadrons/Modules/MContraction/Meson.cc
View File

@ -1,94 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Meson.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MContraction/Meson.hpp>
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
using namespace MContraction;
/******************************************************************************
* Meson implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
Meson::Meson(const std::string name)
: Module<MesonPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> Meson::getInput(void)
{
std::vector<std::string> input = {par().q1, par().q2};
return input;
}
std::vector<std::string> Meson::getOutput(void)
{
std::vector<std::string> output = {getName()};
return output;
}
// execution ///////////////////////////////////////////////////////////////////
void Meson::execute(void)
{
LOG(Message) << "Computing meson contraction '" << getName() << "' using"
<< " quarks '" << par().q1 << "' and '" << par().q2 << "'"
<< std::endl;
XmlWriter writer(par().output);
PropagatorField &q1 = *env().getObject<PropagatorField>(par().q1);
PropagatorField &q2 = *env().getObject<PropagatorField>(par().q2);
LatticeComplex c(env().getGrid());
SpinMatrix g[Ns*Ns], g5;
std::vector<TComplex> buf;
Result result;
g5 = makeGammaProd(Ns*Ns - 1);
result.corr.resize(Ns*Ns);
for (unsigned int i = 0; i < Ns*Ns; ++i)
{
g[i] = makeGammaProd(i);
}
for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
{
result.corr[iSink].resize(Ns*Ns);
for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
{
c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
sliceSum(c, buf, Tp);
result.corr[iSink][iSrc].resize(buf.size());
for (unsigned int t = 0; t < buf.size(); ++t)
{
result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
}
}
}
write(writer, "meson", result);
}

81
extras/Hadrons/Modules/MContraction/Meson.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Meson.hpp Source file: programs/Hadrons/TMeson.hpp
Copyright (C) 2015 Copyright (C) 2015
@ -35,7 +35,7 @@ directory.
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
/****************************************************************************** /******************************************************************************
* Meson * * TMeson *
******************************************************************************/ ******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction) BEGIN_MODULE_NAMESPACE(MContraction)
@ -48,8 +48,12 @@ public:
std::string, output); std::string, output);
}; };
class Meson: public Module<MesonPar> template <typename FImpl1, typename FImpl2>
class TMeson: public Module<MesonPar>
{ {
public:
TYPE_ALIASES(FImpl1, 1);
TYPE_ALIASES(FImpl2, 2);
public: public:
class Result: Serializable class Result: Serializable
{ {
@ -59,9 +63,9 @@ public:
}; };
public: public:
// constructor // constructor
Meson(const std::string name); TMeson(const std::string name);
// destructor // destructor
virtual ~Meson(void) = default; virtual ~TMeson(void) = default;
// dependencies/products // dependencies/products
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -69,6 +73,73 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* TMeson implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl1, typename FImpl2>
TMeson<FImpl1, FImpl2>::TMeson(const std::string name)
: Module<MesonPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl1, typename FImpl2>
std::vector<std::string> TMeson<FImpl1, FImpl2>::getInput(void)
{
std::vector<std::string> input = {par().q1, par().q2};
return input;
}
template <typename FImpl1, typename FImpl2>
std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
{
std::vector<std::string> output = {getName()};
return output;
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl1, typename FImpl2>
void TMeson<FImpl1, FImpl2>::execute(void)
{
LOG(Message) << "Computing meson contraction '" << getName() << "' using"
<< " quarks '" << par().q1 << "' and '" << par().q2 << "'"
<< std::endl;
XmlWriter writer(par().output);
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
LatticeComplex c(env().getGrid());
SpinMatrix g[Ns*Ns], g5;
std::vector<TComplex> buf;
Result result;
g5 = makeGammaProd(Ns*Ns - 1);
result.corr.resize(Ns*Ns);
for (unsigned int i = 0; i < Ns*Ns; ++i)
{
g[i] = makeGammaProd(i);
}
for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
{
result.corr[iSink].resize(Ns*Ns);
for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
{
c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
sliceSum(c, buf, Tp);
result.corr[iSink][iSrc].resize(buf.size());
for (unsigned int t = 0; t < buf.size(); ++t)
{
result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
}
}
}
write(writer, "meson", result);
}
typedef TMeson<FIMPL, FIMPL> Meson;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(Meson, MContraction); MODULE_REGISTER_NS(Meson, MContraction);

71
extras/Hadrons/Modules/MQuark.hpp
View File

@ -1,71 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/MQuark.hpp
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
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.
*******************************************************************************/
#ifndef Hadrons_MQuark_hpp_
#define Hadrons_MQuark_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* MQuark *
******************************************************************************/
class MQuarkPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(MQuarkPar,
std::string, source,
std::string, solver);
};
class MQuark: public Module<MQuarkPar>
{
public:
// constructor
MQuark(const std::string name);
// destructor
virtual ~MQuark(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
unsigned int Ls_;
SolverFn *solver_{nullptr};
};
MODULE_REGISTER(MQuark);
END_HADRONS_NAMESPACE
#endif // Hadrons_MQuark_hpp_

84
extras/Hadrons/Modules/MSolver/RBPrecCG.cc
View File

@ -1,84 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/RBPrecCG.cc
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
using namespace MSolver;
/******************************************************************************
* RBPrecCG implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
RBPrecCG::RBPrecCG(const std::string name)
: Module(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> RBPrecCG::getInput(void)
{
std::vector<std::string> in = {par().action};
return in;
}
std::vector<std::string> RBPrecCG::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void RBPrecCG::setup(void)
{
auto Ls = env().getObjectLs(par().action);
env().registerObject(getName(), 0, Ls);
env().addOwnership(getName(), par().action);
}
// execution ///////////////////////////////////////////////////////////////////
void RBPrecCG::execute(void)
{
auto &mat = *(env().getObject<FMat>(par().action));
auto solver = [&mat, this](LatticeFermion &sol,
const LatticeFermion &source)
{
ConjugateGradient<LatticeFermion> cg(par().residual, 10000);
SchurRedBlackDiagMooeeSolve<LatticeFermion> schurSolver(cg);
schurSolver(mat, source, sol);
};
LOG(Message) << "setting up Schur red-black preconditioned CG for"
<< " action '" << par().action << "' with residual "
<< par().residual << std::endl;
env().setObject(getName(), new SolverFn(solver));
}

68
extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/RBPrecCG.hpp Source file: programs/Hadrons/TRBPrecCG.hpp
Copyright (C) 2016 Copyright (C) 2016
@ -47,13 +47,16 @@ public:
double , residual); double , residual);
}; };
class RBPrecCG: public Module<RBPrecCGPar> template <typename FImpl>
class TRBPrecCG: public Module<RBPrecCGPar>
{ {
public:
TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
RBPrecCG(const std::string name); TRBPrecCG(const std::string name);
// destructor // destructor
virtual ~RBPrecCG(void) = default; virtual ~TRBPrecCG(void) = default;
// dependencies/products // dependencies/products
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -63,6 +66,63 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* TRBPrecCG template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TRBPrecCG<FImpl>::TRBPrecCG(const std::string name)
: Module(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TRBPrecCG<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().action};
return in;
}
template <typename FImpl>
std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TRBPrecCG<FImpl>::setup(void)
{
auto Ls = env().getObjectLs(par().action);
env().registerObject(getName(), 0, Ls);
env().addOwnership(getName(), par().action);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TRBPrecCG<FImpl>::execute(void)
{
auto &mat = *(env().template getObject<FMat>(par().action));
auto solver = [&mat, this](FermionField &sol, const FermionField &source)
{
ConjugateGradient<FermionField> cg(par().residual, 10000);
SchurRedBlackDiagMooeeSolve<FermionField> schurSolver(cg);
schurSolver(mat, source, sol);
};
LOG(Message) << "setting up Schur red-black preconditioned CG for"
<< " action '" << par().action << "' with residual "
<< par().residual << std::endl;
env().setObject(getName(), new SolverFn(solver));
}
typedef TRBPrecCG<FIMPL> RBPrecCG;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(RBPrecCG, MSolver); MODULE_REGISTER_NS(RBPrecCG, MSolver);

75
extras/Hadrons/Modules/MSource/Point.cc
View File

@ -1,75 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Point.cc
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MSource/Point.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MSource;
/******************************************************************************
* Point implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
Point::Point(const std::string name)
: Module<PointPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> Point::getInput(void)
{
std::vector<std::string> in;
return in;
}
std::vector<std::string> Point::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void Point::setup(void)
{
env().registerLattice<PropagatorField>(getName());
}
// execution ///////////////////////////////////////////////////////////////////
void Point::execute(void)
{
std::vector<int> position = strToVec<int>(par().position);
SpinColourMatrix id;
LOG(Message) << "Creating point source at position [" << par().position
<< "]" << std::endl;
PropagatorField &src = *env().createLattice<PropagatorField>(getName());
id = 1.;
src = zero;
pokeSite(id, src, position);
}

63
extras/Hadrons/Modules/MSource/Point.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Point.hpp Source file: programs/Hadrons/TPoint.hpp
Copyright (C) 2016 Copyright (C) 2016
@ -46,7 +46,7 @@ BEGIN_HADRONS_NAMESPACE
*/ */
/****************************************************************************** /******************************************************************************
* Point * * TPoint *
******************************************************************************/ ******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSource) BEGIN_MODULE_NAMESPACE(MSource)
@ -57,13 +57,16 @@ public:
std::string, position); std::string, position);
}; };
class Point: public Module<PointPar> template <typename FImpl>
class TPoint: public Module<PointPar>
{ {
public:
TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
Point(const std::string name); TPoint(const std::string name);
// destructor // destructor
virtual ~Point(void) = default; virtual ~TPoint(void) = default;
// dependency relation // dependency relation
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -73,6 +76,56 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* TPoint template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TPoint<FImpl>::TPoint(const std::string name)
: Module<PointPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TPoint<FImpl>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename FImpl>
std::vector<std::string> TPoint<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPoint<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TPoint<FImpl>::execute(void)
{
std::vector<int> position = strToVec<int>(par().position);
typename SitePropagator::scalar_object id;
LOG(Message) << "Creating point source at position [" << par().position
<< "]" << std::endl;
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
id = 1.;
src = zero;
pokeSite(id, src, position);
}
typedef TPoint<FIMPL> Point;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(Point, MSource); MODULE_REGISTER_NS(Point, MSource);

88
extras/Hadrons/Modules/MSource/Z2.cc
View File

@ -1,88 +0,0 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Z2.cc
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Grid/Hadrons/Modules/MSource/Z2.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MSource;
/******************************************************************************
* Z2 implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
Z2::Z2(const std::string name)
: Module<Z2Par>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> Z2::getInput(void)
{
std::vector<std::string> in;
return in;
}
std::vector<std::string> Z2::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void Z2::setup(void)
{
env().registerLattice<PropagatorField>(getName());
}
// execution ///////////////////////////////////////////////////////////////////
void Z2::execute(void)
{
Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex eta(env().getGrid());
LatticeFermion phi(env().getGrid());
Complex shift(1., 1.);
if (par().tA == par().tB)
{
LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
<< std::endl;
}
else
{
LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
<< par().tB << std::endl;
}
PropagatorField &src = *env().createLattice<PropagatorField>(getName());
LatticeCoordinate(t, Tp);
bernoulli(*env().get4dRng(), eta);
eta = (2.*eta - shift)*(1./::sqrt(2.));
eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);
src = 1.;
src = src*eta;
}

75
extras/Hadrons/Modules/MSource/Z2.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/Z2.hpp Source file: programs/Hadrons/TZ2.hpp
Copyright (C) 2016 Copyright (C) 2016
@ -47,7 +47,7 @@ BEGIN_HADRONS_NAMESPACE
*/ */
/****************************************************************************** /******************************************************************************
* Z2 * * Z2 stochastic source *
******************************************************************************/ ******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSource) BEGIN_MODULE_NAMESPACE(MSource)
@ -59,13 +59,16 @@ public:
unsigned int, tB); unsigned int, tB);
}; };
class Z2: public Module<Z2Par> template <typename FImpl>
class TZ2: public Module<Z2Par>
{ {
public:
TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
Z2(const std::string name); TZ2(const std::string name);
// destructor // destructor
virtual ~Z2(void) = default; virtual ~TZ2(void) = default;
// dependency relation // dependency relation
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
@ -75,6 +78,68 @@ public:
virtual void execute(void); virtual void execute(void);
}; };
/******************************************************************************
* TZ2 template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TZ2<FImpl>::TZ2(const std::string name)
: Module<Z2Par>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TZ2<FImpl>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <typename FImpl>
std::vector<std::string> TZ2<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TZ2<FImpl>::setup(void)
{
env().template registerLattice<PropagatorField>(getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TZ2<FImpl>::execute(void)
{
Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex eta(env().getGrid());
Complex shift(1., 1.);
if (par().tA == par().tB)
{
LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
<< std::endl;
}
else
{
LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
<< par().tB << std::endl;
}
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
LatticeCoordinate(t, Tp);
bernoulli(*env().get4dRng(), eta);
eta = (2.*eta - shift)*(1./::sqrt(2.));
eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);
src = 1.;
src = src*eta;
}
typedef TZ2<FIMPL> Z2;
END_MODULE_NAMESPACE END_MODULE_NAMESPACE
MODULE_REGISTER_NS(Z2, MSource); MODULE_REGISTER_NS(Z2, MSource);

94
extras/Hadrons/Modules/MQuark.cc → extras/Hadrons/Modules/Quark.hpp
View File

@ -1,7 +1,7 @@
/******************************************************************************* /*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: programs/Hadrons/MQuark.cc Source file: programs/Hadrons/TQuark.hpp
Copyright (C) 2015 Copyright (C) 2015
@ -25,29 +25,68 @@ See the full license in the file "LICENSE" in the top level distribution
directory. directory.
*******************************************************************************/ *******************************************************************************/
#include <Grid/Hadrons/Modules/MQuark.hpp> #ifndef Hadrons_Quark_hpp_
#define Hadrons_Quark_hpp_
using namespace Grid; #include <Grid/Hadrons/Global.hpp>
using namespace QCD; #include <Grid/Hadrons/Module.hpp>
using namespace Hadrons; #include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/****************************************************************************** /******************************************************************************
* MQuark implementation * * TQuark *
******************************************************************************/
class QuarkPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(QuarkPar,
std::string, source,
std::string, solver);
};
template <typename FImpl>
class TQuark: public Module<QuarkPar>
{
public:
TYPE_ALIASES(FImpl,);
public:
// constructor
TQuark(const std::string name);
// destructor
virtual ~TQuark(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
unsigned int Ls_;
SolverFn *solver_{nullptr};
};
/******************************************************************************
* TQuark implementation *
******************************************************************************/ ******************************************************************************/
// constructor ///////////////////////////////////////////////////////////////// // constructor /////////////////////////////////////////////////////////////////
MQuark::MQuark(const std::string name) template <typename FImpl>
TQuark<FImpl>::TQuark(const std::string name)
: Module(name) : Module(name)
{} {}
// dependencies/products /////////////////////////////////////////////////////// // dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> MQuark::getInput(void) template <typename FImpl>
std::vector<std::string> TQuark<FImpl>::getInput(void)
{ {
std::vector<std::string> in = {par().source, par().solver}; std::vector<std::string> in = {par().source, par().solver};
return in; return in;
} }
std::vector<std::string> MQuark::getOutput(void) template <typename FImpl>
std::vector<std::string> TQuark<FImpl>::getOutput(void)
{ {
std::vector<std::string> out = {getName(), getName() + "_5d"}; std::vector<std::string> out = {getName(), getName() + "_5d"};
@ -55,31 +94,33 @@ std::vector<std::string> MQuark::getOutput(void)
} }
// setup /////////////////////////////////////////////////////////////////////// // setup ///////////////////////////////////////////////////////////////////////
void MQuark::setup(void) template <typename FImpl>
void TQuark<FImpl>::setup(void)
{ {
Ls_ = env().getObjectLs(par().solver); Ls_ = env().getObjectLs(par().solver);
env().registerLattice<LatticePropagator>(getName()); env().template registerLattice<PropagatorField>(getName());
if (Ls_ > 1) if (Ls_ > 1)
{ {
env().registerLattice<LatticePropagator>(getName() + "_5d", Ls_); env().template registerLattice<PropagatorField>(getName() + "_5d", Ls_);
} }
} }
// execution /////////////////////////////////////////////////////////////////// // execution ///////////////////////////////////////////////////////////////////
void MQuark::execute(void) template <typename FImpl>
void TQuark<FImpl>::execute(void)
{ {
LatticeFermion source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
tmp(env().getGrid());
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
LOG(Message) << "Computing quark propagator '" << getName() << "'" LOG(Message) << "Computing quark propagator '" << getName() << "'"
<< std::endl; << std::endl;
LatticePropagator &prop = *env().createLattice<LatticePropagator>(propName);
LatticePropagator &fullSrc = *env().getObject<LatticePropagator>(par().source); FermionField source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
SolverFn &solver = *env().getObject<SolverFn>(par().solver); tmp(env().getGrid());
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
PropagatorField &prop = *env().template createLattice<PropagatorField>(propName);
PropagatorField &fullSrc = *env().template getObject<PropagatorField>(par().source);
SolverFn &solver = *env().template getObject<SolverFn>(par().solver);
if (Ls_ > 1) if (Ls_ > 1)
{ {
env().createLattice<LatticePropagator>(getName()); env().template createLattice<PropagatorField>(getName());
} }
LOG(Message) << "Inverting using solver '" << par().solver LOG(Message) << "Inverting using solver '" << par().solver
@ -124,7 +165,8 @@ void MQuark::execute(void)
// create 4D propagators from 5D one if necessary // create 4D propagators from 5D one if necessary
if (Ls_ > 1) if (Ls_ > 1)
{ {
LatticePropagator &p4d = *env().getObject<LatticePropagator>(getName()); PropagatorField &p4d =
*env().template getObject<PropagatorField>(getName());
axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0); axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1); axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1);
@ -133,3 +175,11 @@ void MQuark::execute(void)
} }
} }
} }
typedef TQuark<FIMPL> Quark;
MODULE_REGISTER(Quark);
END_HADRONS_NAMESPACE
#endif // Hadrons_Quark_hpp_

13
extras/Hadrons/modules.inc
View File

@ -1,14 +1,7 @@
modules_cc =\ modules_cc =\
Modules/MAction/DWF.cc \
Modules/MAction/Wilson.cc \
Modules/MContraction/Meson.cc \
Modules/MGauge/Load.cc \ Modules/MGauge/Load.cc \
Modules/MGauge/Random.cc \ Modules/MGauge/Random.cc \
Modules/MGauge/Unit.cc \ Modules/MGauge/Unit.cc
Modules/MQuark.cc \
Modules/MSolver/RBPrecCG.cc \
Modules/MSource/Point.cc \
Modules/MSource/Z2.cc
modules_hpp =\ modules_hpp =\
Modules/MAction/DWF.hpp \ Modules/MAction/DWF.hpp \
@ -17,8 +10,8 @@ modules_hpp =\
Modules/MGauge/Load.hpp \ Modules/MGauge/Load.hpp \
Modules/MGauge/Random.hpp \ Modules/MGauge/Random.hpp \
Modules/MGauge/Unit.hpp \ Modules/MGauge/Unit.hpp \
Modules/MQuark.hpp \
Modules/MSolver/RBPrecCG.hpp \ Modules/MSolver/RBPrecCG.hpp \
Modules/MSource/Point.hpp \ Modules/MSource/Point.hpp \
Modules/MSource/Z2.hpp Modules/MSource/Z2.hpp \
Modules/Quark.hpp

4
tests/hadrons/Test_hadrons_mesons.cc
View File

@ -82,10 +82,10 @@ int main(int argc, char *argv[])
solverPar); solverPar);
// propagators // propagators
MQuark::Par quarkPar; Quark::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i]; quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = (flavour[i][0] == 'c') ? "z2" : "pt"; quarkPar.source = (flavour[i][0] == 'c') ? "z2" : "pt";
application.createModule<MQuark>("Q_" + flavour[i], quarkPar); application.createModule<Quark>("Q_" + flavour[i], quarkPar);
} }
for (unsigned int i = 0; i < flavour.size(); ++i) for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i + 1; j < flavour.size(); ++j) for (unsigned int j = i + 1; j < flavour.size(); ++j)