portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Added TwoFlavorsEO

Browse Source 
Had to remove a conformability check in the Derivative of SchurDiff,
see the comments in the file
This commit is contained in:
Guido Cossu
2017-01-20 16:59:31 +00:00
parent f96fac0aee
commit 27dfe816fa
9 changed files with 276 additions and 191 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -247,8 +247,9 @@ template <class Impl>
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
conformable(U._grid, _cbgrid);
conformable(U._grid, V._grid);
conformable(U._grid, mat._grid);
//conformable(U._grid, mat._grid); not general, leaving as a comment (Guido)
// Motivation: look at the SchurDiff operator
assert(V.checkerboard == Even);
assert(U.checkerboard == Odd);
mat.checkerboard = Odd;
@ -260,7 +261,7 @@ template <class Impl>
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
conformable(U._grid, _cbgrid);
conformable(U._grid, V._grid);
conformable(U._grid, mat._grid);
//conformable(U._grid, mat._grid);
assert(V.checkerboard == Odd);
assert(U.checkerboard == Even);

4
lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h
View File

@ -67,6 +67,7 @@ namespace Grid{
// NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE
// it is not conformable with the HMC force field
// Case: Ls vectorised fields
// INHERIT FROM THE Force field instead
GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid);
GaugeField ForceO(forcecb);
@ -76,9 +77,8 @@ namespace Grid{
// X^dag Der_oe MeeInv Meo Y
// Use Mooee as nontrivial but gauge field indept
this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied
this->_Mat.MooeeInv(tmp1,tmp2); // even->even
this->_Mat.MooeeInv(tmp1,tmp2); // even->even
this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo);
// Accumulate X^dag M_oe MeeInv Der_eo Y
this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied
this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even

108
lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h
View File

@ -31,82 +31,89 @@ directory
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
namespace Grid {
namespace QCD {
namespace QCD {
////////////////////////////////////////////////////////////////////////
// Two flavour pseudofermion action for any EO prec dop
////////////////////////////////////////////////////////////////////////
template <class Impl>
class TwoFlavourEvenOddPseudoFermionAction
: public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
////////////////////////////////////////////////////////////////////////
// Two flavour pseudofermion action for any EO prec dop
////////////////////////////////////////////////////////////////////////
template <class Impl>
class TwoFlavourEvenOddPseudoFermionAction
: public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
private:
FermionOperator<Impl> &FermOp; // the basic operator
private:
FermionOperator<Impl> &FermOp; // the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
FermionField PhiOdd; // the pseudo fermion field for this trajectory
FermionField PhiEven; // the pseudo fermion field for this trajectory
FermionField PhiOdd; // the pseudo fermion field for this trajectory
FermionField PhiEven; // the pseudo fermion field for this trajectory
public:
/////////////////////////////////////////////////
// Pass in required objects.
/////////////////////////////////////////////////
TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op,
OperatorFunction<FermionField> &DS,
OperatorFunction<FermionField> &AS)
: FermOp(Op),
DerivativeSolver(DS),
ActionSolver(AS),
PhiEven(Op.FermionRedBlackGrid()),
PhiOdd(Op.FermionRedBlackGrid())
{};
public:
/////////////////////////////////////////////////
// Pass in required objects.
/////////////////////////////////////////////////
TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op,
OperatorFunction<FermionField> &DS,
OperatorFunction<FermionField> &AS)
: FermOp(Op),
DerivativeSolver(DS),
ActionSolver(AS),
PhiEven(Op.FermionRedBlackGrid()),
PhiOdd(Op.FermionRedBlackGrid())
{};
virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";}
virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
return sstream.str();
}
//////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied
//////////////////////////////////////////////////////////////////////////////////////
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
// Phi = McpDag eta
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
RealD scale = std::sqrt(0.5);
FermionField eta (FermOp.FermionGrid());
FermionField etaOdd (FermOp.FermionRedBlackGrid());
FermionField etaEven(FermOp.FermionRedBlackGrid());
gaussian(pRNG,eta);
pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta);
FermOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> PCop(FermOp);
PCop.MpcDag(etaOdd,PhiOdd);
FermOp.MooeeDag(etaEven,PhiEven);
PhiOdd =PhiOdd*scale;
PhiEven=PhiEven*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag (Mdag M)^-1 phi (odd)
// + phi^dag (Mdag M)^-1 phi (even)
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid());
@ -137,7 +144,7 @@ class TwoFlavourEvenOddPseudoFermionAction
//
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
std::cout << GridLogDebug << "Calling deriv" << std::endl;
FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid());
@ -151,10 +158,17 @@ class TwoFlavourEvenOddPseudoFermionAction
X=zero;
DerivativeSolver(Mpc,PhiOdd,X);
std::cout << GridLogDebug << "Calling deriv 2 " << std::endl;
Mpc.Mpc(X,Y);
Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp;
std::cout << GridLogDebug << "Calling deriv 3 " << std::endl;
Mpc.MpcDeriv(tmp , Y, X );
std::cout << GridLogDebug << "Calling deriv 4 " << std::endl;
dSdU=tmp;
std::cout << GridLogDebug << "Calling deriv 5 " << std::endl;
Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp;
std::cout << GridLogDebug << "Calling deriv 6" << std::endl;
// Treat the EE case. (MdagM)^-1 = Minv Minvdag
// Deriv defaults to zero.
// FermOp.MooeeInvDag(PhiOdd,Y);
@ -165,10 +179,10 @@ class TwoFlavourEvenOddPseudoFermionAction
assert(FermOp.ConstEE() == 1);
/*
FermOp.MooeeInvDag(PhiOdd,Y);
FermOp.MooeeInv(Y,X);
FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
FermOp.MooeeInvDag(PhiOdd,Y);
FermOp.MooeeInv(Y,X);
FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
*/
//dSdU = Ta(dSdU);

