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

Added a check that the initial EOFA action agrees with |eta|^2, thus checking the quality of the rational approximation in the heatbath

Browse Source
This commit is contained in:
Christopher Kelly 2021-05-18 13:57:44 -04:00
parent 9f0271039f
commit 86f08c6b9a
3 changed files with 31 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Grid/qcd/action/ActionParams.h
View File

@ -66,7 +66,8 @@ struct StaggeredImplParams {
RealD, tolerance,
int, degree,
int, precision,
int, BoundsCheckFreq);
int, BoundsCheckFreq,
RealD, BoundsCheckTol);
// MaxIter and tolerance, vectors??
@ -77,14 +78,16 @@ struct StaggeredImplParams {
RealD tol = 1.0e-8,
int _degree = 10,
int _precision = 64,
int _BoundsCheckFreq=20)
int _BoundsCheckFreq=20,
double _BoundsCheckTol=1e-6)
: lo(_lo),
hi(_hi),
MaxIter(_maxit),
tolerance(tol),
degree(_degree),
precision(_precision),
BoundsCheckFreq(_BoundsCheckFreq){};
BoundsCheckFreq(_BoundsCheckFreq),
BoundsCheckTol(_BoundsCheckTol){};
};

26
Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
View File

@ -64,6 +64,8 @@ NAMESPACE_BEGIN(Grid);
SchurRedBlackDiagMooeeSolve<FermionField> DerivativeSolverR;
FermionField Phi; // the pseudofermion field for this trajectory
RealD norm2_eta; //|eta|^2 where eta is the random gaussian field used to generate the pseudofermion field
bool initial_action; //true for the first call to S after refresh, for which the identity S = |eta|^2 holds provided the rational approx is good
public:
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
@ -87,7 +89,8 @@ NAMESPACE_BEGIN(Grid);
DerivativeSolverL(DerivCGL, false, true), DerivativeSolverR(DerivCGR, false, true),
Phi(_Lop.FermionGrid()),
param(p),
use_heatbath_forecasting(use_fc)
use_heatbath_forecasting(use_fc),
initial_action(false)
{
AlgRemez remez(param.lo, param.hi, param.precision);
@ -216,8 +219,14 @@ NAMESPACE_BEGIN(Grid);
Lop.RefreshShiftCoefficients(0.0);
Rop.RefreshShiftCoefficients(-1.0);
//Mark that the next call to S is the first after refresh
initial_action = true;
// Bounds check
RealD EtaDagEta = norm2(eta);
norm2_eta = EtaDagEta;
// RealD PhiDagMPhi= norm2(eta);
};
@ -290,6 +299,21 @@ NAMESPACE_BEGIN(Grid);
Rop.Omega(tmp[1], tmp[0], 1, 1);
action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real();
if(initial_action){
//For the first call to S after refresh, S = |eta|^2. We can use this to ensure the rational approx is good
RealD diff = action - norm2_eta;
//S_init = eta^dag M^{-1/2} M M^{-1/2} eta
//S_init - eta^dag eta = eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta
RealD test = sqrt(fabs(diff)/norm2_eta); //test the quality of the rational approx
std::cout << GridLogMessage << action_name() << " initial action " << action << " expect " << norm2_eta << "; diff " << diff << std::endl;
std::cout << GridLogMessage << action_name() << " sqrt( eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta )/sqrt( eta^dag eta ) = " << test << " expect 0 (tol " << param.BoundsCheckTol << ")" << std::endl;
assert( ( test < param.BoundsCheckTol ) && " Initial action check failed" );
initial_action = false;
}
return action;
};

15
tests/forces/Test_mobius_force_eofa.cc
View File

@ -132,21 +132,6 @@ int main (int argc, char** argv)
Meofa.refresh(U, sRNG, RNG5 );
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"