portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Adding PV pseudofermion in prep for DWF HMC.

Browse Source 
Not compiled this yet, but cloned in from BFM.
This commit is contained in:
Peter Boyle 2015-08-17 23:14:48 +01:00
parent 353d66def1
commit 24382d77bb
1 changed files with 132 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

132
lib/qcd/action/pseudofermion/TwoFlavourRatio.h Normal file
View File

@ -0,0 +1,132 @@
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H
namespace Grid{
namespace QCD{
///////////////////////////////////////
// Two flavour ratio
///////////////////////////////////////
template<class Impl>
class TwoFlavourRatioPseudoFermionAction : public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
private:
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
FermionField Phi; // the pseudo fermion field for this trajectory
public:
TwoFlavourRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS
) : NumOp(_NumOp), DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {};
virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi}
//
// phi = Vdag^{-1} Mdag eta
//
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2).
// and must multiply by 0.707....
//
RealD scale = std::sqrt(0.5);
FermionField eta(NumOp.FermionGrid());
gaussian(pRNG,eta);
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp);
DenOp.Mdag(eta,Phi); // Mdag eta
ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta
NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta
Phi=Phi*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag V (Mdag M)^-1 Vdag phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
X=zero;
NumOp.Mdag(Phi,Y); // Vdag phi
ActionSolver(MdagMOp,Y,X); // MdagMinv Vdag phi
MdagMOp.Op(X,Y); // Y=Mdaginv Vdag phi
RealD action = norm2(Y);
return action;
};
//////////////////////////////////////////////////////
// dS/du = phi^dag dV (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi
// + phi^dag V (Mdag M)^-1 dV^dag phi
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
FermionField f1(NumOp.FermionGrid());
GaugeField force(FermOp.GaugeGrid());
X=zero;
//f1=Vdag phi
NumOp.Mdag(phi,f1);
//X = (Mdag M)^-1 V^dag phi
DerivativeSolver(MdagMOp,f1,X);
//Y = (Mdag)^-1 V^dag phi
DenOp.M(X,Y);
// phi^dag V (Mdag M)^-1 dV^dag phi
NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force;
// phi^dag dV (Mdag M)^-1 V^dag phi
NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force;
// - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi
DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force;
DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force;
dSdU = Ta(dSdU);
};
};
}
}
#endif