portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 14:04:32 +00:00
Code Releases Activity

Merge branch 'feature/dirichlet-gparity' of https://github.com/paboyle/Grid into feature/dirichlet-gparity

Browse Source
This commit is contained in:
Peter Boyle
2022-07-01 09:55:43 -04:00
parent 3d8146b596 7319d4e1ad
commit 4df1e0987f
17 changed files with 1920 additions and 310 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -37,6 +37,7 @@ NAMESPACE_BEGIN(Grid);
// These can move into a params header and be given MacroMagic serialisation
struct GparityWilsonImplParams {
Coordinate twists;
//mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
GparityWilsonImplParams() : twists(Nd, 0) {};
};
@@ -66,7 +67,8 @@ struct StaggeredImplParams {
RealD, mdtolerance,
int, degree,
int, precision,
int, BoundsCheckFreq);
int, BoundsCheckFreq,
RealD, BoundsCheckTol);
// MaxIter and tolerance, vectors??
@@ -78,7 +80,8 @@ struct StaggeredImplParams {
int _degree = 10,
int _precision = 64,
int _BoundsCheckFreq=20,
RealD mdtol = 1.0e-6)
RealD mdtol = 1.0e-6,
double _BoundsCheckTol=1e-6)
: lo(_lo),
hi(_hi),
MaxIter(_maxit),
@@ -86,9 +89,52 @@ struct StaggeredImplParams {
mdtolerance(mdtol),
degree(_degree),
precision(_precision),
BoundsCheckFreq(_BoundsCheckFreq),
BoundsCheckTol(_BoundsCheckTol){};
};
/*Action parameters for the generalized rational action
The approximation is for (M^dag M)^{1/inv_pow}
where inv_pow is the denominator of the fractional power.
Default inv_pow=2 for square root, making this equivalent to
the OneFlavourRational action
*/
struct RationalActionParams : Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(RationalActionParams,
int, inv_pow,
RealD, lo, //low eigenvalue bound of rational approx
RealD, hi, //high eigenvalue bound of rational approx
int, MaxIter, //maximum iterations in msCG
RealD, action_tolerance, //msCG tolerance in action evaluation
int, action_degree, //rational approx tolerance in action evaluation
RealD, md_tolerance, //msCG tolerance in MD integration
int, md_degree, //rational approx tolerance in MD integration
int, precision, //precision of floating point arithmetic
int, BoundsCheckFreq); //frequency the approximation is tested (with Metropolis degree/tolerance); 0 disables the check
// constructor
RationalActionParams(int _inv_pow = 2,
RealD _lo = 0.0,
RealD _hi = 1.0,
int _maxit = 1000,
RealD _action_tolerance = 1.0e-8,
int _action_degree = 10,
RealD _md_tolerance = 1.0e-8,
int _md_degree = 10,
int _precision = 64,
int _BoundsCheckFreq=20)
: inv_pow(_inv_pow),
lo(_lo),
hi(_hi),
MaxIter(_maxit),
action_tolerance(_action_tolerance),
action_degree(_action_degree),
md_tolerance(_md_tolerance),
md_degree(_md_degree),
precision(_precision),
BoundsCheckFreq(_BoundsCheckFreq){};
};
NAMESPACE_END(Grid);
#endif

60
Grid/qcd/action/pseudofermion/Bounds.h
View File

@@ -65,13 +65,65 @@ NAMESPACE_BEGIN(Grid);
X=X-Y;
RealD Nd = norm2(X);
std::cout << "************************* "<<std::endl;
std::cout << " noise = "<<Nx<<std::endl;
std::cout << " (MdagM^-1/2)^2 noise = "<<Nz<<std::endl;
std::cout << " MdagM (MdagM^-1/2)^2 noise = "<<Ny<<std::endl;
std::cout << " noise - MdagM (MdagM^-1/2)^2 noise = "<<Nd<<std::endl;
std::cout << " | noise |^2 = "<<Nx<<std::endl;
std::cout << " | (MdagM^-1/2)^2 noise |^2 = "<<Nz<<std::endl;
std::cout << " | MdagM (MdagM^-1/2)^2 noise |^2 = "<<Ny<<std::endl;
std::cout << " | noise - MdagM (MdagM^-1/2)^2 noise |^2 = "<<Nd<<std::endl;
std::cout << " | noise - MdagM (MdagM^-1/2)^2 noise|/|noise| = " << std::sqrt(Nd/Nx) << std::endl;
std::cout << "************************* "<<std::endl;
assert( (std::sqrt(Nd/Nx)<tol) && " InverseSqrtBoundsCheck ");
}
/* For a HermOp = M^dag M, check the approximation of HermOp^{-1/inv_pow}
by computing |X - HermOp * [ Hermop^{-1/inv_pow} ]^{inv_pow} X| < tol
for noise X (aka GaussNoise).
ApproxNegPow should be the rational approximation for X^{-1/inv_pow}
*/
template<class Field> void InversePowerBoundsCheck(int inv_pow,
int MaxIter,double tol,
LinearOperatorBase<Field> &HermOp,
Field &GaussNoise,
MultiShiftFunction &ApproxNegPow)
{
GridBase *FermionGrid = GaussNoise.Grid();
Field X(FermionGrid);
Field Y(FermionGrid);
Field Z(FermionGrid);
Field tmp1(FermionGrid), tmp2(FermionGrid);
X=GaussNoise;
RealD Nx = norm2(X);
ConjugateGradientMultiShift<Field> msCG(MaxIter,ApproxNegPow);
tmp1 = X;
Field* in = &tmp1;
Field* out = &tmp2;
for(int i=0;i<inv_pow;i++){ //apply [ Hermop^{-1/inv_pow} ]^{inv_pow} X = HermOp^{-1} X
msCG(HermOp, *in, *out); //backwards conventions!
if(i!=inv_pow-1) std::swap(in, out);
}
Z = *out;
RealD Nz = norm2(Z);
HermOp.HermOp(Z,Y);
RealD Ny = norm2(Y);
X=X-Y;
RealD Nd = norm2(X);
std::cout << "************************* "<<std::endl;
std::cout << " | noise |^2 = "<<Nx<<std::endl;
std::cout << " | (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2 = "<<Nz<<std::endl;
std::cout << " | MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2 = "<<Ny<<std::endl;
std::cout << " | noise - MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |^2 = "<<Nd<<std::endl;
std::cout << " | noise - MdagM (MdagM^-1/" << inv_pow << ")^" << inv_pow << " noise |/| noise | = "<<std::sqrt(Nd/Nx)<<std::endl;
std::cout << "************************* "<<std::endl;
assert( (std::sqrt(Nd/Nx)<tol) && " InversePowerBoundsCheck ");
}
NAMESPACE_END(Grid);

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

@@ -44,6 +44,10 @@ NAMESPACE_BEGIN(Grid);
// Exact one flavour implementation of DWF determinant ratio //
///////////////////////////////////////////////////////////////
//Note: using mixed prec CG for the heatbath solver in this action class will not work
// because the L, R operators must have their shift coefficients updated throughout the heatbath step
// You will find that the heatbath solver simply won't converge.
// To use mixed precision here use the ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction variant below
template<class Impl>
class ExactOneFlavourRatioPseudoFermionAction : public Action<typename Impl::GaugeField>
{
@@ -57,37 +61,60 @@ NAMESPACE_BEGIN(Grid);
bool use_heatbath_forecasting;
AbstractEOFAFermion<Impl>& Lop; // the basic LH operator
AbstractEOFAFermion<Impl>& Rop; // the basic RH operator
SchurRedBlackDiagMooeeSolve<FermionField> SolverHB;
SchurRedBlackDiagMooeeSolve<FermionField> SolverHBL;
SchurRedBlackDiagMooeeSolve<FermionField> SolverHBR;
SchurRedBlackDiagMooeeSolve<FermionField> SolverL;
SchurRedBlackDiagMooeeSolve<FermionField> SolverR;
SchurRedBlackDiagMooeeSolve<FermionField> DerivativeSolverL;
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:
//Used in the heatbath, refresh the shift coefficients of the L (LorR=0) or R (LorR=1) operator
virtual void heatbathRefreshShiftCoefficients(int LorR, RealD to){
AbstractEOFAFermion<Impl>&op = LorR == 0 ? Lop : Rop;
op.RefreshShiftCoefficients(to);
}
//Use the same solver for L,R in all cases
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
AbstractEOFAFermion<Impl>& _Rop,
OperatorFunction<FermionField>& CG,
Params& p,
bool use_fc=false)
: ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,CG,CG,CG,CG,CG,p,use_fc) {};
: ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,CG,CG,CG,CG,CG,CG,p,use_fc) {};
//Use the same solver for L,R in the heatbath but different solvers elsewhere
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
AbstractEOFAFermion<Impl>& _Rop,
OperatorFunction<FermionField>& HeatbathCG,
OperatorFunction<FermionField>& ActionCGL, OperatorFunction<FermionField>& ActionCGR,
OperatorFunction<FermionField>& DerivCGL , OperatorFunction<FermionField>& DerivCGR,
Params& p,
bool use_fc=false)
: ExactOneFlavourRatioPseudoFermionAction(_Lop,_Rop,HeatbathCG,HeatbathCG, ActionCGL, ActionCGR, DerivCGL,DerivCGR,p,use_fc) {};
//Use different solvers for L,R in all cases
ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop,
AbstractEOFAFermion<Impl>& _Rop,
OperatorFunction<FermionField>& HeatbathCGL, OperatorFunction<FermionField>& HeatbathCGR,
OperatorFunction<FermionField>& ActionCGL, OperatorFunction<FermionField>& ActionCGR,
OperatorFunction<FermionField>& DerivCGL , OperatorFunction<FermionField>& DerivCGR,
Params& p,
bool use_fc=false) :
Lop(_Lop),
Rop(_Rop),
SolverHB(HeatbathCG,false,true),
SolverHBL(HeatbathCGL,false,true), SolverHBR(HeatbathCGR,false,true),
SolverL(ActionCGL, false, true), SolverR(ActionCGR, false, true),
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);
@@ -97,6 +124,8 @@ NAMESPACE_BEGIN(Grid);
PowerNegHalf.Init(remez, param.tolerance, true);
};
const FermionField &getPhi() const{ return Phi; }
virtual std::string action_name() { return "ExactOneFlavourRatioPseudoFermionAction"; }
virtual std::string LogParameters() {
@@ -117,6 +146,19 @@ NAMESPACE_BEGIN(Grid);
else{ for(int s=0; s<Ls; ++s){ axpby_ssp_pminus(out, 0.0, in, 1.0, in, s, s); } }
}
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 (Lop.FermionGrid());
gaussian(pRNG,eta); eta = eta * scale;
refresh(U,eta);
}
// EOFA heatbath: see Eqn. (29) of arXiv:1706.05843
// We generate a Gaussian noise vector \eta, and then compute
// \Phi = M_{\rm EOFA}^{-1/2} * \eta
@@ -124,12 +166,10 @@ NAMESPACE_BEGIN(Grid);
//
// As a check of rational require \Phi^dag M_{EOFA} \Phi == eta^dag M^-1/2^dag M M^-1/2 eta = eta^dag eta
//
virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
{
void refresh(const GaugeField &U, const FermionField &eta) {
Lop.ImportGauge(U);
Rop.ImportGauge(U);
FermionField eta (Lop.FermionGrid());
FermionField CG_src (Lop.FermionGrid());
FermionField CG_soln (Lop.FermionGrid());
FermionField Forecast_src(Lop.FermionGrid());
@@ -140,11 +180,6 @@ NAMESPACE_BEGIN(Grid);
if(use_heatbath_forecasting){ prev_solns.reserve(param.degree); }
ChronoForecast<AbstractEOFAFermion<Impl>, FermionField> Forecast;
// Seed with Gaussian noise vector (var = 0.5)
RealD scale = std::sqrt(0.5);
gaussian(pRNG,eta);
eta = eta * scale;
// \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta
RealD N(PowerNegHalf.norm);
for(int k=0; k<param.degree; ++k){ N += PowerNegHalf.residues[k] / ( 1.0 + PowerNegHalf.poles[k] ); }
@@ -160,15 +195,15 @@ NAMESPACE_BEGIN(Grid);
tmp[1] = Zero();
for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
Lop.RefreshShiftCoefficients(-gamma_l);
heatbathRefreshShiftCoefficients(0, -gamma_l);
if(use_heatbath_forecasting){ // Forecast CG guess using solutions from previous poles
Lop.Mdag(CG_src, Forecast_src);
CG_soln = Forecast(Lop, Forecast_src, prev_solns);
SolverHB(Lop, CG_src, CG_soln);
SolverHBL(Lop, CG_src, CG_soln);
prev_solns.push_back(CG_soln);
} else {
CG_soln = Zero(); // Just use zero as the initial guess
SolverHB(Lop, CG_src, CG_soln);
SolverHBL(Lop, CG_src, CG_soln);
}
Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
tmp[1] = tmp[1] + ( PowerNegHalf.residues[k]*gamma_l*gamma_l*Lop.k ) * tmp[0];
@@ -187,15 +222,15 @@ NAMESPACE_BEGIN(Grid);
if(use_heatbath_forecasting){ prev_solns.clear(); } // empirically, LH solns don't help for RH solves
for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
Rop.RefreshShiftCoefficients(-gamma_l*PowerNegHalf.poles[k]);
heatbathRefreshShiftCoefficients(1, -gamma_l*PowerNegHalf.poles[k]);
if(use_heatbath_forecasting){
Rop.Mdag(CG_src, Forecast_src);
CG_soln = Forecast(Rop, Forecast_src, prev_solns);
SolverHB(Rop, CG_src, CG_soln);
SolverHBR(Rop, CG_src, CG_soln);
prev_solns.push_back(CG_soln);
} else {
CG_soln = Zero();
SolverHB(Rop, CG_src, CG_soln);
SolverHBR(Rop, CG_src, CG_soln);
}
Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
tmp[1] = tmp[1] - ( PowerNegHalf.residues[k]*gamma_l*gamma_l*Rop.k ) * tmp[0];
@@ -205,49 +240,117 @@ NAMESPACE_BEGIN(Grid);
Phi = Phi + tmp[1];
// Reset shift coefficients for energy and force evals
Lop.RefreshShiftCoefficients(0.0);
Rop.RefreshShiftCoefficients(-1.0);
heatbathRefreshShiftCoefficients(0, 0.0);
heatbathRefreshShiftCoefficients(1, -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);
};
void Meofa(const GaugeField& U,const FermionField &phi, FermionField & Mphi)
void Meofa(const GaugeField& U,const FermionField &in, FermionField & out)
{
#if 0
Lop.ImportGauge(U);
Rop.ImportGauge(U);
FermionField spProj_Phi(Lop.FermionGrid());
FermionField mPhi(Lop.FermionGrid());
FermionField spProj_in(Lop.FermionGrid());
std::vector<FermionField> tmp(2, Lop.FermionGrid());
mPhi = phi;
out = in;
// LH term: S = S - k <\Phi| P_{-} \Omega_{-}^{\dagger} H(mf)^{-1} \Omega_{-} P_{-} |\Phi>
spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, tmp[0], -1, 0);
spProj(in, spProj_in, -1, Lop.Ls);
Lop.Omega(spProj_in, tmp[0], -1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = Zero();
SolverL(Lop, tmp[1], tmp[0]);
Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back
Lop.Omega(tmp[1], tmp[0], -1, 1);
mPhi = mPhi - Lop.k * innerProduct(spProj_Phi, tmp[0]).real();
spProj(tmp[0], tmp[1], -1, Lop.Ls);
out = out - Lop.k * tmp[1];
// RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb)
// - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi>
spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, tmp[0], 1, 0);
// - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi>
spProj(in, spProj_in, 1, Rop.Ls);
Rop.Omega(spProj_in, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]);
tmp[0] = Zero();
SolverR(Rop, tmp[1], tmp[0]);
Rop.Dtilde(tmp[0], tmp[1]);
Rop.Omega(tmp[1], tmp[0], 1, 1);
action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real();
#endif
spProj(tmp[0], tmp[1], 1, Rop.Ls);
out = out + Rop.k * tmp[1];
}
//Due to the structure of EOFA, it is no more expensive to compute the inverse of Meofa
//To ensure correctness we can simply reuse the heatbath code but use the rational approx
//f(x) = 1/x which corresponds to alpha_0=0, alpha_1=1, beta_1=0 => gamma_1=1
void MeofaInv(const GaugeField &U, const FermionField &in, FermionField &out) {
Lop.ImportGauge(U);
Rop.ImportGauge(U);
FermionField CG_src (Lop.FermionGrid());
FermionField CG_soln (Lop.FermionGrid());
std::vector<FermionField> tmp(2, Lop.FermionGrid());
// \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta
// = 1 * \eta
out = in;
// LH terms:
// \Phi = \Phi + k \sum_{k=1}^{N_{p}} P_{-} \Omega_{-}^{\dagger} ( H(mf)
// - \gamma_{l} \Delta_{-}(mf,mb) P_{-} )^{-1} \Omega_{-} P_{-} \eta
spProj(in, tmp[0], -1, Lop.Ls);
Lop.Omega(tmp[0], tmp[1], -1, 0);
G5R5(CG_src, tmp[1]);
{
heatbathRefreshShiftCoefficients(0, -1.); //-gamma_1 = -1.
CG_soln = Zero(); // Just use zero as the initial guess
SolverHBL(Lop, CG_src, CG_soln);
Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
tmp[1] = Lop.k * tmp[0];
}
Lop.Omega(tmp[1], tmp[0], -1, 1);
spProj(tmp[0], tmp[1], -1, Lop.Ls);
out = out + tmp[1];
// RH terms:
// \Phi = \Phi - k \sum_{k=1}^{N_{p}} P_{+} \Omega_{+}^{\dagger} ( H(mb)
// - \beta_l\gamma_{l} \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} \eta
spProj(in, tmp[0], 1, Rop.Ls);
Rop.Omega(tmp[0], tmp[1], 1, 0);
G5R5(CG_src, tmp[1]);
{
heatbathRefreshShiftCoefficients(1, 0.); //-gamma_1 * beta_1 = 0
CG_soln = Zero();
SolverHBR(Rop, CG_src, CG_soln);
Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
tmp[1] = - Rop.k * tmp[0];
}
Rop.Omega(tmp[1], tmp[0], 1, 1);
spProj(tmp[0], tmp[1], 1, Rop.Ls);
out = out + tmp[1];
// Reset shift coefficients for energy and force evals
heatbathRefreshShiftCoefficients(0, 0.0);
heatbathRefreshShiftCoefficients(1, -1.0);
};
// EOFA action: see Eqn. (10) of arXiv:1706.05843
virtual RealD S(const GaugeField& U)
{
@@ -271,7 +374,7 @@ NAMESPACE_BEGIN(Grid);
action -= Lop.k * innerProduct(spProj_Phi, tmp[0]).real();
// RH term: S = S + k <\Phi| P_{+} \Omega_{+}^{\dagger} ( H(mb)
// - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{-} P_{-} |\Phi>
// - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} |\Phi>
spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]);
@@ -281,6 +384,26 @@ 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
//If approximate solution
//S_init - eta^dag eta = eta^dag ( [M^{-1/2}+\delta M^{-1/2}] M [M^{-1/2}+\delta M^{-1/2}] - 1 ) eta
// \approx eta^dag ( \delta M^{-1/2} M^{1/2} + M^{1/2}\delta M^{-1/2} ) eta
// We divide out |eta|^2 to remove source scaling but the tolerance on this check should still be somewhat higher than the actual approx tolerance
RealD test = 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() << "[ eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta ]/|eta^2| = " << test << " expect 0 (tol " << param.BoundsCheckTol << ")" << std::endl;
assert( ( test < param.BoundsCheckTol ) && " Initial action check failed" );
initial_action = false;
}
return action;
};
@@ -329,6 +452,40 @@ NAMESPACE_BEGIN(Grid);
};
};
template<class ImplD, class ImplF>
class ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction : public ExactOneFlavourRatioPseudoFermionAction<ImplD>{
public:
INHERIT_IMPL_TYPES(ImplD);
typedef OneFlavourRationalParams Params;
private:
AbstractEOFAFermion<ImplF>& LopF; // the basic LH operator
AbstractEOFAFermion<ImplF>& RopF; // the basic RH operator
public:
virtual std::string action_name() { return "ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction"; }
//Used in the heatbath, refresh the shift coefficients of the L (LorR=0) or R (LorR=1) operator
virtual void heatbathRefreshShiftCoefficients(int LorR, RealD to){
AbstractEOFAFermion<ImplF> &op = LorR == 0 ? LopF : RopF;
op.RefreshShiftCoefficients(to);
this->ExactOneFlavourRatioPseudoFermionAction<ImplD>::heatbathRefreshShiftCoefficients(LorR,to);
}
ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction(AbstractEOFAFermion<ImplF>& _LopF,
AbstractEOFAFermion<ImplF>& _RopF,
AbstractEOFAFermion<ImplD>& _LopD,
AbstractEOFAFermion<ImplD>& _RopD,
OperatorFunction<FermionField>& HeatbathCGL, OperatorFunction<FermionField>& HeatbathCGR,
OperatorFunction<FermionField>& ActionCGL, OperatorFunction<FermionField>& ActionCGR,
OperatorFunction<FermionField>& DerivCGL , OperatorFunction<FermionField>& DerivCGR,
Params& p,
bool use_fc=false) :
LopF(_LopF), RopF(_RopF), ExactOneFlavourRatioPseudoFermionAction<ImplD>(_LopD, _RopD, HeatbathCGL, HeatbathCGR, ActionCGL, ActionCGR, DerivCGL, DerivCGR, p, use_fc){}
};
NAMESPACE_END(Grid);
#endif

372
Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h Normal file
View File

@@ -0,0 +1,372 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <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_GENERAL_EVEN_ODD_RATIONAL_RATIO_H
#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_H
NAMESPACE_BEGIN(Grid);
/////////////////////////////////////////////////////////
// Generic rational approximation for ratios of operators
/////////////////////////////////////////////////////////
/* S_f = -log( det( [M^dag M]/[V^dag V] )^{1/inv_pow} )
= chi^dag ( [M^dag M]/[V^dag V] )^{-1/inv_pow} chi\
= chi^dag ( [V^dag V]^{-1/2} [M^dag M] [V^dag V]^{-1/2} )^{-1/inv_pow} chi\
= chi^dag [V^dag V]^{1/(2*inv_pow)} [M^dag M]^{-1/inv_pow} [V^dag V]^{1/(2*inv_pow)} chi\
S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
BIG WARNING:
Here V^dag V is referred to in this code as the "numerator" operator and M^dag M is the *denominator* operator.
this refers to their position in the pseudofermion action, which is the *inverse* of what appears in the determinant
Thus for DWF the numerator operator is the Pauli-Villars operator
Here P/Q \sim R_{1/(2*inv_pow)} ~ (V^dagV)^{1/(2*inv_pow)}
Here N/D \sim R_{-1/inv_pow} ~ (M^dagM)^{-1/inv_pow}
*/
template<class Impl>
class GeneralEvenOddRatioRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
typedef RationalActionParams Params;
Params param;
//For action evaluation
MultiShiftFunction ApproxPowerAction ; //rational approx for X^{1/inv_pow}
MultiShiftFunction ApproxNegPowerAction; //rational approx for X^{-1/inv_pow}
MultiShiftFunction ApproxHalfPowerAction; //rational approx for X^{1/(2*inv_pow)}
MultiShiftFunction ApproxNegHalfPowerAction; //rational approx for X^{-1/(2*inv_pow)}
//For the MD integration
MultiShiftFunction ApproxPowerMD ; //rational approx for X^{1/inv_pow}
MultiShiftFunction ApproxNegPowerMD; //rational approx for X^{-1/inv_pow}
MultiShiftFunction ApproxHalfPowerMD; //rational approx for X^{1/(2*inv_pow)}
MultiShiftFunction ApproxNegHalfPowerMD; //rational approx for X^{-1/(2*inv_pow)}
private:
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
FermionField PhiEven; // the pseudo fermion field for this trajectory
FermionField PhiOdd; // the pseudo fermion field for this trajectory
//Generate the approximation to x^{1/inv_pow} (->approx) and x^{-1/inv_pow} (-> approx_inv) by an approx_degree degree rational approximation
//CG_tolerance is used to issue a warning if the approximation error is larger than the tolerance of the CG and is otherwise just stored in the MultiShiftFunction for use by the multi-shift
static void generateApprox(MultiShiftFunction &approx, MultiShiftFunction &approx_inv, int inv_pow, int approx_degree, double CG_tolerance, AlgRemez &remez){
std::cout<<GridLogMessage << "Generating degree "<< approx_degree<<" approximation for x^(1/" << inv_pow << ")"<<std::endl;
double error = remez.generateApprox(approx_degree,1,inv_pow);
if(error > CG_tolerance)
std::cout<<GridLogMessage << "WARNING: Remez approximation has a larger error " << error << " than the CG tolerance " << CG_tolerance << "! Try increasing the number of poles" << std::endl;
approx.Init(remez, CG_tolerance,false);
approx_inv.Init(remez, CG_tolerance,true);
}
protected:
static constexpr bool Numerator = true;
static constexpr bool Denominator = false;
//Allow derived classes to override the multishift CG
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionField &in, FermionField &out){
SchurDifferentiableOperator<Impl> schurOp(numerator ? NumOp : DenOp);
ConjugateGradientMultiShift<FermionField> msCG(MaxIter, approx);
msCG(schurOp,in, out);
}
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionField &in, std::vector<FermionField> &out_elems, FermionField &out){
SchurDifferentiableOperator<Impl> schurOp(numerator ? NumOp : DenOp);
ConjugateGradientMultiShift<FermionField> msCG(MaxIter, approx);
msCG(schurOp,in, out_elems, out);
}
//Allow derived classes to override the gauge import
virtual void ImportGauge(const GaugeField &U){
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
}
public:
GeneralEvenOddRatioRationalPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
const Params & p
) :
NumOp(_NumOp),
DenOp(_DenOp),
PhiOdd (_NumOp.FermionRedBlackGrid()),
PhiEven(_NumOp.FermionRedBlackGrid()),
param(p)
{
std::cout<<GridLogMessage << action_name() << " initialize: starting" << std::endl;
AlgRemez remez(param.lo,param.hi,param.precision);
//Generate approximations for action eval
generateApprox(ApproxPowerAction, ApproxNegPowerAction, param.inv_pow, param.action_degree, param.action_tolerance, remez);
generateApprox(ApproxHalfPowerAction, ApproxNegHalfPowerAction, 2*param.inv_pow, param.action_degree, param.action_tolerance, remez);
//Generate approximations for MD
if(param.md_degree != param.action_degree){ //note the CG tolerance is unrelated to the stopping condition of the Remez algorithm
generateApprox(ApproxPowerMD, ApproxNegPowerMD, param.inv_pow, param.md_degree, param.md_tolerance, remez);
generateApprox(ApproxHalfPowerMD, ApproxNegHalfPowerMD, 2*param.inv_pow, param.md_degree, param.md_tolerance, remez);
}else{
std::cout<<GridLogMessage << "Using same rational approximations for MD as for action evaluation" << std::endl;
ApproxPowerMD = ApproxPowerAction;
ApproxNegPowerMD = ApproxNegPowerAction;
for(int i=0;i<ApproxPowerMD.tolerances.size();i++)
ApproxNegPowerMD.tolerances[i] = ApproxPowerMD.tolerances[i] = param.md_tolerance; //used for multishift
ApproxHalfPowerMD = ApproxHalfPowerAction;
ApproxNegHalfPowerMD = ApproxNegHalfPowerAction;
for(int i=0;i<ApproxPowerMD.tolerances.size();i++)
ApproxNegHalfPowerMD.tolerances[i] = ApproxHalfPowerMD.tolerances[i] = param.md_tolerance;
}
std::cout<<GridLogMessage << action_name() << " initialize: complete" << std::endl;
};
virtual std::string action_name(){return "GeneralEvenOddRatioRationalPseudoFermionAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] Power : 1/" << param.inv_pow << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] High :" << param.hi << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Tolerance (Action) :" << param.action_tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree (Action) :" << param.action_degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Tolerance (MD) :" << param.md_tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree (MD) :" << param.md_degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl;
return sstream.str();
}
//Access the fermion field
const FermionField &getPhiOdd() const{ return PhiOdd; }
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
std::cout<<GridLogMessage << action_name() << " refresh: starting" << std::endl;
FermionField eta(NumOp.FermionGrid());
// P(eta) \propto e^{- eta^dag eta}
//
// The gaussian function draws from P(x) \propto e^{- x^2 / 2 } [i.e. sigma=1]
// Thus eta = x/sqrt{2} = x * sqrt(1/2)
RealD scale = std::sqrt(0.5);
gaussian(pRNG,eta); eta=eta*scale;
refresh(U,eta);
}
//Allow for manual specification of random field for testing
void refresh(const GaugeField &U, const FermionField &eta) {
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//
// P(phi) = e^{- phi^dag (VdagV)^1/(2*inv_pow) (MdagM)^-1/inv_pow (VdagV)^1/(2*inv_pow) phi}
// = e^{- phi^dag (VdagV)^1/(2*inv_pow) (MdagM)^-1/(2*inv_pow) (MdagM)^-1/(2*inv_pow) (VdagV)^1/(2*inv_pow) phi}
//
// Phi = (VdagV)^-1/(2*inv_pow) Mdag^{1/(2*inv_pow)} eta
std::cout<<GridLogMessage << action_name() << " refresh: starting" << std::endl;
FermionField etaOdd (NumOp.FermionRedBlackGrid());
FermionField etaEven(NumOp.FermionRedBlackGrid());
FermionField tmp(NumOp.FermionRedBlackGrid());
pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta);
ImportGauge(U);
// MdagM^1/(2*inv_pow) eta
std::cout<<GridLogMessage << action_name() << " refresh: doing (M^dag M)^{1/" << 2*param.inv_pow << "} eta" << std::endl;
multiShiftInverse(Denominator, ApproxHalfPowerAction, param.MaxIter, etaOdd, tmp);
// VdagV^-1/(2*inv_pow) MdagM^1/(2*inv_pow) eta
std::cout<<GridLogMessage << action_name() << " refresh: doing (V^dag V)^{-1/" << 2*param.inv_pow << "} ( (M^dag M)^{1/" << 2*param.inv_pow << "} eta)" << std::endl;
multiShiftInverse(Numerator, ApproxNegHalfPowerAction, param.MaxIter, tmp, PhiOdd);
assert(NumOp.ConstEE() == 1);
assert(DenOp.ConstEE() == 1);
PhiEven = Zero();
std::cout<<GridLogMessage << action_name() << " refresh: starting" << std::endl;
};
//////////////////////////////////////////////////////
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
std::cout<<GridLogMessage << action_name() << " compute action: starting" << std::endl;
ImportGauge(U);
FermionField X(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid());
// VdagV^1/(2*inv_pow) Phi
std::cout<<GridLogMessage << action_name() << " compute action: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
multiShiftInverse(Numerator, ApproxHalfPowerAction, param.MaxIter, PhiOdd,X);
// MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi
std::cout<<GridLogMessage << action_name() << " compute action: doing (M^dag M)^{-1/" << 2*param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
multiShiftInverse(Denominator, ApproxNegHalfPowerAction, param.MaxIter, X,Y);
// Randomly apply rational bounds checks.
int rcheck = rand();
auto grid = NumOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( param.BoundsCheckFreq != 0 && (r % param.BoundsCheckFreq)==0 ) {
std::cout<<GridLogMessage << action_name() << " compute action: doing bounds check" << std::endl;
FermionField gauss(NumOp.FermionRedBlackGrid());
gauss = PhiOdd;
SchurDifferentiableOperator<Impl> MdagM(DenOp);
std::cout<<GridLogMessage << action_name() << " compute action: checking high bounds" << std::endl;
HighBoundCheck(MdagM,gauss,param.hi);
std::cout<<GridLogMessage << action_name() << " compute action: full approximation" << std::endl;
InversePowerBoundsCheck(param.inv_pow,param.MaxIter,param.action_tolerance*100,MdagM,gauss,ApproxNegPowerAction);
std::cout<<GridLogMessage << action_name() << " compute action: bounds check complete" << std::endl;
}
// Phidag VdagV^1/(2*inv_pow) MdagM^-1/(2*inv_pow) MdagM^-1/(2*inv_pow) VdagV^1/(2*inv_pow) Phi
RealD action = norm2(Y);
std::cout<<GridLogMessage << action_name() << " compute action: complete" << std::endl;
return action;
};
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//
// Here, M is some 5D operator and V is the Pauli-Villars field
// N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
//
// Need
// dS_f/dU = chi^dag d[P/Q] N/D P/Q chi
// + chi^dag P/Q d[N/D] P/Q chi
// + chi^dag P/Q N/D d[P/Q] chi
//
// P/Q is expressed as partial fraction expansion:
//
// a0 + \sum_k ak/(V^dagV + bk)
//
// d[P/Q] is then
//
// \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
//
// and similar for N/D.
//
// Need
// MpvPhi_k = [Vdag V + bk]^{-1} chi
// MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi
//
// MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
// MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
//
// MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
//
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
std::cout<<GridLogMessage << action_name() << " deriv: starting" << std::endl;
const int n_f = ApproxNegPowerMD.poles.size();
const int n_pv = ApproxHalfPowerMD.poles.size();
std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid());
std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid());
std::vector<FermionField> MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid());
FermionField MpvPhi(NumOp.FermionRedBlackGrid());
FermionField MfMpvPhi(NumOp.FermionRedBlackGrid());
FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid());
GaugeField tmp(NumOp.GaugeGrid());
ImportGauge(U);
std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} Phi" << std::endl;
multiShiftInverse(Numerator, ApproxHalfPowerMD, param.MaxIter, PhiOdd,MpvPhi_k,MpvPhi);
std::cout<<GridLogMessage << action_name() << " deriv: doing (M^dag M)^{-1/" << param.inv_pow << "} ( (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
multiShiftInverse(Denominator, ApproxNegPowerMD, param.MaxIter, MpvPhi,MfMpvPhi_k,MfMpvPhi);
std::cout<<GridLogMessage << action_name() << " deriv: doing (V^dag V)^{1/" << 2*param.inv_pow << "} ( (M^dag M)^{-1/" << param.inv_pow << "} (V^dag V)^{1/" << 2*param.inv_pow << "} Phi)" << std::endl;
multiShiftInverse(Numerator, ApproxHalfPowerMD, param.MaxIter, MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
SchurDifferentiableOperator<Impl> MdagM(DenOp);
SchurDifferentiableOperator<Impl> VdagV(NumOp);
RealD ak;
dSdU = Zero();
// With these building blocks
//
// dS/dU =
// \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1)
// + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2)
// -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3)
//(1)
std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (1)" << std::endl;
for(int k=0;k<n_f;k++){
ak = ApproxNegPowerMD.residues[k];
MdagM.Mpc(MfMpvPhi_k[k],Y);
MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp;
MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp;
}
//(2)
//(3)
std::cout<<GridLogMessage << action_name() << " deriv: doing dS/dU part (2)+(3)" << std::endl;
for(int k=0;k<n_pv;k++){
ak = ApproxHalfPowerMD.residues[k];
VdagV.Mpc(MpvPhi_k[k],Y);
VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp;
VdagV.MpcDeriv (tmp,Y,MpvMfMpvPhi_k[k]); dSdU=dSdU+ak*tmp;
VdagV.Mpc(MpvMfMpvPhi_k[k],Y); // V as we take Ydag
VdagV.MpcDeriv (tmp,Y, MpvPhi_k[k]); dSdU=dSdU+ak*tmp;
VdagV.MpcDagDeriv(tmp,MpvPhi_k[k], Y); dSdU=dSdU+ak*tmp;
}
//dSdU = Ta(dSdU);
std::cout<<GridLogMessage << action_name() << " deriv: complete" << std::endl;
};
};
NAMESPACE_END(Grid);
#endif

93
Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h Normal file
View File

@@ -0,0 +1,93 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <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_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
NAMESPACE_BEGIN(Grid);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Generic rational approximation for ratios of operators utilizing the mixed precision multishift algorithm
// cf. GeneralEvenOddRational.h for details
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class ImplD, class ImplF>
class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> {
private:
typedef typename ImplD::FermionField FermionFieldD;
typedef typename ImplF::FermionField FermionFieldF;
FermionOperator<ImplD> & NumOpD;
FermionOperator<ImplD> & DenOpD;
FermionOperator<ImplF> & NumOpF;
FermionOperator<ImplF> & DenOpF;
Integer ReliableUpdateFreq;
protected:
//Allow derived classes to override the multishift CG
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){
SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
msCG(schurOpD, in, out);
}
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &out){
SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
msCG(schurOpD, in, out_elems, out);
}
//Allow derived classes to override the gauge import
virtual void ImportGauge(const typename ImplD::GaugeField &Ud){
typename ImplF::GaugeField Uf(NumOpF.GaugeGrid());
precisionChange(Uf, Ud);
NumOpD.ImportGauge(Ud);
DenOpD.ImportGauge(Ud);
NumOpF.ImportGauge(Uf);
DenOpF.ImportGauge(Uf);
}
public:
GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD> &_NumOpD, FermionOperator<ImplD> &_DenOpD,
FermionOperator<ImplF> &_NumOpF, FermionOperator<ImplF> &_DenOpF,
const RationalActionParams & p, Integer _ReliableUpdateFreq
) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p),
ReliableUpdateFreq(_ReliableUpdateFreq), NumOpD(_NumOpD), DenOpD(_DenOpD), NumOpF(_NumOpF), DenOpF(_DenOpF){}
virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";}
};
NAMESPACE_END(Grid);
#endif

263
Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h
View File

@@ -40,262 +40,31 @@ NAMESPACE_BEGIN(Grid);
// Here N/D \sim R_{-1/2} ~ (M^dagM)^{-1/2}
template<class Impl>
class OneFlavourEvenOddRatioRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
class OneFlavourEvenOddRatioRationalPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<Impl> {
public:
INHERIT_IMPL_TYPES(Impl);
typedef OneFlavourRationalParams Params;
Params param;
MultiShiftFunction PowerHalf ;
MultiShiftFunction PowerNegHalf;
MultiShiftFunction PowerQuarter;
MultiShiftFunction PowerNegQuarter;
MultiShiftFunction MDPowerNegHalf;
MultiShiftFunction MDPowerQuarter;
private:
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
FermionField PhiEven; // the pseudo fermion field for this trajectory
FermionField PhiOdd; // the pseudo fermion field for this trajectory
FermionField Noise; // spare noise field for bounds check
static RationalActionParams transcribe(const Params &in){
RationalActionParams out;
out.inv_pow = 2;
out.lo = in.lo;
out.hi = in.hi;
out.MaxIter = in.MaxIter;
out.action_tolerance = out.md_tolerance = in.tolerance;
out.action_degree = out.md_degree = in.degree;
out.precision = in.precision;
out.BoundsCheckFreq = in.BoundsCheckFreq;
return out;
}
public:
OneFlavourEvenOddRatioRationalPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
Params & p
const Params & p
) :
NumOp(_NumOp),
DenOp(_DenOp),
PhiOdd (_NumOp.FermionRedBlackGrid()),
PhiEven(_NumOp.FermionRedBlackGrid()),
Noise(_NumOp.FermionRedBlackGrid()),
param(p)
{
AlgRemez remez(param.lo,param.hi,param.precision);
GeneralEvenOddRatioRationalPseudoFermionAction<Impl>(_NumOp, _DenOp, transcribe(p)){}
// MdagM^(+- 1/2)
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(param.degree,1,2);
PowerHalf.Init(remez,param.tolerance,false);
PowerNegHalf.Init(remez,param.tolerance,true);
MDPowerNegHalf.Init(remez,param.mdtolerance,true);
// MdagM^(+- 1/4)
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl;
remez.generateApprox(param.degree,1,4);
PowerQuarter.Init(remez,param.tolerance,false);
MDPowerQuarter.Init(remez,param.mdtolerance,false);
PowerNegQuarter.Init(remez,param.tolerance,true);
};
virtual std::string action_name(){
std::stringstream sstream;
sstream<< "OneFlavourEvenOddRatioRationalPseudoFermionAction det("<< DenOp.Mass() << ") / det("<<NumOp.Mass()<<")";
return sstream.str();
}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] Low :" << param.lo << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] High :" << param.hi << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Max iterations :" << param.MaxIter << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Tolerance :" << param.tolerance << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Degree :" << param.degree << std::endl;
sstream << GridLogMessage << "["<<action_name()<<"] Precision :" << param.precision << std::endl;
return sstream.str();
}
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//
// P(phi) = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/2 (VdagV)^1/4 phi}
// = e^{- phi^dag (VdagV)^1/4 (MdagM)^-1/4 (MdagM)^-1/4 (VdagV)^1/4 phi}
//
// Phi = (VdagV)^-1/4 Mdag^{1/4} 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).
RealD scale = std::sqrt(0.5);
FermionField eta(NumOp.FermionGrid());
FermionField etaOdd (NumOp.FermionRedBlackGrid());
FermionField etaEven(NumOp.FermionRedBlackGrid());
FermionField tmp(NumOp.FermionRedBlackGrid());
gaussian(pRNG,eta); eta=eta*scale;
pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta);
Noise = etaOdd;
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
// MdagM^1/4 eta
SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerQuarter);
msCG_M(MdagM,etaOdd,tmp);
// VdagV^-1/4 MdagM^1/4 eta
SchurDifferentiableOperator<Impl> VdagV(NumOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerNegQuarter);
msCG_V(VdagV,tmp,PhiOdd);
assert(NumOp.ConstEE() == 1);
assert(DenOp.ConstEE() == 1);
PhiEven = Zero();
};
//////////////////////////////////////////////////////
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField X(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid());
// VdagV^1/4 Phi
SchurDifferentiableOperator<Impl> VdagV(NumOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,PowerQuarter);
msCG_V(VdagV,PhiOdd,X);
// MdagM^-1/4 VdagV^1/4 Phi
SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,PowerNegQuarter);
msCG_M(MdagM,X,Y);
// Randomly apply rational bounds checks.
auto grid = NumOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( (r%param.BoundsCheckFreq)==0 ) {
FermionField gauss(NumOp.FermionRedBlackGrid());
gauss = Noise;
HighBoundCheck(MdagM,gauss,param.hi);
InverseSqrtBoundsCheck(param.MaxIter,param.tolerance*100,MdagM,gauss,PowerNegHalf);
ChebyBoundsCheck(MdagM,Noise,param.lo,param.hi);
}
// Phidag VdagV^1/4 MdagM^-1/4 MdagM^-1/4 VdagV^1/4 Phi
RealD action = norm2(Y);
return action;
};
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
//
// Here, M is some 5D operator and V is the Pauli-Villars field
// N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the denom term
//
// Need
// dS_f/dU = chi^dag d[P/Q] N/D P/Q chi
// + chi^dag P/Q d[N/D] P/Q chi
// + chi^dag P/Q N/D d[P/Q] chi
//
// P/Q is expressed as partial fraction expansion:
//
// a0 + \sum_k ak/(V^dagV + bk)
//
// d[P/Q] is then
//
// \sum_k -ak [V^dagV+bk]^{-1} [ dV^dag V + V^dag dV ] [V^dag V + bk]^{-1}
//
// and similar for N/D.
//
// Need
// MpvPhi_k = [Vdag V + bk]^{-1} chi
// MpvPhi = {a0 + \sum_k ak [Vdag V + bk]^{-1} }chi
//
// MfMpvPhi_k = [MdagM+bk]^{-1} MpvPhi
// MfMpvPhi = {a0 + \sum_k ak [Mdag M + bk]^{-1} } MpvPhi
//
// MpvMfMpvPhi_k = [Vdag V + bk]^{-1} MfMpvchi
//
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
const int n_f = MDPowerNegHalf.poles.size();
const int n_pv = MDPowerQuarter.poles.size();
std::vector<FermionField> MpvPhi_k (n_pv,NumOp.FermionRedBlackGrid());
std::vector<FermionField> MpvMfMpvPhi_k(n_pv,NumOp.FermionRedBlackGrid());
std::vector<FermionField> MfMpvPhi_k (n_f ,NumOp.FermionRedBlackGrid());
FermionField MpvPhi(NumOp.FermionRedBlackGrid());
FermionField MfMpvPhi(NumOp.FermionRedBlackGrid());
FermionField MpvMfMpvPhi(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid());
GaugeField tmp(NumOp.GaugeGrid());
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> VdagV(NumOp);
SchurDifferentiableOperator<Impl> MdagM(DenOp);
ConjugateGradientMultiShift<FermionField> msCG_V(param.MaxIter,MDPowerQuarter);
ConjugateGradientMultiShift<FermionField> msCG_M(param.MaxIter,MDPowerNegHalf);
msCG_V(VdagV,PhiOdd,MpvPhi_k,MpvPhi);
msCG_M(MdagM,MpvPhi,MfMpvPhi_k,MfMpvPhi);
msCG_V(VdagV,MfMpvPhi,MpvMfMpvPhi_k,MpvMfMpvPhi);
RealD ak;
dSdU = Zero();
// With these building blocks
//
// dS/dU =
// \sum_k -ak MfMpvPhi_k^dag [ dM^dag M + M^dag dM ] MfMpvPhi_k (1)
// + \sum_k -ak MpvMfMpvPhi_k^\dag [ dV^dag V + V^dag dV ] MpvPhi_k (2)
// -ak MpvPhi_k^dag [ dV^dag V + V^dag dV ] MpvMfMpvPhi_k (3)
//(1)
for(int k=0;k<n_f;k++){
ak = MDPowerNegHalf.residues[k];
MdagM.Mpc(MfMpvPhi_k[k],Y);
MdagM.MpcDagDeriv(tmp , MfMpvPhi_k[k], Y ); dSdU=dSdU+ak*tmp;
MdagM.MpcDeriv(tmp , Y, MfMpvPhi_k[k] ); dSdU=dSdU+ak*tmp;
}
//(2)
//(3)
for(int k=0;k<n_pv;k++){
ak = MDPowerQuarter.residues[k];
VdagV.Mpc(MpvPhi_k[k],Y);
VdagV.MpcDagDeriv(tmp,MpvMfMpvPhi_k[k],Y); dSdU=dSdU+ak*tmp;
VdagV.MpcDeriv (tmp,Y,MpvMfMpvPhi_k[k]); dSdU=dSdU+ak*tmp;
VdagV.Mpc(MpvMfMpvPhi_k[k],Y); // V as we take Ydag
VdagV.MpcDeriv (tmp,Y, MpvPhi_k[k]); dSdU=dSdU+ak*tmp;
VdagV.MpcDagDeriv(tmp,MpvPhi_k[k], Y); dSdU=dSdU+ak*tmp;
}
//dSdU = Ta(dSdU);
};
virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";}
};
NAMESPACE_END(Grid);

2
Grid/qcd/action/pseudofermion/PseudoFermion.h
View File

@@ -40,6 +40,8 @@ directory
#include <Grid/qcd/action/pseudofermion/OneFlavourRational.h>
#include <Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h>
#include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
#include <Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h>
#include <Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h>
#include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h>
#include <Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h>

28
Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
View File

@@ -88,15 +88,9 @@ NAMESPACE_BEGIN(Grid);
}
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
const FermionField &getPhiOdd() const{ return PhiOdd; }
// P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi}
//
// NumOp == V
// DenOp == M
//
// Take phi_o = Vpcdag^{-1} Mpcdag eta_o ; eta_o = Mpcdag^{-1} Vpcdag Phi
//
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.
@@ -104,12 +98,22 @@ NAMESPACE_BEGIN(Grid);
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());
gaussian(pRNG,eta);
pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta);
@@ -128,10 +132,6 @@ NAMESPACE_BEGIN(Grid);
// Even det factors
DenOp.MooeeDag(etaEven,tmp);
NumOp.MooeeInvDag(tmp,PhiEven);
PhiOdd =PhiOdd*scale;
PhiEven=PhiEven*scale;
};
//////////////////////////////////////////////////////

12
Grid/qcd/hmc/GenericHMCrunner.h
View File

@@ -151,12 +151,22 @@ public:
Resources.GetCheckPointer()->CheckpointRestore(Parameters.StartTrajectory, U,
Resources.GetSerialRNG(),
Resources.GetParallelRNG());
} else if (Parameters.StartingType == "CheckpointStartReseed") {
// Same as CheckpointRestart but reseed the RNGs using the fixed integer seeding used for ColdStart and HotStart
// Useful for creating new evolution streams from an existing stream
// WARNING: Unfortunately because the checkpointer doesn't presently allow us to separately restore the RNG and gauge fields we have to load
// an existing RNG checkpoint first; make sure one is available and named correctly
Resources.GetCheckPointer()->CheckpointRestore(Parameters.StartTrajectory, U,
Resources.GetSerialRNG(),
Resources.GetParallelRNG());
Resources.SeedFixedIntegers();
} else {
// others
std::cout << GridLogError << "Unrecognized StartingType\n";
std::cout
<< GridLogError
<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart, CheckpointStartReseed]\n";
exit(1);
}
}

2
Grid/qcd/hmc/HMCModules.h
View File

@@ -80,7 +80,9 @@ public:
std::cout << GridLogError << "Seeds not initialized" << std::endl;
exit(1);
}
std::cout << GridLogMessage << "Reseeding serial RNG with seed vector " << SerialSeeds << std::endl;
sRNG_.SeedFixedIntegers(SerialSeeds);
std::cout << GridLogMessage << "Reseeding parallel RNG with seed vector " << ParallelSeeds << std::endl;
pRNG_->SeedFixedIntegers(ParallelSeeds);
}
};

6
Grid/qcd/hmc/integrators/Integrator.h
View File

@@ -334,15 +334,19 @@ public:
void refresh(Field& U, GridSerialRNG & sRNG, GridParallelRNG& pRNG)
{
assert(P.Grid() == U.Grid());
std::cout << GridLogIntegrator << "Integrator refresh\n";
std::cout << GridLogIntegrator << "Integrator refresh" << std::endl;
std::cout << GridLogIntegrator << "Generating momentum" << std::endl;
FieldImplementation::generate_momenta(P, sRNG, pRNG);
// Update the smeared fields, can be implemented as observer
// necessary to keep the fields updated even after a reject
// of the Metropolis
std::cout << GridLogIntegrator << "Updating smeared fields" << std::endl;
Smearer.set_Field(U);
// Set the (eventual) representations gauge fields
std::cout << GridLogIntegrator << "Updating representations" << std::endl;
Representations.update(U);
// The Smearer is attached to a pointer of the gauge field

446
tests/forces/Test_gpdwf_force_1f_2f.cc Normal file
View File

@@ -0,0 +1,446 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./forces/Test_gpdwf_force_1f_2f.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
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 */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
//Here we test the G-parity action and force between the 1f (doubled-lattice) and 2f approaches
void copyConjGauge(LatticeGaugeFieldD &Umu_1f, const LatticeGaugeFieldD &Umu_2f, const int nu){
GridBase* UGrid_2f = Umu_2f.Grid();
GridBase* UGrid_1f = Umu_1f.Grid();
Replicate(Umu_2f,Umu_1f);
int L_2f = UGrid_2f->FullDimensions()[nu];
int L_1f = UGrid_1f->FullDimensions()[nu];
assert(L_1f == 2 * L_2f);
//Coordinate grid for reference
LatticeInteger xcoor_1f(UGrid_1f);
LatticeCoordinate(xcoor_1f,nu);
//Copy-conjugate the gauge field
//First C-shift the lattice by Lx/2
{
LatticeGaugeField Umu_shift = conjugate( Cshift(Umu_1f,nu,L_2f) );
Umu_1f = where( xcoor_1f >= Integer(L_2f), Umu_shift, Umu_1f );
//We use the in built APBC
//Make the gauge field antiperiodic in nu-direction
//decltype(PeekIndex<LorentzIndex>(Umu_1f,nu)) Unu(UGrid_1f);
//Unu = PeekIndex<LorentzIndex>(Umu_1f,nu);
//Unu = where(xcoor_1f == Integer(2*L_2f-1), -Unu, Unu);
//PokeIndex<LorentzIndex>(Umu_1f,Unu,nu);
}
}
template<typename FermionField2f, typename FermionField1f>
void convertFermion1f_from_2f(FermionField1f &out_1f, const FermionField2f &in_2f, const int nu, bool is_4d){
GridBase* FGrid_1f = out_1f.Grid();
GridBase* FGrid_2f = in_2f.Grid();
int nuoff = is_4d ? 0 : 1; //s in 0 direction
int L_2f = FGrid_2f->FullDimensions()[nu+nuoff];
int L_1f = FGrid_1f->FullDimensions()[nu+nuoff];
assert(L_1f == 2 * L_2f);
auto in_f0_2fgrid = PeekIndex<GparityFlavourIndex>(in_2f,0); //flavor 0 on 2f Grid
FermionField1f in_f0_1fgrid(FGrid_1f);
Replicate(in_f0_2fgrid, in_f0_1fgrid); //has flavor 0 on both halves
auto in_f1_2fgrid = PeekIndex<GparityFlavourIndex>(in_2f,1); //flavor 1 on 2f Grid
FermionField1f in_f1_1fgrid(FGrid_1f);
Replicate(in_f1_2fgrid, in_f1_1fgrid); //has flavor 1 on both halves
LatticeInteger xcoor_1f(FGrid_1f);
LatticeCoordinate(xcoor_1f,nu+nuoff);
out_1f = where(xcoor_1f < L_2f, in_f0_1fgrid, in_f1_1fgrid);
}
template<typename GparityAction, typename StandardAction>
class RatioActionSetupBase{
protected:
TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplD> *pf_1f;
TwoFlavourEvenOddRatioPseudoFermionAction<GparityWilsonImplD> *pf_2f;
GparityAction* action_2f;
GparityAction* action_PV_2f;
StandardAction* action_1f;
StandardAction* action_PV_1f;
ConjugateGradient<typename StandardAction::FermionField> CG_1f;
ConjugateGradient<typename GparityAction::FermionField> CG_2f;
RatioActionSetupBase(): CG_1f(1.0e-8,10000), CG_2f(1.0e-8,10000){}
void setupPseudofermion(){
pf_1f = new TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplD>(*action_PV_1f, *action_1f, CG_1f, CG_1f);
pf_2f = new TwoFlavourEvenOddRatioPseudoFermionAction<GparityWilsonImplD>(*action_PV_2f, *action_2f, CG_2f, CG_2f);
}
public:
GparityAction & action2f(){ return *action_2f; }
StandardAction & action1f(){ return *action_1f; }
void refreshAction(LatticeGaugeField &Umu_2f, typename GparityAction::FermionField &eta_2f,
LatticeGaugeField &Umu_1f, typename StandardAction::FermionField &eta_1f){
pf_1f->refresh(Umu_1f, eta_1f);
pf_2f->refresh(Umu_2f, eta_2f);
//Compare PhiOdd
RealD norm_1f = norm2(pf_1f->getPhiOdd());
RealD norm_2f = norm2(pf_2f->getPhiOdd());
std::cout << "Test PhiOdd 2f: " << norm_2f << " 1f: " << norm_1f << std::endl;
}
void computeAction(RealD &S_2f, RealD &S_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){
S_1f = pf_1f->S(Umu_1f);
S_2f = pf_2f->S(Umu_2f);
}
void computeDeriv(LatticeGaugeField &deriv_2f, LatticeGaugeField &deriv_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){
pf_1f->deriv(Umu_1f, deriv_1f);
pf_2f->deriv(Umu_2f, deriv_2f);
}
};
template<typename GparityAction, typename StandardAction>
struct setupAction{};
template<>
struct setupAction<GparityWilsonTMFermionD, WilsonTMFermionD>: public RatioActionSetupBase<GparityWilsonTMFermionD, WilsonTMFermionD>{
typedef GparityWilsonTMFermionD GparityAction;
typedef WilsonTMFermionD StandardAction;
setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f,
GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f,
LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): RatioActionSetupBase(){
RealD mass=-1.8;
//Use same DSDR twists as https://arxiv.org/pdf/1208.4412.pdf
RealD epsilon_f = 0.02; //numerator (in determinant)
RealD epsilon_b = 0.5;
std::vector<int> twists(Nd,0);
twists[nu] = 1; //GPBC in y
twists[3] = 1; //APBC
GparityAction::ImplParams params_2f; params_2f.twists = twists;
action_2f = new GparityWilsonTMFermionD(Umu_2f,*UGrid_2f,*UrbGrid_2f, mass, epsilon_f, params_2f);
action_PV_2f = new GparityWilsonTMFermionD(Umu_2f,*UGrid_2f,*UrbGrid_2f, mass, epsilon_b, params_2f);
DomainWallFermionD::ImplParams params_1f;
params_1f.boundary_phases[nu] = -1;
params_1f.boundary_phases[3] = -1;
action_1f = new WilsonTMFermionD(Umu_1f,*UGrid_1f,*UrbGrid_1f, mass, epsilon_f, params_1f);
action_PV_1f = new WilsonTMFermionD(Umu_1f,*UGrid_1f,*UrbGrid_1f, mass, epsilon_b, params_1f);
setupPseudofermion();
}
static bool is4d(){ return true; }
};
template<>
struct setupAction<GparityDomainWallFermionD, DomainWallFermionD>: public RatioActionSetupBase<GparityDomainWallFermionD, DomainWallFermionD>{
typedef GparityDomainWallFermionD GparityAction;
typedef DomainWallFermionD StandardAction;
setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f,
GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f,
LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): RatioActionSetupBase(){
RealD mass=0.01;
RealD M5=1.8;
std::vector<int> twists(Nd,0);
twists[nu] = 1; //GPBC in y
twists[3] = 1; //APBC
GparityDomainWallFermionD::ImplParams params_2f; params_2f.twists = twists;
action_2f = new GparityDomainWallFermionD(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5,params_2f);
action_PV_2f = new GparityDomainWallFermionD(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,1.0,M5,params_2f);
DomainWallFermionD::ImplParams params_1f;
params_1f.boundary_phases[nu] = -1;
params_1f.boundary_phases[3] = -1;
action_1f = new DomainWallFermionD(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5,params_1f);
action_PV_1f = new DomainWallFermionD(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,1.0,M5,params_1f);
setupPseudofermion();
}
static bool is4d(){ return false; }
};
//For EOFA we need a different pseudofermion type
template<>
struct setupAction<GparityDomainWallEOFAFermionD, DomainWallEOFAFermionD>{
typedef GparityDomainWallEOFAFermionD GparityAction;
typedef DomainWallEOFAFermionD StandardAction;
ExactOneFlavourRatioPseudoFermionAction<WilsonImplD> *pf_1f;
ExactOneFlavourRatioPseudoFermionAction<GparityWilsonImplD> *pf_2f;
GparityAction* action_2f;
GparityAction* action_PV_2f;
StandardAction* action_1f;
StandardAction* action_PV_1f;
ConjugateGradient<typename StandardAction::FermionField> CG_1f;
ConjugateGradient<typename GparityAction::FermionField> CG_2f;
public:
GparityAction & action2f(){ return *action_2f; }
StandardAction & action1f(){ return *action_1f; }
void refreshAction(LatticeGaugeField &Umu_2f, typename GparityAction::FermionField &eta_2f,
LatticeGaugeField &Umu_1f, typename StandardAction::FermionField &eta_1f){
pf_1f->refresh(Umu_1f, eta_1f);
pf_2f->refresh(Umu_2f, eta_2f);
//Compare PhiOdd
RealD norm_1f = norm2(pf_1f->getPhi());
RealD norm_2f = norm2(pf_2f->getPhi());
std::cout << "Test Phi 2f: " << norm_2f << " 1f: " << norm_1f << std::endl;
}
void computeAction(RealD &S_2f, RealD &S_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){
S_1f = pf_1f->S(Umu_1f);
S_2f = pf_2f->S(Umu_2f);
}
void computeDeriv(LatticeGaugeField &deriv_2f, LatticeGaugeField &deriv_1f, LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f){
pf_1f->deriv(Umu_1f, deriv_1f);
pf_2f->deriv(Umu_2f, deriv_2f);
}
setupAction(GridCartesian* UGrid_2f, GridRedBlackCartesian* UrbGrid_2f, GridCartesian* FGrid_2f, GridRedBlackCartesian* FrbGrid_2f,
GridCartesian* UGrid_1f, GridRedBlackCartesian* UrbGrid_1f, GridCartesian* FGrid_1f, GridRedBlackCartesian* FrbGrid_1f,
LatticeGaugeField &Umu_2f, LatticeGaugeField &Umu_1f, int nu): CG_1f(1.0e-8,10000), CG_2f(1.0e-8,10000){
RealD mass=0.01;
RealD M5=1.8;
std::vector<int> twists(Nd,0);
twists[nu] = 1; //GPBC in y
twists[3] = 1; //APBC
GparityAction::ImplParams params_2f; params_2f.twists = twists;
action_2f = new GparityAction(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f, mass, mass, 1.0, 0.0, -1, M5, params_2f);
action_PV_2f = new GparityAction(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f, 1.0, mass, 1.0, -1.0, 1, M5, params_2f); //cf Test_dwf_gpforce_eofa.cc
StandardAction::ImplParams params_1f;
params_1f.boundary_phases[nu] = -1;
params_1f.boundary_phases[3] = -1;
action_1f = new StandardAction(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f, mass, mass, 1.0, 0.0, -1, M5, params_1f);
action_PV_1f = new StandardAction(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f, 1.0, mass, 1.0, -1.0, 1, M5, params_1f);
OneFlavourRationalParams RationalParams(0.95, 100.0, 5000, 1.0e-12, 12);
pf_1f = new ExactOneFlavourRatioPseudoFermionAction<WilsonImplD>(*action_1f, *action_PV_1f, CG_1f, CG_1f, CG_1f, CG_1f, CG_1f, RationalParams, true);
pf_2f = new ExactOneFlavourRatioPseudoFermionAction<GparityWilsonImplD>(*action_2f, *action_PV_2f, CG_2f, CG_2f, CG_2f, CG_2f, CG_2f, RationalParams, true);
}
static bool is4d(){ return false; }
};
template<typename GparityAction, typename StandardAction>
void runTest(int argc, char** argv){
Grid_init(&argc,&argv);
const int nu = 1;
Coordinate latt_2f = GridDefaultLatt();
Coordinate latt_1f = latt_2f;
latt_1f[nu] *= 2;
Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
const int Ls=8;
GridCartesian * UGrid_1f = SpaceTimeGrid::makeFourDimGrid(latt_1f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_1f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_1f);
GridCartesian * FGrid_1f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_1f);
GridRedBlackCartesian * FrbGrid_1f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_1f);
GridCartesian * UGrid_2f = SpaceTimeGrid::makeFourDimGrid(latt_2f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_2f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_2f);
GridCartesian * FGrid_2f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_2f);
GridRedBlackCartesian * FrbGrid_2f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_2f);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu_2f(UGrid_2f);
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
LatticeGaugeField Umu_1f(UGrid_1f);
copyConjGauge(Umu_1f, Umu_2f, nu);
typedef typename GparityAction::FermionField GparityFermionField;
typedef typename StandardAction::FermionField StandardFermionField;
setupAction<GparityAction, StandardAction> setup(UGrid_2f, UrbGrid_2f, FGrid_2f, FrbGrid_2f,
UGrid_1f, UrbGrid_1f, FGrid_1f, FrbGrid_1f,
Umu_2f, Umu_1f, nu);
GridBase* FGrid_2f_a = setup.action2f().FermionGrid();
GridBase* FGrid_1f_a = setup.action1f().FermionGrid();
GridBase* FrbGrid_2f_a = setup.action2f().FermionRedBlackGrid();
GridBase* FrbGrid_1f_a = setup.action1f().FermionRedBlackGrid();
bool is_4d = setup.is4d();
//Check components by doing an inversion
{
setup.action2f().ImportGauge(Umu_2f);
setup.action1f().ImportGauge(Umu_1f);
GparityFermionField src_2f(FGrid_2f_a);
gaussian(is_4d ? RNG4_2f : RNG5_2f, src_2f);
StandardFermionField src_1f(FGrid_1f_a);
convertFermion1f_from_2f(src_1f, src_2f, nu, is_4d);
StandardFermionField src_o_1f(FrbGrid_1f_a);
StandardFermionField result_o_1f(FrbGrid_1f_a);
pickCheckerboard(Odd,src_o_1f,src_1f);
result_o_1f=Zero();
SchurDiagMooeeOperator<StandardAction,StandardFermionField> HermOpEO_1f(setup.action1f());
ConjugateGradient<StandardFermionField> CG_1f(1.0e-8,10000);
CG_1f(HermOpEO_1f,src_o_1f,result_o_1f);
GparityFermionField src_o_2f(FrbGrid_2f_a);
GparityFermionField result_o_2f(FrbGrid_2f_a);
pickCheckerboard(Odd,src_o_2f,src_2f);
result_o_2f=Zero();
SchurDiagMooeeOperator<GparityAction,GparityFermionField> HermOpEO_2f(setup.action2f());
ConjugateGradient<GparityFermionField> CG_2f(1.0e-8,10000);
CG_2f(HermOpEO_2f,src_o_2f,result_o_2f);
RealD norm_1f = norm2(result_o_1f);
RealD norm_2f = norm2(result_o_2f);
std::cout << "Test fermion inversion 2f: " << norm_2f << " 1f: " << norm_1f << std::endl;
}
//Generate eta
RealD scale = std::sqrt(0.5);
GparityFermionField eta_2f(FGrid_2f_a);
gaussian(is_4d ? RNG4_2f : RNG5_2f,eta_2f); eta_2f = eta_2f * scale;
StandardFermionField eta_1f(FGrid_1f_a);
convertFermion1f_from_2f(eta_1f, eta_2f, nu, is_4d);
setup.refreshAction(Umu_2f, eta_2f, Umu_1f, eta_1f);
//Initial action is just |eta^2|
RealD S_1f, S_2f;
setup.computeAction(S_2f, S_1f, Umu_2f, Umu_1f);
std::cout << "Test Initial action 2f: " << S_2f << " 1f: " << S_1f << " diff: " << S_2f - S_1f << std::endl;
//Do a random gauge field refresh
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
copyConjGauge(Umu_1f, Umu_2f, nu);
//Compute the action again
setup.computeAction(S_2f, S_1f, Umu_2f, Umu_1f);
std::cout << "Test Action after gauge field randomize 2f: " << S_2f << " 1f: " << S_1f << " diff: " << S_2f - S_1f << std::endl;
//Compute the derivative and test the conjugate relation
LatticeGaugeField deriv_2f(UGrid_2f);
LatticeGaugeField deriv_1f(UGrid_1f);
setup.computeDeriv(deriv_2f, deriv_1f, Umu_2f, Umu_1f);
//Have to combine the two forces on the 1f by symmetrizing under the complex conjugate
{
RealD norm2_pre = norm2(deriv_1f);
LatticeGaugeField deriv_1f_shift = conjugate( Cshift(deriv_1f, nu, latt_2f[nu]) );
deriv_1f = deriv_1f + deriv_1f_shift;
std::cout << "Test combine/symmetrize forces on 1f lattice, dS/dU : " << norm2_pre << " -> " << norm2(deriv_1f) << std::endl;
}
LatticeGaugeField deriv_1f_from_2f(UGrid_1f);
copyConjGauge(deriv_1f_from_2f, deriv_2f, nu);
std::cout << "Test copy-conj 2f dS/dU to obtain equivalent 1f force : " << norm2(deriv_2f) << " -> " << norm2(deriv_1f_from_2f) << std::endl;
LatticeGaugeField diff_deriv_1f = deriv_1f - deriv_1f_from_2f;
std::cout << "Test dS/dU 1f constructed from 2f derivative: " << norm2(deriv_1f_from_2f) << " dS/dU 1f actual: " << norm2(deriv_1f) << " Norm of difference: " << norm2(diff_deriv_1f) << std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}
int main (int argc, char ** argv)
{
std::string action = "DWF";
for(int i=1;i<argc;i++){
if(std::string(argv[i]) == "--action"){
action = argv[i+1];
}
}
if(action == "DWF"){
runTest<GparityDomainWallFermionD, DomainWallFermionD>(argc, argv);
}else if(action == "EOFA"){
runTest<GparityDomainWallEOFAFermionD, DomainWallEOFAFermionD>(argc, argv);
}else if(action == "DSDR"){
runTest<GparityWilsonTMFermionD, WilsonTMFermionD>(argc,argv);
}else{
assert(0);
}
}

42
tests/forces/Test_mobius_force_eofa.cc
View File

@@ -89,7 +89,49 @@ int main (int argc, char** argv)
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
//Check the rational approximation
{
RealD scale = std::sqrt(0.5);
LatticeFermion eta (Lop.FermionGrid());
gaussian(RNG5,eta); eta = eta * scale;
Meofa.refresh(U, eta);
//Phi = M^{-1/2} eta
//M is Hermitian
//(Phi, M Phi) = eta^\dagger M^{-1/2} M M^{-1/2} eta = eta^\dagger eta
LatticeFermion phi = Meofa.getPhi();
LatticeFermion Mphi(FGrid);
Meofa.Meofa(U, phi, Mphi);
std::cout << "Computing inner product" << std::endl;
ComplexD inner = innerProduct(phi, Mphi);
ComplexD test = inner - norm2(eta);
std::cout << "(phi, Mphi) - (eta,eta): " << test << " expect 0" << std::endl;
assert(test.real() < 1e-8);
assert(test.imag() < 1e-8);
//Another test is to use heatbath twice to apply M^{-1/2} to Phi then apply M
// M Phi'
//= M M^{-1/2} Phi
//= M M^{-1/2} M^{-1/2} eta
//= eta
Meofa.refresh(U, phi);
LatticeFermion phi2 = Meofa.getPhi();
LatticeFermion test2(FGrid);
Meofa.Meofa(U, phi2, test2);
test2 = test2 - eta;
RealD test2_norm = norm2(test2);
std::cout << "|M M^{-1/2} M^{-1/2} eta - eta|^2 = " << test2_norm << " expect 0" << std::endl;
assert( test2_norm < 1e-8 );
}
Meofa.refresh(U, sRNG, RNG5 );
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"

233
tests/forces/Test_mobius_gparity_eofa_mixed.cc Normal file
View File

@@ -0,0 +1,233 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/forces/Test_mobius_gparity_eofa_mixed.cc
Copyright (C) 2017
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: David Murphy <dmurphy@phys.columbia.edu>
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 */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
;
typedef GparityWilsonImplD FermionImplPolicyD;
typedef GparityMobiusEOFAFermionD FermionActionD;
typedef typename FermionActionD::FermionField FermionFieldD;
typedef GparityWilsonImplF FermionImplPolicyF;
typedef GparityMobiusEOFAFermionF FermionActionF;
typedef typename FermionActionF::FermionField FermionFieldF;
NAMESPACE_BEGIN(Grid);
template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF>
class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
public:
typedef typename FermionOperatorD::FermionField FieldD;
typedef typename FermionOperatorF::FermionField FieldF;
using OperatorFunction<FieldD>::operator();
RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations;
Integer MaxOuterIterations;
GridBase* SinglePrecGrid4; //Grid for single-precision fields
GridBase* SinglePrecGrid5; //Grid for single-precision fields
RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
FermionOperatorF &FermOpF;
FermionOperatorD &FermOpD;;
SchurOperatorF &LinOpF;
SchurOperatorD &LinOpD;
Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
MixedPrecisionConjugateGradientOperatorFunction(RealD tol,
Integer maxinnerit,
Integer maxouterit,
GridBase* _sp_grid4,
GridBase* _sp_grid5,
FermionOperatorF &_FermOpF,
FermionOperatorD &_FermOpD,
SchurOperatorF &_LinOpF,
SchurOperatorD &_LinOpD):
LinOpF(_LinOpF),
LinOpD(_LinOpD),
FermOpF(_FermOpF),
FermOpD(_FermOpD),
Tolerance(tol),
InnerTolerance(tol),
MaxInnerIterations(maxinnerit),
MaxOuterIterations(maxouterit),
SinglePrecGrid4(_sp_grid4),
SinglePrecGrid5(_sp_grid5),
OuterLoopNormMult(100.)
{
};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
precisionChange(FermOpF.Umu, FermOpD.Umu);
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
MPCG.InnerTolerance = InnerTolerance;
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
};
NAMESPACE_END(Grid);
int main (int argc, char** argv)
{
Grid_init(&argc, &argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate mpi_layout = GridDefaultMpi();
const int Ls = 8;
GridCartesian *UGridD = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian *UrbGridD = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridD);
GridCartesian *FGridD = SpaceTimeGrid::makeFiveDimGrid(Ls, UGridD);
GridRedBlackCartesian *FrbGridD = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGridD);
GridCartesian *UGridF = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian *UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF);
GridCartesian *FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls, UGridF);
GridRedBlackCartesian *FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGridF);
std::vector<int> seeds4({1,2,3,5});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGridD); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGridD); RNG4.SeedFixedIntegers(seeds4);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
LatticeGaugeFieldD Ud(UGridD);
SU<Nc>::HotConfiguration(RNG4,Ud);
LatticeGaugeFieldF Uf(UGridF);
precisionChange(Uf, Ud);
RealD b = 2.5;
RealD c = 1.5;
RealD mf = 0.01;
RealD mb = 1.0;
RealD M5 = 1.8;
FermionActionD::ImplParams params;
params.twists[0] = 1; //GPBC in X
params.twists[Nd-1] = 1; //APRD in T
std::vector<int> gtwists(4,0);
gtwists[0] = 1;
ConjugateGimplD::setDirections(gtwists);
FermionActionD LopD(Ud, *FGridD, *FrbGridD, *UGridD, *UrbGridD, mf, mf, mb, 0.0, -1, M5, b, c, params);
FermionActionD RopD(Ud, *FGridD, *FrbGridD, *UGridD, *UrbGridD, mb, mf, mb, -1.0, 1, M5, b, c, params);
FermionActionF LopF(Uf, *FGridF, *FrbGridF, *UGridF, *UrbGridF, mf, mf, mb, 0.0, -1, M5, b, c, params);
FermionActionF RopF(Uf, *FGridF, *FrbGridF, *UGridF, *UrbGridF, mb, mf, mb, -1.0, 1, M5, b, c, params);
OneFlavourRationalParams OFRp(0.95, 100.0, 5000, 1.0e-12, 12);
ConjugateGradient<FermionFieldD> CG(1.0e-10, 10000);
typedef SchurDiagMooeeOperator<FermionActionD,FermionFieldD> EOFAschuropD;
typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> EOFAschuropF;
EOFAschuropD linopL_D(LopD);
EOFAschuropD linopR_D(RopD);
EOFAschuropF linopL_F(LopF);
EOFAschuropF linopR_F(RopF);
typedef MixedPrecisionConjugateGradientOperatorFunction<FermionActionD, FermionActionF, EOFAschuropD, EOFAschuropF> EOFA_mxCG;
EOFA_mxCG MCG_L(1e-10, 10000, 1000, UGridF, FrbGridF, LopF, LopD, linopL_F, linopL_D);
MCG_L.InnerTolerance = 1e-5;
EOFA_mxCG MCG_R(1e-10, 10000, 1000, UGridF, FrbGridF, RopF, RopD, linopR_F, linopR_D);
MCG_R.InnerTolerance = 1e-5;
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicyD> MeofaD(LopD, RopD, CG, CG, CG, CG, CG, OFRp, true);
ExactOneFlavourRatioMixedPrecHeatbathPseudoFermionAction<FermionImplPolicyD, FermionImplPolicyF> MeofaMx(LopF, RopF, LopD, RopD, MCG_L, MCG_R, MCG_L, MCG_R, MCG_L, MCG_R, OFRp, true);
FermionFieldD eta(FGridD);
gaussian(RNG5, eta);
MeofaD.refresh(Ud, eta);
MeofaMx.refresh(Ud, eta);
FermionFieldD diff_phi(FGridD);
diff_phi = MeofaD.getPhi() - MeofaMx.getPhi();
RealD n = norm2(diff_phi);
std::cout << GridLogMessage << "Phi(double)=" << norm2(MeofaD.getPhi()) << " Phi(mixed)=" << norm2(MeofaMx.getPhi()) << " diff=" << n << std::endl;
assert(n < 1e-8);
RealD Sd = MeofaD.S(Ud);
RealD Smx = MeofaMx.S(Ud);
std::cout << GridLogMessage << "Initial action double=" << Sd << " mixed=" << Smx << " diff=" << Sd-Smx << std::endl;
assert(fabs(Sd-Smx) < 1e-6);
SU<Nc>::HotConfiguration(RNG4,Ud);
precisionChange(Uf, Ud);
Sd = MeofaD.S(Ud);
Smx = MeofaMx.S(Ud);
std::cout << GridLogMessage << "After randomizing U, action double=" << Sd << " mixed=" << Smx << " diff=" << Sd-Smx << std::endl;
assert(fabs(Sd-Smx) < 1e-6);
std::cout << GridLogMessage << "Done" << std::endl;
Grid_finalize();
}

257
tests/hmc/Test_action_dwf_gparity2fvs1f.cc Normal file
View File

@@ -0,0 +1,257 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hmc/Test_action_dwf_gparity2fvs1f.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
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 */
#include <Grid/Grid.h>
using namespace Grid;
template<typename FermionField2f, typename FermionField1f>
void copy2fTo1fFermionField(FermionField1f &out, const FermionField2f &in, int gpdir){
auto f0_halfgrid = PeekIndex<GparityFlavourIndex>(in,0); //on 2f Grid
FermionField1f f0_fullgrid_dbl(out.Grid());
Replicate(f0_halfgrid, f0_fullgrid_dbl); //double it up to live on the 1f Grid
auto f1_halfgrid = PeekIndex<GparityFlavourIndex>(in,1);
FermionField1f f1_fullgrid_dbl(out.Grid());
Replicate(f1_halfgrid, f1_fullgrid_dbl);
const Coordinate &dim_2f = in.Grid()->GlobalDimensions();
const Coordinate &dim_1f = out.Grid()->GlobalDimensions();
//We have to be careful for 5d fields; the s-direction is placed before the x,y,z,t and so we need to shift gpdir by 1
std::cout << "gpdir " << gpdir << std::endl;
gpdir+=1;
std::cout << "gpdir for 5D fields " << gpdir << std::endl;
std::cout << "dim_2f " << dim_2f << std::endl;
std::cout << "dim_1f " << dim_1f << std::endl;
assert(dim_1f[gpdir] == 2*dim_2f[gpdir]);
LatticeInteger xcoor_1f(out.Grid()); //5d lattice integer
LatticeCoordinate(xcoor_1f,gpdir);
int L = dim_2f[gpdir];
out = where(xcoor_1f < L, f0_fullgrid_dbl, f1_fullgrid_dbl);
}
//Both have the same field type
void copy2fTo1fGaugeField(LatticeGaugeField &out, const LatticeGaugeField &in, int gpdir){
LatticeGaugeField U_dbl(out.Grid());
Replicate(in, U_dbl);
LatticeGaugeField Uconj_dbl = conjugate( U_dbl );
const Coordinate &dim_2f = in.Grid()->GlobalDimensions();
LatticeInteger xcoor_1f(out.Grid());
LatticeCoordinate(xcoor_1f,gpdir);
int L = dim_2f[gpdir];
out = where(xcoor_1f < L, U_dbl, Uconj_dbl);
}
std::ostream & operator<<(std::ostream &os, const Coordinate &x){
os << "(";
for(int i=0;i<x.size();i++) os << x[i] << (i<x.size()-1 ? " " : "");
os << ")";
return os;
}
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
int Ls = 16;
Coordinate latt_2f = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd, vComplexD::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
int mu = 0; //Gparity direction
Coordinate latt_1f = latt_2f;
latt_1f[mu] *= 2;
GridCartesian * UGrid_1f = SpaceTimeGrid::makeFourDimGrid(latt_1f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_1f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_1f);
GridCartesian * FGrid_1f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_1f);
GridRedBlackCartesian * FrbGrid_1f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_1f);
GridCartesian * UGrid_2f = SpaceTimeGrid::makeFourDimGrid(latt_2f, simd_layout, mpi_layout);
GridRedBlackCartesian * UrbGrid_2f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_2f);
GridCartesian * FGrid_2f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_2f);
GridRedBlackCartesian * FrbGrid_2f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_2f);
std::cout << "SIMD layout " << simd_layout << std::endl;
std::cout << "MPI layout " << mpi_layout << std::endl;
std::cout << "2f dimensions " << latt_2f << std::endl;
std::cout << "1f dimensions " << latt_1f << std::endl;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4);
std::cout << "Generating hot 2f gauge configuration" << std::endl;
LatticeGaugeField Umu_2f(UGrid_2f);
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
std::cout << "Copying 2f->1f gauge field" << std::endl;
LatticeGaugeField Umu_1f(UGrid_1f);
copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu);
typedef GparityWilsonImplR FermionImplPolicy2f;
typedef GparityDomainWallFermionR FermionAction2f;
typedef typename FermionAction2f::FermionField FermionField2f;
typedef WilsonImplR FermionImplPolicy1f;
typedef DomainWallFermionR FermionAction1f;
typedef typename FermionAction1f::FermionField FermionField1f;
std::cout << "Generating eta 2f" << std::endl;
FermionField2f eta_2f(FGrid_2f);
gaussian(RNG5_2f, eta_2f);
RealD scale = std::sqrt(0.5);
eta_2f=eta_2f*scale;
std::cout << "Copying 2f->1f eta" << std::endl;
FermionField1f eta_1f(FGrid_1f);
copy2fTo1fFermionField(eta_1f, eta_2f, mu);
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.032;
Real pv_mass = 1.0;
RealD M5 = 1.8;
//Setup the Dirac operators
std::cout << "Initializing Dirac operators" << std::endl;
FermionAction2f::ImplParams Params_2f;
Params_2f.twists[mu] = 1;
Params_2f.twists[Nd-1] = 1; //APBC in time direction
//note 'Num' and 'Den' here refer to the determinant ratio, not the operator ratio in the pseudofermion action where the two are inverted
//to my mind the Pauli Villars and 'denominator' are synonymous but the Grid convention has this as the 'Numerator' operator in the RHMC implementation
FermionAction2f NumOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, light_mass,M5,Params_2f);
FermionAction2f DenOp_2f(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f, *UrbGrid_2f, pv_mass, M5,Params_2f);
FermionAction1f::ImplParams Params_1f;
Params_1f.boundary_phases[mu] = -1; //antiperiodic in doubled lattice in GP direction
Params_1f.boundary_phases[Nd-1] = -1;
FermionAction1f NumOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, light_mass,M5,Params_1f);
FermionAction1f DenOp_1f(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f, *UrbGrid_1f, pv_mass, M5,Params_1f);
//Test the replication routines by running a CG on eta
double StoppingCondition = 1e-10;
double MaxCGIterations = 30000;
ConjugateGradient<FermionField2f> CG_2f(StoppingCondition,MaxCGIterations);
ConjugateGradient<FermionField1f> CG_1f(StoppingCondition,MaxCGIterations);
NumOp_1f.ImportGauge(Umu_1f);
NumOp_2f.ImportGauge(Umu_2f);
FermionField1f test_1f(FGrid_1f);
FermionField2f test_2f(FGrid_2f);
MdagMLinearOperator<FermionAction1f, FermionField1f> Linop_1f(NumOp_1f);
MdagMLinearOperator<FermionAction2f, FermionField2f> Linop_2f(NumOp_2f);
CG_1f(Linop_1f, eta_1f, test_1f);
CG_2f(Linop_2f, eta_2f, test_2f);
RealD test_1f_norm = norm2(test_1f);
RealD test_2f_norm = norm2(test_2f);
std::cout << "Verification of replication routines: " << test_1f_norm << " " << test_2f_norm << " " << test_1f_norm - test_2f_norm << std::endl;
#if 1
typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy2f> Action2f;
typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy1f> Action1f;
RationalActionParams rational_params;
rational_params.inv_pow = 2;
rational_params.lo = 1e-5;
rational_params.hi = 32;
rational_params.md_degree = 16;
rational_params.action_degree = 16;
Action2f action_2f(DenOp_2f, NumOp_2f, rational_params);
Action1f action_1f(DenOp_1f, NumOp_1f, rational_params);
#else
typedef TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy2f> Action2f;
typedef TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy1f> Action1f;
Action2f action_2f(DenOp_2f, NumOp_2f, CG_2f, CG_2f);
Action1f action_1f(DenOp_1f, NumOp_1f, CG_1f, CG_1f);
#endif
std::cout << "Action refresh" << std::endl;
action_2f.refresh(Umu_2f, eta_2f);
action_1f.refresh(Umu_1f, eta_1f);
std::cout << "Action compute post heatbath" << std::endl;
RealD S_2f = action_2f.S(Umu_2f);
RealD S_1f = action_1f.S(Umu_1f);
std::cout << "Action comparison post heatbath" << std::endl;
std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl;
//Change the gauge field between refresh and action eval else the matrix and inverse matrices all cancel and we just get |eta|^2
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
copy2fTo1fGaugeField(Umu_1f, Umu_2f, mu);
//Now compute the action with the new gauge field
std::cout << "Action compute post gauge field update" << std::endl;
S_2f = action_2f.S(Umu_2f);
S_1f = action_1f.S(Umu_1f);
std::cout << "Action comparison post gauge field update" << std::endl;
std::cout << S_2f << " " << S_1f << " " << S_2f-S_1f << std::endl;
Grid_finalize();
} // main

4
tests/hmc/Test_hmc_GparityIwasakiGauge.cc
View File

@@ -58,7 +58,7 @@ int main(int argc, char **argv) {
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 5;
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@@ -79,7 +79,7 @@ int main(int argc, char **argv) {
// that have a complex construction
// standard
RealD beta = 2.6 ;
const int nu = 3;
const int nu = 1;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);

125
tests/solver/Test_eofa_inv.cc Normal file
View File

@@ -0,0 +1,125 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/solver/Test_eofa_inv.cc
Copyright (C) 2017
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: David Murphy <dmurphy@phys.columbia.edu>
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 */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
;
int main (int argc, char** argv)
{
Grid_init(&argc, &argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
const int Ls = 8;
GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
// Want a different conf at every run
// First create an instance of an engine.
std::random_device rnd_device;
// Specify the engine and distribution.
std::mt19937 mersenne_engine(rnd_device());
std::uniform_int_distribution<int> dist(1, 100);
auto gen = std::bind(dist, mersenne_engine);
std::vector<int> seeds4(4);
generate(begin(seeds4), end(seeds4), gen);
//std::vector<int> seeds4({1,2,3,5});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
LatticeFermion phi (FGrid); gaussian(RNG5, phi);
LatticeFermion Mphi (FGrid);
LatticeFermion MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD b = 2.5;
RealD c = 1.5;
RealD mf = 0.01;
RealD mb = 1.0;
RealD M5 = 1.8;
MobiusEOFAFermionR Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5, b, c);
OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-10, 12);
ConjugateGradient<LatticeFermion> CG(1.0e-10, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
//Random field
LatticeFermion eta(FGrid);
gaussian(RNG5,eta);
//Check left inverse
LatticeFermion Meta(FGrid);
Meofa.Meofa(U, eta, Meta);
LatticeFermion MinvMeta(FGrid);
Meofa.MeofaInv(U, Meta, MinvMeta);
LatticeFermion diff = MinvMeta - eta;
std::cout << GridLogMessage << "eta: " << norm2(eta) << " M*eta: " << norm2(Meta) << " M^{-1}*M*eta: " << norm2(MinvMeta) << " M^{-1}*M*eta - eta: " << norm2(diff) << " (expect 0)" << std::endl;
assert(norm2(diff) < 1e-8);
//Check right inverse
LatticeFermion MinvEta(FGrid);
Meofa.MeofaInv(U, eta, MinvEta);
LatticeFermion MMinvEta(FGrid);
Meofa.Meofa(U, MinvEta, MMinvEta);
diff = MMinvEta - eta;
std::cout << GridLogMessage << "eta: " << norm2(eta) << " M^{-1}*eta: " << norm2(MinvEta) << " M*M^{-1}*eta: " << norm2(MMinvEta) << " M*M^{-1}*eta - eta: " << norm2(diff) << " (expect 0)" << std::endl;
assert(norm2(diff) < 1e-8);
std::cout << GridLogMessage << "Done" << std::endl;
Grid_finalize();
}