8
lib/qcd/hmc/HMCModules.h
View File

@ -72,8 +72,12 @@ public:
class GridModule {
public:
GridCartesian* get_full() { return grid_.get(); }
GridRedBlackCartesian* get_rb() { return rbgrid_.get(); }
GridCartesian* get_full() {
std::cout << GridLogDebug << "Getting cartesian in module"<< std::endl;
return grid_.get(); }
GridRedBlackCartesian* get_rb() {
std::cout << GridLogDebug << "Getting rb-cartesian in module"<< std::endl;
return rbgrid_.get(); }
void set_full(GridCartesian* grid) { grid_.reset(grid); }
void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); }

107
lib/qcd/modules/ActionModules.h
View File

@ -40,11 +40,11 @@ namespace Grid {
// Actions
//////////////////////////////////////////////
template <class Product, class Resource>
template <class Product, class R>
class ActionModuleBase: public HMCModuleBase<Product>{
public:
typedef Resource Res;
virtual void acquireResource(Resource& ){};
typedef R Resource;
virtual void acquireResource(R& ){};
};
@ -193,57 +193,95 @@ class DBW2GModule: public ActionModule<DBW2GaugeAction<Impl>, BetaGaugeActionPar
/////////////////////////////////////////
template <class Impl >
class TwoFlavourFModule: public ActionModule<TwoFlavourPseudoFermionAction<Impl>, NoParameters> {
typedef ActionModule<TwoFlavourPseudoFermionAction<Impl>, NoParameters> ActionBase;
template <class Impl, template <typename> class FermionA >
class PseudoFermionModuleBase: public ActionModule<FermionA<Impl>, NoParameters> {
protected:
typedef ActionModule<FermionA<Impl>, NoParameters> ActionBase;
using ActionBase::ActionBase; // for constructors
std::unique_ptr<HMCModuleBase<OperatorFunction<typename Impl::FermionField> > >solver_mod;
std::unique_ptr<FermionOperatorModuleBase<FermionOperator<Impl>> > fop_mod;
public:
typedef std::unique_ptr<FermionOperatorModuleBase<FermionOperator<Impl>> > operator_type;
typedef std::unique_ptr<HMCModuleBase<OperatorFunction<typename Impl::FermionField> > > solver_type;
virtual void acquireResource(typename ActionBase::Res& GridMod){
fop_mod->AddGridPair(GridMod);
template <class ReaderClass>
void getFermionOperator(Reader<ReaderClass>& Reader, operator_type &fo, std::string section_name){
auto &FOFactory = HMC_FermionOperatorModuleFactory<fermionop_string, Impl, ReaderClass>::getInstance();
Reader.push(section_name);
std::string op_name;
read(Reader,"name", op_name);
fo = FOFactory.create(op_name, Reader);
Reader.pop();
}
template <class ReaderClass>
TwoFlavourFModule(Reader<ReaderClass>& Reader) : ActionBase(Reader){
std::cout << "Constructing TwoFlavourFModule" << std::endl;
void getSolverOperator(Reader<ReaderClass>& Reader, solver_type &so, std::string section_name){
auto& SolverFactory = HMC_SolverModuleFactory<solver_string, typename Impl::FermionField, ReaderClass>::getInstance();
Reader.push("Solver");
Reader.push(section_name);
std::string solv_name;
read(Reader,"name", solv_name);
solver_mod = SolverFactory.create(solv_name, Reader);
std::cout << "Registered types " << std::endl;
std::cout << SolverFactory.getBuilderList() << std::endl;
solver_mod->print_parameters();
Reader.pop();
so = SolverFactory.create(solv_name, Reader);
Reader.pop();
}
};
template <class Impl >
class TwoFlavourFModule: public PseudoFermionModuleBase<Impl, TwoFlavourPseudoFermionAction>{
typedef PseudoFermionModuleBase<Impl, TwoFlavourPseudoFermionAction> Base;
using Base::Base;
auto &FOFactory = HMC_FermionOperatorModuleFactory<fermionop_string, Impl, ReaderClass>::getInstance();
Reader.push("Operator");
std::string op_name;
read(Reader,"name", op_name);
fop_mod = FOFactory.create(op_name, Reader);
std::cout << "Registered types " << std::endl;
std::cout << FOFactory.getBuilderList() << std::endl;
typename Base::operator_type fop_mod;
typename Base::solver_type solver_mod;
fop_mod->print_parameters();
Reader.pop();
public:
virtual void acquireResource(typename Base::Resource& GridMod){
fop_mod->AddGridPair(GridMod);
}
// constructor
template <class ReaderClass>
TwoFlavourFModule(Reader<ReaderClass>& R): PseudoFermionModuleBase<Impl, TwoFlavourPseudoFermionAction>(R) {
this->getSolverOperator(R, solver_mod, "Solver");
this->getFermionOperator(R, fop_mod, "Operator");
R.pop();
}
};
private:
// acquire resource
virtual void initialize() {
this->ActionPtr.reset(new TwoFlavourPseudoFermionAction<Impl>(*fop_mod->getPtr(), *solver_mod->getPtr(), *solver_mod->getPtr()));
// here temporarily assuming that the force and action solver are the same
this->ActionPtr.reset(new TwoFlavourPseudoFermionAction<Impl>(*(this->fop_mod->getPtr()), *(this->solver_mod->getPtr()), *(this->solver_mod->getPtr())));
}
};
// very similar, I could have templated this but it is overkilling
template <class Impl >
class TwoFlavourEOFModule: public PseudoFermionModuleBase<Impl, TwoFlavourEvenOddPseudoFermionAction>{
typedef PseudoFermionModuleBase<Impl, TwoFlavourEvenOddPseudoFermionAction> Base;
using Base::Base;
typename Base::operator_type fop_mod;
typename Base::solver_type solver_mod;
public:
virtual void acquireResource(typename Base::Resource& GridMod){
fop_mod->AddGridPair(GridMod);
}
// constructor
template <class ReaderClass>
TwoFlavourEOFModule(Reader<ReaderClass>& R): PseudoFermionModuleBase<Impl, TwoFlavourEvenOddPseudoFermionAction>(R) {
this->getSolverOperator(R, solver_mod, "Solver");
this->getFermionOperator(R, fop_mod, "Operator");
R.pop();
}
// acquire resource
virtual void initialize() {
// here temporarily assuming that the force and action solver are the same
this->ActionPtr.reset(new TwoFlavourEvenOddPseudoFermionAction<Impl>(*(this->fop_mod->getPtr()), *(this->solver_mod->getPtr()), *(this->solver_mod->getPtr())));
}
};
@ -308,6 +346,9 @@ static Registrar<QCD::PlaqPlusRectangleGMod, HMC_LGTActionModuleFactory<gauge_st
// FIXME more general implementation
static Registrar<QCD::TwoFlavourFModule<QCD::WilsonImplR> , HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __TwoFlavourFmodXMLInit("TwoFlavours");
static Registrar<QCD::TwoFlavourEOFModule<QCD::WilsonImplR> , HMC_LGTActionModuleFactory<gauge_string, XmlReader> > __TwoFlavourEOFmodXMLInit("TwoFlavoursEvenOdd");
// add here the registration for other implementations and readers

5
lib/qcd/modules/FermionOperatorModules.h
View File

@ -124,8 +124,9 @@ class WilsonFermionModule: public FermionOperatorModule<WilsonFermion, FermionIm
// acquire resource
virtual void initialize(){
typename FermionImpl::GaugeField U(this->GridRefs[0].get().get_full());
this->FOPtr.reset(new WilsonFermion<FermionImpl>(U, *(this->GridRefs[0].get().get_full()), *(this->GridRefs[0].get().get_rb()), this->Par_.mass));
auto &GridMod = this->GridRefs[0].get();
typename FermionImpl::GaugeField U(GridMod.get_full());
this->FOPtr.reset(new WilsonFermion<FermionImpl>(U, *(GridMod.get_full()), *(GridMod.get_rb()), this->Par_.mass));
}
};

13
lib/qcd/modules/Modules.h
View File

@ -43,16 +43,6 @@ class NoParameters{};
/*
Base class for modules with parameters
*/
class ObjectInfo: Serializable {
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ObjectInfo,
std::string, name);
};
template < class P >
class Parametrized{
public:
@ -101,9 +91,6 @@ class Parametrized<NoParameters>{
/*
Lowest level abstract module class
*/
template <class Prod>
class HMCModuleBase {
public: