portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Namespace, indentation

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:44:06 +00:00
parent be612b3931
commit 6b7e82f1a9
1 changed files with 122 additions and 123 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -26,164 +26,163 @@ with this program; if not, write to the Free Software Foundation, Inc.,
See the full license in the file "LICENSE" in the top level distribution See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
namespace Grid { NAMESPACE_BEGIN(Grid);
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);
private:
FermionOperator<Impl> &FermOp; // the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
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())
{};
virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";}
////////////////////////////////////////////////////////////////////////
// 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
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
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())
{};
virtual std::string action_name(){return "TwoFlavourEvenOddPseudoFermionAction";}
virtual std::string LogParameters(){ virtual std::string LogParameters(){
std::stringstream sstream; std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl; sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
return sstream.str(); return sstream.str();
} }
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied // Push the gauge field in to the dops. Assume any BC's and smearing already applied
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi} // P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
// Phi = McpDag eta // Phi = McpDag eta
// P(eta) = e^{- eta^dag eta} // P(eta) = e^{- eta^dag eta}
// //
// e^{x^2/2 sig^2} => sig^2 = 0.5. // e^{x^2/2 sig^2} => sig^2 = 0.5.
RealD scale = std::sqrt(0.5); RealD scale = std::sqrt(0.5);
FermionField eta (FermOp.FermionGrid()); FermionField eta (FermOp.FermionGrid());
FermionField etaOdd (FermOp.FermionRedBlackGrid()); FermionField etaOdd (FermOp.FermionRedBlackGrid());
FermionField etaEven(FermOp.FermionRedBlackGrid()); FermionField etaEven(FermOp.FermionRedBlackGrid());
gaussian(pRNG,eta); gaussian(pRNG,eta);
pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta); pickCheckerboard(Odd,etaOdd,eta);
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> PCop(FermOp); SchurDifferentiableOperator<Impl> PCop(FermOp);
PCop.MpcDag(etaOdd,PhiOdd); PCop.MpcDag(etaOdd,PhiOdd);
FermOp.MooeeDag(etaEven,PhiEven); FermOp.MooeeDag(etaEven,PhiEven);
PhiOdd =PhiOdd*scale; PhiOdd =PhiOdd*scale;
PhiEven=PhiEven*scale; PhiEven=PhiEven*scale;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
// S = phi^dag (Mdag M)^-1 phi (odd) // S = phi^dag (Mdag M)^-1 phi (odd)
// + phi^dag (Mdag M)^-1 phi (even) // + phi^dag (Mdag M)^-1 phi (even)
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid()); FermionField X(FermOp.FermionRedBlackGrid());
FermionField Y(FermOp.FermionRedBlackGrid()); FermionField Y(FermOp.FermionRedBlackGrid());
SchurDifferentiableOperator<Impl> PCop(FermOp); SchurDifferentiableOperator<Impl> PCop(FermOp);
X=zero; X=zero;
ActionSolver(PCop,PhiOdd,X); ActionSolver(PCop,PhiOdd,X);
PCop.Op(X,Y); PCop.Op(X,Y);
RealD action = norm2(Y); RealD action = norm2(Y);
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
// Only really clover term that creates this. // Only really clover term that creates this.
FermOp.MooeeInvDag(PhiEven,Y); FermOp.MooeeInvDag(PhiEven,Y);
action = action + norm2(Y); action = action + norm2(Y);
std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl; std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl;
return action; return action;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
// //
// dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi // dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi
// = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi // = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi
// //
// = - Ydag dM X - Xdag dMdag Y // = - Ydag dM X - Xdag dMdag Y
// //
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) { virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid()); FermionField X(FermOp.FermionRedBlackGrid());
FermionField Y(FermOp.FermionRedBlackGrid()); FermionField Y(FermOp.FermionRedBlackGrid());
GaugeField tmp(FermOp.GaugeGrid()); GaugeField tmp(FermOp.GaugeGrid());
SchurDifferentiableOperator<Impl> Mpc(FermOp); SchurDifferentiableOperator<Impl> Mpc(FermOp);
// Our conventions really make this UdSdU; We do not differentiate wrt Udag here. // Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
// So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
X=zero; X=zero;
DerivativeSolver(Mpc,PhiOdd,X); DerivativeSolver(Mpc,PhiOdd,X);
Mpc.Mpc(X,Y); Mpc.Mpc(X,Y);
Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp;
Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp;
// Treat the EE case. (MdagM)^-1 = Minv Minvdag // Treat the EE case. (MdagM)^-1 = Minv Minvdag
// Deriv defaults to zero. // Deriv defaults to zero.
// FermOp.MooeeInvDag(PhiOdd,Y); // FermOp.MooeeInvDag(PhiOdd,Y);
// FermOp.MooeeInv(Y,X); // FermOp.MooeeInv(Y,X);
// FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; // FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
// FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; // FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
assert(FermOp.ConstEE() == 1); assert(FermOp.ConstEE() == 1);
/* /*
FermOp.MooeeInvDag(PhiOdd,Y); FermOp.MooeeInvDag(PhiOdd,Y);
FermOp.MooeeInv(Y,X); FermOp.MooeeInv(Y,X);
FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp; FermOp.MeeDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp; FermOp.MeeDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
*/ */
//dSdU = Ta(dSdU); //dSdU = Ta(dSdU);
}; };
}; };
} NAMESPACE_END(Grid);
}
#endif #endif