mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-23 18:22:02 +01:00
fd933420c6
Added a bounds-check function for the RHMC with arbitrary power
Added a pseudofermion action for the rational ratio with an arbitrary power and a mixed-precision variant of the same. The existing one-flavor rational ratio class now uses the general class under the hood
To support testing of the two-flavor even-odd ratio pseudofermion, separated the functionality of generating the random field and performing the heatbath step, and added a method to obtain the pseudofermion field
Added a new HMC runner start type: CheckpointStartReseed, which reseeds the RNG from scratch, allowing for the creation of new evolution streams from an existing checkpoint. Added log output of seeds used when the RNG is seeded.
EOFA changes:
To support mixed-precision inversion, generalized the class to maintain a separate solver for the L and R operators in the heatbath (separate solvers are already implemented for the other stages)
To support mixed-precision, the action of setting the operator shift coefficients is now maintained in a virtual function. A derived class for mixed-precision solvers ensures the coefficients are applied to both the double and single-prec operators
The ||^2 of the random source is now stored by the heatbath and compared to the initial action when it is computed. These should be equal but may differ if the rational bounds are not chosen correctly, hence serving as a useful and free test
Fixed calculation of M_eofa (previously incomplete and #if'd out)
Added functionality to compute M_eofa^-1 to complement the calculation of M_eofa (both are equally expensive!)
To support testing, separated the functionality of generating the random field and performing the heatbath step, and added a method to obtain the pseudofermion field
Added a test program which computes the G-parity force using the 1 and 2 flavor implementations and compares the result. Test supports DWF, EOFA and DSDR actions, chosen by a command line option.
The Mobius EOFA force test now also checks the rational approximation used for the heatbath
Added a test program for the mixed precision EOFA compared to the double-prec implementation,
G-parity HMC test now applied GPBC in the y direction and not the t direction (GPBC in t are no longer supported) and checkpoints after every configuration
Added a test program which computes the two-flavor G-parity action (via RHMC) with both the 1 and 2 flavor implementations and checks they agree
Added a test program to check the implementation of M_eofa^{-1}
220 lines
8.0 KiB
C++
220 lines
8.0 KiB
C++
/*************************************************************************************
|
|
|
|
Grid physics library, www.github.com/paboyle/Grid
|
|
|
|
Source file: ./lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
|
|
|
|
Copyright (C) 2015
|
|
|
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
|
Author: paboyle <paboyle@ph.ed.ac.uk>
|
|
|
|
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 QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H
|
|
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_RATIO_H
|
|
|
|
NAMESPACE_BEGIN(Grid);
|
|
|
|
///////////////////////////////////////
|
|
// Two flavour ratio
|
|
///////////////////////////////////////
|
|
template<class Impl>
|
|
class TwoFlavourEvenOddRatioPseudoFermionAction : 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;
|
|
OperatorFunction<FermionField> &HeatbathSolver;
|
|
|
|
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
|
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
|
|
|
public:
|
|
TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp,
|
|
FermionOperator<Impl> &_DenOp,
|
|
OperatorFunction<FermionField> & DS,
|
|
OperatorFunction<FermionField> & AS ) :
|
|
TwoFlavourEvenOddRatioPseudoFermionAction(_NumOp,_DenOp, DS,AS,AS) {};
|
|
|
|
TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp,
|
|
FermionOperator<Impl> &_DenOp,
|
|
OperatorFunction<FermionField> & DS,
|
|
OperatorFunction<FermionField> & AS, OperatorFunction<FermionField> & HS) :
|
|
NumOp(_NumOp),
|
|
DenOp(_DenOp),
|
|
DerivativeSolver(DS),
|
|
ActionSolver(AS),
|
|
HeatbathSolver(HS),
|
|
PhiEven(_NumOp.FermionRedBlackGrid()),
|
|
PhiOdd(_NumOp.FermionRedBlackGrid())
|
|
{
|
|
conformable(_NumOp.FermionGrid(), _DenOp.FermionGrid());
|
|
conformable(_NumOp.FermionRedBlackGrid(), _DenOp.FermionRedBlackGrid());
|
|
conformable(_NumOp.GaugeGrid(), _DenOp.GaugeGrid());
|
|
conformable(_NumOp.GaugeRedBlackGrid(), _DenOp.GaugeRedBlackGrid());
|
|
};
|
|
|
|
virtual std::string action_name(){
|
|
std::stringstream sstream;
|
|
sstream<<"TwoFlavourEvenOddRatioPseudoFermionAction det("<<DenOp.Mass()<<") / det("<<NumOp.Mass()<<")";
|
|
return sstream.str();
|
|
}
|
|
|
|
virtual std::string LogParameters(){
|
|
std::stringstream sstream;
|
|
sstream<< GridLogMessage << "["<<action_name()<<"] -- No further parameters "<<std::endl;
|
|
return sstream.str();
|
|
}
|
|
|
|
|
|
const FermionField &getPhiOdd() const{ return PhiOdd; }
|
|
|
|
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
|
|
// P(eta_o) = e^{- eta_o^dag eta_o}
|
|
//
|
|
// e^{x^2/2 sig^2} => sig^2 = 0.5.
|
|
//
|
|
RealD scale = std::sqrt(0.5);
|
|
|
|
FermionField eta (NumOp.FermionGrid());
|
|
gaussian(pRNG,eta); eta = eta * scale;
|
|
|
|
refresh(U,eta);
|
|
}
|
|
|
|
void refresh(const GaugeField &U, const FermionField &eta) {
|
|
|
|
// P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi}
|
|
//
|
|
// NumOp == V
|
|
// DenOp == M
|
|
//
|
|
FermionField etaOdd (NumOp.FermionRedBlackGrid());
|
|
FermionField etaEven(NumOp.FermionRedBlackGrid());
|
|
FermionField tmp (NumOp.FermionRedBlackGrid());
|
|
|
|
pickCheckerboard(Even,etaEven,eta);
|
|
pickCheckerboard(Odd,etaOdd,eta);
|
|
|
|
NumOp.ImportGauge(U);
|
|
DenOp.ImportGauge(U);
|
|
|
|
SchurDifferentiableOperator<Impl> Mpc(DenOp);
|
|
SchurDifferentiableOperator<Impl> Vpc(NumOp);
|
|
|
|
// Odd det factors
|
|
Mpc.MpcDag(etaOdd,PhiOdd);
|
|
tmp=Zero();
|
|
HeatbathSolver(Vpc,PhiOdd,tmp);
|
|
Vpc.Mpc(tmp,PhiOdd);
|
|
|
|
// Even det factors
|
|
DenOp.MooeeDag(etaEven,tmp);
|
|
NumOp.MooeeInvDag(tmp,PhiEven);
|
|
};
|
|
|
|
//////////////////////////////////////////////////////
|
|
// S = phi^dag V (Mdag M)^-1 Vdag phi
|
|
//////////////////////////////////////////////////////
|
|
virtual RealD S(const GaugeField &U) {
|
|
|
|
NumOp.ImportGauge(U);
|
|
DenOp.ImportGauge(U);
|
|
|
|
SchurDifferentiableOperator<Impl> Mpc(DenOp);
|
|
SchurDifferentiableOperator<Impl> Vpc(NumOp);
|
|
|
|
FermionField X(NumOp.FermionRedBlackGrid());
|
|
FermionField Y(NumOp.FermionRedBlackGrid());
|
|
|
|
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
|
|
X=Zero();
|
|
ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
|
|
//Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
|
|
// Multiply by Ydag
|
|
RealD action = real(innerProduct(Y,X));
|
|
|
|
//RealD action = norm2(Y);
|
|
|
|
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
|
|
// Only really clover term that creates this. Leave the EE portion as a future to do to make most
|
|
// rapid progresss on DWF for now.
|
|
//
|
|
NumOp.MooeeDag(PhiEven,X);
|
|
DenOp.MooeeInvDag(X,Y);
|
|
action = 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);
|
|
|
|
SchurDifferentiableOperator<Impl> Mpc(DenOp);
|
|
SchurDifferentiableOperator<Impl> Vpc(NumOp);
|
|
|
|
FermionField X(NumOp.FermionRedBlackGrid());
|
|
FermionField Y(NumOp.FermionRedBlackGrid());
|
|
|
|
// This assignment is necessary to be compliant with the HMC grids
|
|
GaugeField force(dSdU.Grid());
|
|
|
|
//Y=Vdag phi
|
|
//X = (Mdag M)^-1 V^dag phi
|
|
//Y = (Mdag)^-1 V^dag phi
|
|
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
|
|
X=Zero();
|
|
DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
|
|
Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
|
|
|
|
// phi^dag V (Mdag M)^-1 dV^dag phi
|
|
Vpc.MpcDagDeriv(force , X, PhiOdd ); dSdU = force;
|
|
|
|
// phi^dag dV (Mdag M)^-1 V^dag phi
|
|
Vpc.MpcDeriv(force , PhiOdd, X ); 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
|
|
Mpc.MpcDeriv(force,Y,X); dSdU = dSdU-force;
|
|
Mpc.MpcDagDeriv(force,X,Y); dSdU = dSdU-force;
|
|
|
|
// FIXME No force contribution from EvenEven assumed here
|
|
// Needs a fix for clover.
|
|
assert(NumOp.ConstEE() == 1);
|
|
assert(DenOp.ConstEE() == 1);
|
|
|
|
dSdU = -dSdU;
|
|
|
|
};
|
|
};
|
|
NAMESPACE_END(Grid);
|
|
#endif
|