portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-13 01:05:36 +00:00
Code Releases Activity

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

Browse Source
This commit is contained in:
Peter Boyle 2022-10-04 17:41:48 -04:00
commit c82b164f6b
35 changed files with 731 additions and 749 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h
View File

@ -154,6 +154,7 @@ public:
// dynamic sized arrays on stack; 2d is a pain with vector // dynamic sized arrays on stack; 2d is a pain with vector
RealD bs[nshift]; RealD bs[nshift];
RealD rsq[nshift]; RealD rsq[nshift];
RealD rsqf[nshift];
RealD z[nshift][2]; RealD z[nshift][2];
int converged[nshift]; int converged[nshift];
@ -190,6 +191,7 @@ public:
for(int s=0;s<nshift;s++){ for(int s=0;s<nshift;s++){
rsq[s] = cp * mresidual[s] * mresidual[s]; rsq[s] = cp * mresidual[s] * mresidual[s];
rsqf[s] =rsq[s];
std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift "<< s <<" target resid "<<rsq[s]<<std::endl; std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift "<< s <<" target resid "<<rsq[s]<<std::endl;
ps_d[s] = src_d; ps_d[s] = src_d;
} }
@ -198,7 +200,15 @@ public:
r_d = p_d; r_d = p_d;
//MdagM+m[0] //MdagM+m[0]
precisionChangeFast(p_f,p_d);
Linop_f.HermOpAndNorm(p_f,mmp_f,d,qq); // mmp = MdagM p d=real(dot(p, mmp)), qq=norm2(mmp)
precisionChangeFast(tmp_d,mmp_f);
Linop_d.HermOpAndNorm(p_d,mmp_d,d,qq); // mmp = MdagM p d=real(dot(p, mmp)), qq=norm2(mmp) Linop_d.HermOpAndNorm(p_d,mmp_d,d,qq); // mmp = MdagM p d=real(dot(p, mmp)), qq=norm2(mmp)
tmp_d = tmp_d - mmp_d;
std::cout << " Testing operators match "<<norm2(mmp_d)<<" f "<<norm2(mmp_f)<<" diff "<< norm2(tmp_d)<<std::endl;
// assert(norm2(tmp_d)< 1.0e-4);
axpy(mmp_d,mass[0],p_d,mmp_d); axpy(mmp_d,mass[0],p_d,mmp_d);
RealD rn = norm2(p_d); RealD rn = norm2(p_d);
d += rn*mass[0]; d += rn*mass[0];
@ -301,6 +311,7 @@ public:
//Perform reliable update if necessary; otherwise update residual from single-prec mmp //Perform reliable update if necessary; otherwise update residual from single-prec mmp
c = axpy_norm(r_d,b,mmp_d,r_d); c = axpy_norm(r_d,b,mmp_d,r_d);
AXPYTimer.Stop(); AXPYTimer.Stop();
if(k % ReliableUpdateFreq == 0){ if(k % ReliableUpdateFreq == 0){
@ -328,7 +339,7 @@ public:
RealD css = c * z[s][iz]* z[s][iz]; RealD css = c * z[s][iz]* z[s][iz];
if(css<rsq[s]){ if(css<rsqf[s]){
if ( ! converged[s] ) if ( ! converged[s] )
std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<<" Shift "<<s<<" has converged"<<std::endl; std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
converged[s]=1; converged[s]=1;

1
Grid/lattice/Lattice_reduction.h
View File

@ -91,7 +91,6 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites)
for(int i=0;i<nthread;i++){ for(int i=0;i<nthread;i++){
ssum = ssum+sumarray[i]; ssum = ssum+sumarray[i];
} }
return ssum; return ssum;
} }
/* /*

18
Grid/lattice/Lattice_transfer.h
View File

@ -1083,20 +1083,18 @@ vectorizeFromRevLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
template<class VobjOut, class VobjIn> template<class VobjOut, class VobjIn>
void precisionChangeFast(Lattice<VobjOut> &out, const Lattice<VobjIn> &in) void precisionChangeFast(Lattice<VobjOut> &out, const Lattice<VobjIn> &in)
{ {
typedef typename VobjOut::scalar_object SobjOut; typedef typename VobjOut::vector_type Vout;
typedef typename VobjIn::scalar_object SobjIn; typedef typename VobjIn::vector_type Vin;
const int N = sizeof(VobjOut)/sizeof(Vout);
conformable(out.Grid(),in.Grid()); conformable(out.Grid(),in.Grid());
out.Checkerboard() = in.Checkerboard(); out.Checkerboard() = in.Checkerboard();
int nsimd = out.Grid()->Nsimd(); int nsimd = out.Grid()->Nsimd();
autoView( out_v , out, AcceleratorWrite); autoView( out_v , out, AcceleratorWrite);
autoView( in_v , in, AcceleratorRead); autoView( in_v , in, AcceleratorRead);
accelerator_for(idx,out.Grid()->oSites(),nsimd,{ accelerator_for(idx,out.Grid()->oSites(),1,{
auto itmp = coalescedRead(in_v[idx]); Vout *vout = (Vout *)&out_v[idx];
auto otmp = coalescedRead(out_v[idx]); Vin *vin = (Vin *)&in_v[idx];
#ifdef GRID_SIMT precisionChange(vout,vin,N);
otmp=itmp;
#endif
coalescedWrite(out_v[idx],otmp);
}); });
} }
//Convert a Lattice from one precision to another //Convert a Lattice from one precision to another
@ -1116,7 +1114,7 @@ void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in)
int ndim = out.Grid()->Nd(); int ndim = out.Grid()->Nd();
int out_nsimd = out_grid->Nsimd(); int out_nsimd = out_grid->Nsimd();
int in_nsimd = in_grid->Nsimd();
std::vector<Coordinate > out_icoor(out_nsimd); std::vector<Coordinate > out_icoor(out_nsimd);
for(int lane=0; lane < out_nsimd; lane++){ for(int lane=0; lane < out_nsimd; lane++){

53
Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h
View File

@ -36,15 +36,19 @@ NAMESPACE_BEGIN(Grid);
// cf. GeneralEvenOddRational.h for details // cf. GeneralEvenOddRational.h for details
///////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class ImplD, class ImplF> template<class ImplD, class ImplF, class ImplD2>
class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> { class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> {
private: private:
typedef typename ImplD2::FermionField FermionFieldD2;
typedef typename ImplD::FermionField FermionFieldD; typedef typename ImplD::FermionField FermionFieldD;
typedef typename ImplF::FermionField FermionFieldF; typedef typename ImplF::FermionField FermionFieldF;
FermionOperator<ImplD> & NumOpD; FermionOperator<ImplD> & NumOpD;
FermionOperator<ImplD> & DenOpD; FermionOperator<ImplD> & DenOpD;
FermionOperator<ImplD2> & NumOpD2;
FermionOperator<ImplD2> & DenOpD2;
FermionOperator<ImplF> & NumOpF; FermionOperator<ImplF> & NumOpF;
FermionOperator<ImplF> & DenOpF; FermionOperator<ImplF> & DenOpF;
@ -53,37 +57,70 @@ NAMESPACE_BEGIN(Grid);
//Allow derived classes to override the multishift CG //Allow derived classes to override the multishift CG
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){ virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){
SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD); #if 0
SchurDifferentiableOperator<ImplD> schurOp(numerator ? NumOp : DenOp);
ConjugateGradientMultiShift<FermionFieldD> msCG(MaxIter, approx);
msCG(schurOp,in, out);
#else
SchurDifferentiableOperator<ImplD2> schurOpD2(numerator ? NumOpD2 : DenOpD2);
SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF); SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); ConjugateGradientMultiShiftMixedPrec<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
msCG(schurOpD, in, out); precisionChange(inD2,in);
std::cout << "msCG single solve "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
msCG(schurOpD2, inD2, outD2);
precisionChange(out,outD2);
#endif
} }
virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &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<ImplD2> schurOpD2(numerator ? NumOpD2 : DenOpD2);
SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF); SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
msCG(schurOpD, in, out_elems, out); FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
std::vector<FermionFieldD2> out_elemsD2(out_elems.size(),NumOpD2.FermionRedBlackGrid());
ConjugateGradientMultiShiftMixedPrec<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
precisionChange(inD2,in);
std::cout << "msCG in "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
msCG(schurOpD2, inD2, out_elemsD2, outD2);
precisionChange(out,outD2);
for(int i=0;i<out_elems.size();i++){
precisionChange(out_elems[i],out_elemsD2[i]);
}
} }
//Allow derived classes to override the gauge import //Allow derived classes to override the gauge import
virtual void ImportGauge(const typename ImplD::GaugeField &Ud){ virtual void ImportGauge(const typename ImplD::GaugeField &Ud){
typename ImplF::GaugeField Uf(NumOpF.GaugeGrid()); typename ImplF::GaugeField Uf(NumOpF.GaugeGrid());
typename ImplD2::GaugeField Ud2(NumOpD2.GaugeGrid());
precisionChange(Uf, Ud); precisionChange(Uf, Ud);
precisionChange(Ud2, Ud);
std::cout << "Importing "<<norm2(Ud)<<" "<< norm2(Uf)<<" " << norm2(Ud2)<<std::endl;
NumOpD.ImportGauge(Ud); NumOpD.ImportGauge(Ud);
DenOpD.ImportGauge(Ud); DenOpD.ImportGauge(Ud);
NumOpF.ImportGauge(Uf); NumOpF.ImportGauge(Uf);
DenOpF.ImportGauge(Uf); DenOpF.ImportGauge(Uf);
NumOpD2.ImportGauge(Ud2);
DenOpD2.ImportGauge(Ud2);
} }
public: public:
GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD> &_NumOpD, FermionOperator<ImplD> &_DenOpD, GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD> &_NumOpD, FermionOperator<ImplD> &_DenOpD,
FermionOperator<ImplF> &_NumOpF, FermionOperator<ImplF> &_DenOpF, FermionOperator<ImplF> &_NumOpF, FermionOperator<ImplF> &_DenOpF,
FermionOperator<ImplD2> &_NumOpD2, FermionOperator<ImplD2> &_DenOpD2,
const RationalActionParams & p, Integer _ReliableUpdateFreq const RationalActionParams & p, Integer _ReliableUpdateFreq
) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p), ) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p),
ReliableUpdateFreq(_ReliableUpdateFreq), NumOpD(_NumOpD), DenOpD(_DenOpD), NumOpF(_NumOpF), DenOpF(_DenOpF){} ReliableUpdateFreq(_ReliableUpdateFreq),
NumOpD(_NumOpD), DenOpD(_DenOpD),
NumOpF(_NumOpF), DenOpF(_DenOpF),
NumOpD2(_NumOpD2), DenOpD2(_DenOpD2)
{}
virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";} virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";}
}; };

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

@ -67,8 +67,9 @@ NAMESPACE_BEGIN(Grid);
virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";}
}; };
template<class Impl,class ImplF> template<class Impl,class ImplF,class ImplD2>
class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF> { class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction
: public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF,ImplD2> {
public: public:
typedef OneFlavourRationalParams Params; typedef OneFlavourRationalParams Params;
private: private:
@ -90,9 +91,11 @@ NAMESPACE_BEGIN(Grid);
FermionOperator<Impl> &_DenOp, FermionOperator<Impl> &_DenOp,
FermionOperator<ImplF> &_NumOpF, FermionOperator<ImplF> &_NumOpF,
FermionOperator<ImplF> &_DenOpF, FermionOperator<ImplF> &_DenOpF,
FermionOperator<ImplD2> &_NumOpD2,
FermionOperator<ImplD2> &_DenOpD2,
const Params & p, Integer ReliableUpdateFreq const Params & p, Integer ReliableUpdateFreq
) : ) :
GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF>(_NumOp, _DenOp,_NumOpF, _DenOpF, transcribe(p),ReliableUpdateFreq){} GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF,ImplD2>(_NumOp, _DenOp,_NumOpF, _DenOpF,_NumOpD2, _DenOpD2, transcribe(p),ReliableUpdateFreq){}
virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";}
}; };

17
Grid/qcd/hmc/HMC.h
View File

@ -53,6 +53,7 @@ struct HMCparameters: Serializable {
Integer, Trajectories, /* @brief Number of sweeps in this run */ Integer, Trajectories, /* @brief Number of sweeps in this run */
bool, MetropolisTest, bool, MetropolisTest,
Integer, NoMetropolisUntil, Integer, NoMetropolisUntil,
bool, PerformRandomShift, /* @brief Randomly shift the gauge configuration at the start of a trajectory */
std::string, StartingType, std::string, StartingType,
IntegratorParameters, MD) IntegratorParameters, MD)
@ -63,6 +64,7 @@ struct HMCparameters: Serializable {
StartTrajectory = 0; StartTrajectory = 0;
Trajectories = 10; Trajectories = 10;
StartingType = "HotStart"; StartingType = "HotStart";
PerformRandomShift = true;
///////////////////////////////// /////////////////////////////////
} }
@ -83,6 +85,7 @@ struct HMCparameters: Serializable {
std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n"; std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n"; std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n"; std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
std::cout << GridLogMessage << "[HMC parameters] Doing random shift : " << std::boolalpha << PerformRandomShift << "\n";
std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n"; std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
MD.print_parameters(); MD.print_parameters();
} }
@ -95,6 +98,7 @@ private:
const HMCparameters Params; const HMCparameters Params;
typedef typename IntegratorType::Field Field; typedef typename IntegratorType::Field Field;
typedef typename IntegratorType::FieldImplementation FieldImplementation;
typedef std::vector< HmcObservable<Field> * > ObsListType; typedef std::vector< HmcObservable<Field> * > ObsListType;
//pass these from the resource manager //pass these from the resource manager
@ -138,11 +142,16 @@ private:
GridBase *Grid = U.Grid(); GridBase *Grid = U.Grid();
if(Params.PerformRandomShift){
////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////
// Mainly for DDHMC perform a random translation of U modulo volume // Mainly for DDHMC perform a random translation of U modulo volume
////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << GridLogMessage << "--------------------------------------------------\n"; std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Random shifting gauge field by ["; std::cout << GridLogMessage << "Random shifting gauge field by [";
std::vector<typename FieldImplementation::GaugeLinkField> Umu(Grid->Nd(), U.Grid());
for(int mu=0;mu<Grid->Nd();mu++) Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
for(int d=0;d<Grid->Nd();d++) { for(int d=0;d<Grid->Nd();d++) {
int L = Grid->GlobalDimensions()[d]; int L = Grid->GlobalDimensions()[d];
@ -155,9 +164,15 @@ private:
if(d<Grid->Nd()-1) std::cout <<","; if(d<Grid->Nd()-1) std::cout <<",";
else std::cout <<"]\n"; else std::cout <<"]\n";
U = Cshift(U,d,shift); //shift all fields together in a way that respects the gauge BCs
for(int mu=0; mu < Grid->Nd(); mu++)
Umu[mu] = FieldImplementation::CshiftLink(Umu[mu],d,shift);
} }
for(int mu=0;mu<Grid->Nd();mu++) PokeIndex<LorentzIndex>(U,Umu[mu],mu);
std::cout << GridLogMessage << "--------------------------------------------------\n"; std::cout << GridLogMessage << "--------------------------------------------------\n";
}
TheIntegrator.reset_timer(); TheIntegrator.reset_timer();

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

@ -63,10 +63,10 @@ public:
}; };
/*! @brief Class for Molecular Dynamics management */ /*! @brief Class for Molecular Dynamics management */
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy> template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy>
class Integrator { class Integrator {
protected: protected:
typedef FieldImplementation_ FieldImplementation;
typedef typename FieldImplementation::Field MomentaField; //for readability typedef typename FieldImplementation::Field MomentaField; //for readability
typedef typename FieldImplementation::Field Field; typedef typename FieldImplementation::Field Field;

15
Grid/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid);
* P 1/2 P 1/2 * P 1/2 P 1/2
*/ */
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> > template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class LeapFrog : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> class LeapFrog : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{ {
public: public:
typedef FieldImplementation_ FieldImplementation;
typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy> Algorithm; typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy> Algorithm;
INHERIT_FIELD_TYPES(FieldImplementation); INHERIT_FIELD_TYPES(FieldImplementation);
@ -135,13 +136,14 @@ public:
} }
}; };
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> > template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class MinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> class MinimumNorm2 : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{ {
private: private:
const RealD lambda = 0.1931833275037836; const RealD lambda = 0.1931833275037836;
public: public:
typedef FieldImplementation_ FieldImplementation;
INHERIT_FIELD_TYPES(FieldImplementation); INHERIT_FIELD_TYPES(FieldImplementation);
MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm) MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
@ -192,8 +194,8 @@ public:
} }
}; };
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> > template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy = Representations<FundamentalRepresentation> >
class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy> class ForceGradient : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy>
{ {
private: private:
const RealD lambda = 1.0 / 6.0; const RealD lambda = 1.0 / 6.0;
@ -202,6 +204,7 @@ private:
const RealD theta = 0.0; const RealD theta = 0.0;
public: public:
typedef FieldImplementation_ FieldImplementation;
INHERIT_FIELD_TYPES(FieldImplementation); INHERIT_FIELD_TYPES(FieldImplementation);
// Looks like dH scales as dt^4. tested wilson/wilson 2 level. // Looks like dH scales as dt^4. tested wilson/wilson 2 level.

15
Grid/qcd/observables/topological_charge.h
View File

@ -31,15 +31,16 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
struct TopologySmearingParameters : Serializable { struct TopologySmearingParameters : Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters, GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters,
int, steps,
float, step_size,
int, meas_interval, int, meas_interval,
float, maxTau); float, init_step_size,
float, maxTau,
float, tolerance);
TopologySmearingParameters(int s = 0, float ss = 0.0f, int mi = 0, float mT = 0.0f): TopologySmearingParameters(float ss = 0.0f, int mi = 0, float mT = 0.0f, float tol = 1e-4):
steps(s), step_size(ss), meas_interval(mi), maxTau(mT){} init_step_size(ss), meas_interval(mi), maxTau(mT), tolerance(tol){}
template < class ReaderClass > template < class ReaderClass >
TopologySmearingParameters(Reader<ReaderClass>& Reader){ TopologySmearingParameters(Reader<ReaderClass>& Reader){
@ -97,8 +98,8 @@ public:
if (Pars.do_smearing){ if (Pars.do_smearing){
// using wilson flow by default here // using wilson flow by default here
WilsonFlow<PeriodicGimplR> WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval); WilsonFlowAdaptive<PeriodicGimplR> WF(Pars.Smearing.init_step_size, Pars.Smearing.maxTau, Pars.Smearing.tolerance, Pars.Smearing.meas_interval);
WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau); WF.smear(Usmear, U);
Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear); Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear);
std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1) std::cout << GridLogMessage << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
<< "T0 : [ " << traj << " ] "<< T0 << std::endl; << "T0 : [ " << traj << " ] "<< T0 << std::endl;

288
Grid/qcd/smearing/WilsonFlow.h
View File

@ -33,27 +33,25 @@ directory
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
template <class Gimpl> template <class Gimpl>
class WilsonFlow: public Smear<Gimpl>{ class WilsonFlowBase: public Smear<Gimpl>{
public: public:
//Store generic measurements to take during smearing process using std::function //Store generic measurements to take during smearing process using std::function
typedef std::function<void(int, RealD, const typename Gimpl::GaugeField &)> FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field typedef std::function<void(int, RealD, const typename Gimpl::GaugeField &)> FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field
private: protected:
unsigned int Nstep;
RealD epsilon; //for regular smearing this is the time step, for adaptive it is the initial time step
std::vector< std::pair<int, FunctionType> > functions; //The int maps to the measurement frequency std::vector< std::pair<int, FunctionType> > functions; //The int maps to the measurement frequency
mutable WilsonGaugeAction<Gimpl> SG; mutable WilsonGaugeAction<Gimpl> SG;
//Evolve the gauge field by 1 step and update tau
void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const;
//Evolve the gauge field by 1 step and update tau and the current time step eps
void evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps, RealD maxTau) const;
public: public:
INHERIT_GIMPL_TYPES(Gimpl) INHERIT_GIMPL_TYPES(Gimpl)
explicit WilsonFlowBase(unsigned int meas_interval =1):
SG(WilsonGaugeAction<Gimpl>(3.0)) {
// WilsonGaugeAction with beta 3.0
setDefaultMeasurements(meas_interval);
}
void resetActions(){ functions.clear(); } void resetActions(){ functions.clear(); }
void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); } void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); }
@ -64,34 +62,11 @@ public:
//and output to stdout //and output to stdout
void setDefaultMeasurements(int topq_meas_interval = 1); void setDefaultMeasurements(int topq_meas_interval = 1);
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1): void derivative(GaugeField&, const GaugeField&, const GaugeField&) const override{
Nstep(Nstep),
epsilon(epsilon),
SG(WilsonGaugeAction<Gimpl>(3.0)) {
// WilsonGaugeAction with beta 3.0
assert(epsilon > 0.0);
LogMessage();
setDefaultMeasurements(interval);
}
void LogMessage() {
std::cout << GridLogMessage
<< "[WilsonFlow] Nstep : " << Nstep << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] epsilon : " << epsilon << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
}
virtual void smear(GaugeField&, const GaugeField&) const;
virtual void derivative(GaugeField&, const GaugeField&, const GaugeField&) const {
assert(0); assert(0);
// undefined for WilsonFlow // undefined for WilsonFlow
} }
void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau) const;
//Compute t^2 <E(t)> for time t from the plaquette //Compute t^2 <E(t)> for time t from the plaquette
static RealD energyDensityPlaquette(const RealD t, const GaugeField& U); static RealD energyDensityPlaquette(const RealD t, const GaugeField& U);
@ -115,82 +90,63 @@ public:
std::vector<RealD> flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1); std::vector<RealD> flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1);
}; };
//Basic iterative Wilson flow
template <class Gimpl>
class WilsonFlow: public WilsonFlowBase<Gimpl>{
private:
int Nstep; //number of steps
RealD epsilon; //step size
//Evolve the gauge field by 1 step of size eps and update tau
void evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const;
public:
INHERIT_GIMPL_TYPES(Gimpl)
//Integrate the Wilson flow for Nstep steps of size epsilon
WilsonFlow(const RealD epsilon, const int Nstep, unsigned int meas_interval = 1): WilsonFlowBase<Gimpl>(meas_interval), Nstep(Nstep), epsilon(epsilon){}
void smear(GaugeField& out, const GaugeField& in) const override;
};
//Wilson flow with adaptive step size
template <class Gimpl>
class WilsonFlowAdaptive: public WilsonFlowBase<Gimpl>{
private:
RealD init_epsilon; //initial step size
RealD maxTau; //integrate to t=maxTau
RealD tolerance; //integration error tolerance
//Evolve the gauge field by 1 step and update tau and the current time step eps
//
//If the step size eps is too large that a significant integration error results,
//the gauge field (U) and tau will not be updated and the function will return 0; eps will be adjusted to a smaller
//value for the next iteration.
//
//For a successful integration step the function will return 1
int evolve_step_adaptive(typename Gimpl::GaugeField&U, RealD &tau, RealD &eps) const;
public:
INHERIT_GIMPL_TYPES(Gimpl)
WilsonFlowAdaptive(const RealD init_epsilon, const RealD maxTau, const RealD tolerance, unsigned int meas_interval = 1):
WilsonFlowBase<Gimpl>(meas_interval), init_epsilon(init_epsilon), maxTau(maxTau), tolerance(tolerance){}
void smear(GaugeField& out, const GaugeField& in) const override;
};
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// Implementations // Implementations
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{ RealD WilsonFlowBase<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
GaugeField Z(U.Grid());
GaugeField tmp(U.Grid());
SG.deriv(U, Z);
Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0;
SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1
Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0;
SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
tau += epsilon;
}
template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps, RealD maxTau) const{
if (maxTau - tau < eps){
eps = maxTau-tau;
}
//std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
GaugeField Z(U.Grid());
GaugeField Zprime(U.Grid());
GaugeField tmp(U.Grid()), Uprime(U.Grid());
Uprime = U;
SG.deriv(U, Z);
Zprime = -Z;
Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*eps); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0;
SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Zprime += 2.0*tmp;
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1
Gimpl::update_field(Z, U, -2.0*eps); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0;
SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
Gimpl::update_field(Z, U, -2.0*eps); // V(t+e) = exp(ep*Z)*W2
// Ramos
Gimpl::update_field(Zprime, Uprime, -2.0*eps); // V'(t+e) = exp(ep*Z')*W0
// Compute distance as norm^2 of the difference
GaugeField diffU = U - Uprime;
RealD diff = norm2(diffU);
// adjust integration step
tau += eps;
//std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
eps = eps*0.95*std::pow(1e-4/diff,1./3.);
//std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
}
template <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityPlaquette(const RealD t, const GaugeField& U){
static WilsonGaugeAction<Gimpl> SG(3.0); static WilsonGaugeAction<Gimpl> SG(3.0);
return 2.0 * t * t * SG.S(U)/U.Grid()->gSites(); return 2.0 * t * t * SG.S(U)/U.Grid()->gSites();
} }
//Compute t^2 <E(t)> for time from the 1x1 cloverleaf form //Compute t^2 <E(t)> for time from the 1x1 cloverleaf form
template <class Gimpl> template <class Gimpl>
RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){ RealD WilsonFlowBase<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField& U){
typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz; typedef typename Gimpl::GaugeField GaugeLorentz;
@ -215,7 +171,7 @@ RealD WilsonFlow<Gimpl>::energyDensityCloverleaf(const RealD t, const GaugeField
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){
std::vector<RealD> out; std::vector<RealD> out;
resetActions(); resetActions();
addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){
@ -227,13 +183,13 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(GaugeFie
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){
GaugeField V(U); GaugeField V(U);
return flowMeasureEnergyDensityPlaquette(V,U, measure_interval); return flowMeasureEnergyDensityPlaquette(V,U, measure_interval);
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){
std::vector<RealD> out; std::vector<RealD> out;
resetActions(); resetActions();
addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){
@ -245,16 +201,52 @@ std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(GaugeFi
} }
template <class Gimpl> template <class Gimpl>
std::vector<RealD> WilsonFlow<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){
GaugeField V(U); GaugeField V(U);
return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval); return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval);
} }
template <class Gimpl>
void WilsonFlowBase<Gimpl>::setDefaultMeasurements(int topq_meas_interval){
addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << t << " " << energyDensityPlaquette(t,U) << std::endl;
});
addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " << step << " " << WilsonLoops<Gimpl>::TopologicalCharge(U) << std::endl;
});
}
//#define WF_TIMING
template <class Gimpl>
void WilsonFlow<Gimpl>::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{
GaugeField Z(U.Grid());
GaugeField tmp(U.Grid());
this->SG.deriv(U, Z);
Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*epsilon); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0;
this->SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1
Gimpl::update_field(Z, U, -2.0*epsilon); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0;
this->SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
tau += epsilon;
}
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
std::cout << GridLogMessage
<< "[WilsonFlow] Nstep : " << Nstep << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] epsilon : " << epsilon << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
out = in; out = in;
RealD taus = 0.; RealD taus = 0.;
for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement
@ -266,36 +258,92 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
std::cout << "Time to evolve " << diff.count() << " s\n"; std::cout << "Time to evolve " << diff.count() << " s\n";
#endif #endif
//Perform measurements //Perform measurements
for(auto const &meas : functions) for(auto const &meas : this->functions)
if( step % meas.first == 0 ) meas.second(step,taus,out); if( step % meas.first == 0 ) meas.second(step,taus,out);
} }
} }
template <class Gimpl> template <class Gimpl>
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau) const{ int WilsonFlowAdaptive<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, RealD &tau, RealD &eps) const{
if (maxTau - tau < eps){
eps = maxTau-tau;
}
//std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
GaugeField Z(U.Grid());
GaugeField Zprime(U.Grid());
GaugeField tmp(U.Grid()), Uprime(U.Grid()), Usave(U.Grid());
Uprime = U;
Usave = U;
this->SG.deriv(U, Z);
Zprime = -Z;
Z *= 0.25; // Z0 = 1/4 * F(U)
Gimpl::update_field(Z, U, -2.0*eps); // U = W1 = exp(ep*Z0)*W0
Z *= -17.0/8.0;
this->SG.deriv(U, tmp); Z += tmp; // -17/32*Z0 +Z1
Zprime += 2.0*tmp;
Z *= 8.0/9.0; // Z = -17/36*Z0 +8/9*Z1
Gimpl::update_field(Z, U, -2.0*eps); // U_= W2 = exp(ep*Z)*W1
Z *= -4.0/3.0;
this->SG.deriv(U, tmp); Z += tmp; // 4/3*(17/36*Z0 -8/9*Z1) +Z2
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
Gimpl::update_field(Z, U, -2.0*eps); // V(t+e) = exp(ep*Z)*W2
// Ramos arXiv:1301.4388
Gimpl::update_field(Zprime, Uprime, -2.0*eps); // V'(t+e) = exp(ep*Z')*W0
// Compute distance using Ramos' definition
GaugeField diffU = U - Uprime;
RealD max_dist = 0;
for(int mu=0;mu<Nd;mu++){
typename Gimpl::GaugeLinkField diffU_mu = PeekIndex<LorentzIndex>(diffU, mu);
RealD dist_mu = sqrt( maxLocalNorm2(diffU_mu) ) /Nc/Nc; //maximize over sites
max_dist = std::max(max_dist, dist_mu); //maximize over mu
}
int ret;
if(max_dist < tolerance) {
tau += eps;
ret = 1;
} else {
U = Usave;
ret = 0;
}
eps = eps*0.95*std::pow(tolerance/max_dist,1./3.);
std::cout << GridLogMessage << "Adaptive smearing : Distance: "<< max_dist <<" Step successful: " << ret << " New epsilon: " << eps << std::endl;
return ret;
}
template <class Gimpl>
void WilsonFlowAdaptive<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
std::cout << GridLogMessage
<< "[WilsonFlow] initial epsilon : " << init_epsilon << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] full trajectory : " << maxTau << std::endl;
std::cout << GridLogMessage
<< "[WilsonFlow] tolerance : " << tolerance << std::endl;
out = in; out = in;
RealD taus = 0.; RealD taus = 0.;
RealD eps = epsilon; RealD eps = init_epsilon;
unsigned int step = 0; unsigned int step = 0;
do{ do{
step++; int step_success = evolve_step_adaptive(out, taus, eps);
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; step += step_success; //step will not be incremented if the integration step fails
evolve_step_adaptive(out, taus, eps, maxTau);
//Perform measurements //Perform measurements
for(auto const &meas : functions) if(step_success)
for(auto const &meas : this->functions)
if( step % meas.first == 0 ) meas.second(step,taus,out); if( step % meas.first == 0 ) meas.second(step,taus,out);
} while (taus < maxTau); } while (taus < maxTau);
} }
template <class Gimpl>
void WilsonFlow<Gimpl>::setDefaultMeasurements(int topq_meas_interval){
addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << t << " " << energyDensityPlaquette(t,U) << std::endl;
});
addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " << step << " " << WilsonLoops<Gimpl>::TopologicalCharge(U) << std::endl;
});
}
NAMESPACE_END(Grid); NAMESPACE_END(Grid);

30
Grid/qcd/utils/CovariantCshift.h
View File

@ -227,26 +227,38 @@ namespace ConjugateBC {
//shift = -1 //shift = -1
//Out(x) = U_\mu(x-mu) | x_\mu != 0 //Out(x) = U_\mu(x-mu) | x_\mu != 0
// = U*_\mu(L-1) | x_\mu == 0 // = U*_\mu(L-1) | x_\mu == 0
//shift = 2
//Out(x) = U_\mu(x+2\hat\mu) | x_\mu < L-2
// = U*_\mu(1) | x_\mu == L-1
// = U*_\mu(0) | x_\mu == L-2
//shift = -2
//Out(x) = U_\mu(x-2mu) | x_\mu > 1
// = U*_\mu(L-2) | x_\mu == 0
// = U*_\mu(L-1) | x_\mu == 1
//etc
template<class gauge> Lattice<gauge> template<class gauge> Lattice<gauge>
CshiftLink(const Lattice<gauge> &Link, int mu, int shift) CshiftLink(const Lattice<gauge> &Link, int mu, int shift)
{ {
GridBase *grid = Link.Grid(); GridBase *grid = Link.Grid();
int Lmu = grid->GlobalDimensions()[mu] - 1; int Lmu = grid->GlobalDimensions()[mu];
assert(abs(shift) < Lmu && "Invalid shift value");
Lattice<iScalar<vInteger>> coor(grid); Lattice<iScalar<vInteger>> coor(grid);
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
Lattice<gauge> tmp(grid); Lattice<gauge> tmp(grid);
if(shift == 1){ if(shift > 0){
tmp = Cshift(Link, mu, 1); tmp = Cshift(Link, mu, shift);
tmp = where(coor == Lmu, conjugate(tmp), tmp); tmp = where(coor >= Lmu-shift, conjugate(tmp), tmp);
return tmp; return tmp;
}else if(shift == -1){ }else if(shift < 0){
tmp = Link; tmp = Link;
tmp = where(coor == Lmu, conjugate(tmp), tmp); tmp = where(coor >= Lmu+shift, conjugate(tmp), tmp);
return Cshift(tmp, mu, -1); return Cshift(tmp, mu, shift);
}else assert(0 && "Invalid shift value"); }
return tmp; //shuts up the compiler fussing about the return type
//shift == 0
return Link;
} }
} }

4
Grid/qcd/utils/GaugeFix.h
View File

@ -72,12 +72,12 @@ public:
//Fix the gauge field Umu //Fix the gauge field Umu
//0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf //0 < alpha < 1 is related to the step size, cf https://arxiv.org/pdf/1405.5812.pdf
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
GridBase *grid = Umu.Grid(); GridBase *grid = Umu.Grid();
GaugeMat xform(grid); GaugeMat xform(grid);
SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge); SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge);
} }
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
//Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform //Fix the gauge field Umu and also return the gauge transformation from the original gauge field, xform
GridBase *grid = Umu.Grid(); GridBase *grid = Umu.Grid();

2
Grid/simd/Grid_sse4.h
View File

@ -35,7 +35,7 @@ Author: neo <cossu@post.kek.jp>
*/ */
// Time-stamp: <2015-06-16 23:27:54 neo> // Time-stamp: <2015-06-16 23:27:54 neo>
//---------------------------------------------------------------------- //----------------------------------------------------------------------
#include <immintrin.h>
#include <pmmintrin.h> #include <pmmintrin.h>
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);

7
Grid/simd/Simd.h
View File

@ -258,6 +258,13 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexD &o){
stream<<">"; stream<<">";
return stream; return stream;
} }
inline std::ostream& operator<< (std::ostream& stream, const vComplexD2 &o){
stream<<"<";
stream<<o.v[0];
stream<<o.v[1];
stream<<">";
return stream;
}
inline std::ostream& operator<< (std::ostream& stream, const vRealF &o){ inline std::ostream& operator<< (std::ostream& stream, const vRealF &o){
int nn=vRealF::Nsimd(); int nn=vRealF::Nsimd();

2
HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc
View File

@ -492,6 +492,7 @@ NAMESPACE_END(Grid);
int main(int argc, char **argv) { int main(int argc, char **argv) {
#if 0
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run // here make a routine to print all the relevant information on the run
@ -915,4 +916,5 @@ int main(int argc, char **argv) {
std::cout << GridLogMessage << " Done" << std::endl; std::cout << GridLogMessage << " Done" << std::endl;
Grid_finalize(); Grid_finalize();
return 0; return 0;
#endif
} // main } // main

4
HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc
View File

@ -491,6 +491,7 @@ NAMESPACE_END(Grid);
int main(int argc, char **argv) { int main(int argc, char **argv) {
#if 0
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run // here make a routine to print all the relevant information on the run
@ -549,7 +550,7 @@ int main(int argc, char **argv) {
typedef typename FermionActionF::Impl_t FermionImplPolicyF; typedef typename FermionActionF::Impl_t FermionImplPolicyF;
typedef typename FermionActionF::FermionField FermionFieldF; typedef typename FermionActionF::FermionField FermionFieldF;
typedef GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicyD,FermionImplPolicyF> MixedPrecRHMC; typedef GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicyD,FermionImplPolicyF,FermionImplPolicyD2> MixedPrecRHMC;
typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicyD> DoublePrecRHMC; typedef GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicyD> DoublePrecRHMC;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@ -870,4 +871,5 @@ int main(int argc, char **argv) {
std::cout << GridLogMessage << " Done" << std::endl; std::cout << GridLogMessage << " Done" << std::endl;
Grid_finalize(); Grid_finalize();
return 0; return 0;
#endif
} // main } // main

6
HMC/Mobius2p1f_DD_RHMC_96I.cc
View File

@ -286,8 +286,8 @@ int main(int argc, char **argv) {
// ii) Break low bound, how rapidly? // ii) Break low bound, how rapidly?
// iii) Run lanczos // iii) Run lanczos
// iv) Have CG return spectral range estimate // iv) Have CG return spectral range estimate
FermionField vec(StrangeOp.FermionDedBlackGrid()); FermionField vec(StrangeOp.FermionRedBlackGrid());
FermionField res(StrangeOp.FermionDedBlackGrid()); FermionField res(StrangeOp.FermionRedBlackGrid());
vec = 1; // Fill with any old junk vec = 1; // Fill with any old junk
std::cout << "Bounds check on strange operator mass "<< StrangeOp.Mass()<<std::endl; std::cout << "Bounds check on strange operator mass "<< StrangeOp.Mass()<<std::endl;
@ -327,7 +327,7 @@ int main(int argc, char **argv) {
auto grid4= GridPtr; auto grid4= GridPtr;
auto rbgrid4= GridRBPtr; auto rbgrid4= GridRBPtr;
auto rbgrid = StrangeOp.FermionDedBlackGrid(); auto rbgrid = StrangeOp.FermionRedBlackGrid();
auto grid = StrangeOp.FermionGrid(); auto grid = StrangeOp.FermionGrid();
if(1){ if(1){
const int Nstop = 5; const int Nstop = 5;

492
HMC/Mobius2p1f_DD_RHMC_96I_mixedmshift.cc
View File

@ -1,492 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_EODWFRatio.cc
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@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>
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.)
{
/* Debugging instances of objects; references are stored
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl;
*/
};
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);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
RealD delta=1.e-4;
std::cout << GridLogMessage << "Calling reliable update Conjugate Gradient" <<std::endl;
ConjugateGradientReliableUpdate<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations*MaxOuterIterations,delta,SinglePrecGrid5,LinOpF,LinOpD);
MPCG(src,psi);
}
};
NAMESPACE_END(Grid);
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// Typedefs to simplify notation
typedef WilsonImplR FermionImplPolicy;
typedef WilsonImplF FermionImplPolicyF;
typedef MobiusFermionD FermionAction;
typedef MobiusFermionF FermionActionF;
typedef typename FermionAction::FermionField FermionField;
typedef typename FermionActionF::FermionField FermionFieldF;
typedef Grid::XmlReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2");
MD.MDsteps = 6;
MD.trajL = 1.0;
HMCparameters HMCparams;
HMCparams.StartTrajectory = 1077;
HMCparams.Trajectories = 1;
HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
// HMCparams.StartingType =std::string("ColdStart");
HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams);
// Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_DDHMC_lat";
CPparams.rng_prefix = "ckpoint_DDHMC_rng";
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
const int Ls = 12;
RealD M5 = 1.8;
RealD b = 1.5;
RealD c = 0.5;
Real beta = 2.31;
// Real light_mass = 5.4e-4;
Real light_mass = 7.8e-4;
Real strange_mass = 0.02132;
Real pv_mass = 1.0;
// std::vector<Real> hasenbusch({ light_mass, 3.8e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass });
std::vector<Real> hasenbusch({ light_mass, 5e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass });
// FIXME:
// Same in MC and MD
// Need to mix precision too
OneFlavourRationalParams SFRp; // Strange
SFRp.lo = 4.0e-3;
SFRp.hi = 90.0;
SFRp.MaxIter = 60000;
SFRp.tolerance= 1.0e-8;
SFRp.mdtolerance= 1.0e-6;
SFRp.degree = 12;
SFRp.precision= 50;
SFRp.BoundsCheckFreq=0;
OneFlavourRationalParams OFRp; // Up/down
OFRp.lo = 2.0e-5;
OFRp.hi = 90.0;
OFRp.MaxIter = 60000;
OFRp.tolerance= 1.0e-8;
OFRp.mdtolerance= 1.0e-6;
// OFRp.degree = 20; converges
// OFRp.degree = 16;
OFRp.degree = 12;
OFRp.precision= 80;
OFRp.BoundsCheckFreq=0;
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
typedef SchurDiagMooeeOperator<FermionAction ,FermionField > LinearOperatorD;
typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusFermionD,MobiusFermionF,LinearOperatorD,LinearOperatorF> MxPCG;
////////////////////////////////////////////////////////////////
// Domain decomposed
////////////////////////////////////////////////////////////////
Coordinate latt4 = GridPtr->GlobalDimensions();
Coordinate mpi = GridPtr->ProcessorGrid();
Coordinate shm;
GlobalSharedMemory::GetShmDims(mpi,shm);
Coordinate CommDim(Nd);
for(int d=0;d<Nd;d++) CommDim[d]= (mpi[d]/shm[d])>1 ? 1 : 0;
Coordinate NonDirichlet(Nd+1,0);
Coordinate Dirichlet(Nd+1,0);
Dirichlet[1] = CommDim[0]*latt4[0]/mpi[0] * shm[0];
Dirichlet[2] = CommDim[1]*latt4[1]/mpi[1] * shm[1];
Dirichlet[3] = CommDim[2]*latt4[2]/mpi[2] * shm[2];
Dirichlet[4] = CommDim[3]*latt4[3]/mpi[3] * shm[3];
Coordinate Block4(Nd);
Block4[0] = Dirichlet[1];
Block4[1] = Dirichlet[2];
Block4[2] = Dirichlet[3];
Block4[3] = Dirichlet[4];
int Width=3;
TheHMC.Resources.SetMomentumFilter(new DDHMCFilter<WilsonImplR::Field>(Block4,Width));
//////////////////////////
// Fermion Grids
//////////////////////////
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd());
auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt4,simdF,mpi);
auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF);
auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF);
auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF);
IwasakiGaugeActionR GaugeAction(beta);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
LatticeGaugeFieldF UF(GridPtrF);
std::cout << GridLogMessage << " Running the HMC "<< std::endl;
TheHMC.ReadCommandLine(argc,argv); // params on CML or from param file
TheHMC.initializeGaugeFieldAndRNGs(U);
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
FermionAction::ImplParams ParamsDir(boundary);
FermionActionF::ImplParams ParamsF(boundary);
FermionActionF::ImplParams ParamsDirF(boundary);
Params.dirichlet=NonDirichlet;
ParamsF.dirichlet=NonDirichlet;
ParamsDir.dirichlet=Dirichlet;
ParamsDirF.dirichlet=Dirichlet;
// double StoppingCondition = 1e-14;
// double MDStoppingCondition = 1e-9;
double StoppingCondition = 1e-10;
double MDStoppingCondition = 1e-7;
double MDStoppingConditionLoose = 1e-6;
double MaxCGIterations = 300000;
ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations);
ConjugateGradient<FermionField> MDCG(MDStoppingCondition,MaxCGIterations);
////////////////////////////////////
// Collect actions
////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(4);
ActionLevel<HMCWrapper::Field> Level3(8);
////////////////////////////////////
// Strange action
////////////////////////////////////
FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params);
FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params);
FermionAction StrangeOpDir (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, ParamsDir);
FermionAction StrangePauliVillarsOpDir(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, ParamsDir);
FermionActionF StrangeOpF (UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,strange_mass,M5,b,c, ParamsF);
FermionActionF StrangePauliVillarsOpF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,pv_mass, M5,b,c, ParamsF);
FermionActionF StrangeOpDirF (UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,strange_mass,M5,b,c, ParamsDirF);
FermionActionF StrangePauliVillarsOpDirF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,pv_mass, M5,b,c, ParamsDirF);
#if 1
OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF> StrangePseudoFermionBdy(StrangeOpDir,StrangeOp,
StrangeOpDirF,StrangeOpF,
SFRp,500);
OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF> StrangePseudoFermionLocal(StrangePauliVillarsOpDir,StrangeOpDir,
StrangePauliVillarsOpDirF,StrangeOpDirF,
SFRp,500);
OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF> StrangePseudoFermionPVBdy(StrangePauliVillarsOp,StrangePauliVillarsOpDir,
StrangePauliVillarsOpF,StrangePauliVillarsOpDirF,
SFRp,500);
#else
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermionBdy(StrangeOpDir,StrangeOp,SFRp);
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermionLocal(StrangePauliVillarsOpDir,StrangeOpDir,SFRp);
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermionPVBdy(StrangePauliVillarsOp,StrangePauliVillarsOpDir,SFRp);
#endif
Level1.push_back(&StrangePseudoFermionBdy); // ok
Level2.push_back(&StrangePseudoFermionLocal);
Level1.push_back(&StrangePseudoFermionPVBdy); //ok
////////////////////////////////////
// up down action
////////////////////////////////////
std::vector<Real> light_den;
std::vector<Real> light_num;
std::vector<int> dirichlet_den;
std::vector<int> dirichlet_num;
int n_hasenbusch = hasenbusch.size();
light_den.push_back(light_mass); dirichlet_den.push_back(0);
for(int h=0;h<n_hasenbusch;h++){
light_den.push_back(hasenbusch[h]); dirichlet_den.push_back(1);
}
for(int h=0;h<n_hasenbusch;h++){
light_num.push_back(hasenbusch[h]); dirichlet_num.push_back(1);
}
light_num.push_back(pv_mass); dirichlet_num.push_back(0);
std::vector<FermionAction *> Numerators;
std::vector<FermionActionF *> NumeratorsF;
std::vector<FermionAction *> Denominators;
std::vector<FermionActionF *> DenominatorsF;
std::vector<TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy> *> Quotients;
#define MIXED_PRECISION
#ifdef MIXED_PRECISION
std::vector<OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF> *> Bdys;
#else
std::vector<OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> *> Bdys;
#endif
std::vector<MxPCG *> ActionMPCG;
std::vector<MxPCG *> MPCG;
typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
typedef SchurDiagMooeeOperator<FermionAction ,FermionField > LinearOperatorD;
std::vector<LinearOperatorD *> LinOpD;
std::vector<LinearOperatorF *> LinOpF;
for(int h=0;h<n_hasenbusch+1;h++){
std::cout << GridLogMessage
<< " 2f quotient Action ";
std::cout << "det D("<<light_den[h]<<")";
if ( dirichlet_den[h] ) std::cout << "^dirichlet ";
std::cout << "/ det D("<<light_num[h]<<")";
if ( dirichlet_num[h] ) std::cout << "^dirichlet ";
std::cout << std::endl;
FermionAction::ImplParams ParamsNum(boundary);
FermionAction::ImplParams ParamsDen(boundary);
FermionActionF::ImplParams ParamsNumF(boundary);
FermionActionF::ImplParams ParamsDenF(boundary);
if ( dirichlet_num[h]==1) ParamsNum.dirichlet = Dirichlet;
else ParamsNum.dirichlet = NonDirichlet;
Numerators.push_back (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, ParamsNum));
if ( dirichlet_den[h]==1) ParamsDen.dirichlet = Dirichlet;
else ParamsDen.dirichlet = NonDirichlet;
Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, ParamsDen));
ParamsDenF.dirichlet = ParamsDen.dirichlet;
DenominatorsF.push_back(new FermionActionF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_den[h],M5,b,c, ParamsDenF));
ParamsNumF.dirichlet = ParamsNum.dirichlet;
NumeratorsF.push_back (new FermionActionF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_num[h],M5,b,c, ParamsNumF));
LinOpD.push_back(new LinearOperatorD(*Denominators[h]));
LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h]));
double conv = MDStoppingCondition;
if (h<3) conv= MDStoppingConditionLoose; // Relax on first two hasenbusch factors
const int MX_inner = 5000;
MPCG.push_back(new MxPCG(conv,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
ActionMPCG.push_back(new MxPCG(StoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
if(h!=0) {
// Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],MDCG,CG));
Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],CG));
} else {
#ifdef MIXED_PRECISION
Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF>(
*Numerators[h],*Denominators[h],
*NumeratorsF[h],*DenominatorsF[h],
OFRp, 500) );
Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF>(
*Numerators[h],*Denominators[h],
*NumeratorsF[h],*DenominatorsF[h],
OFRp, 500) );
#else
Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
#endif
}
}
int nquo=Quotients.size();
Level1.push_back(Bdys[0]);
Level1.push_back(Bdys[1]);
for(int h=0;h<nquo-1;h++){
Level2.push_back(Quotients[h]);
}
Level2.push_back(Quotients[nquo-1]);
/////////////////////////////////////////////////////////////
// Gauge action
/////////////////////////////////////////////////////////////
Level3.push_back(&GaugeAction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
TheHMC.TheAction.push_back(Level3);
std::cout << GridLogMessage << " Action complete "<< std::endl;
/////////////////////////////////////////////////////////////
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

6
benchmarks/Benchmark_dwf.cc
View File

@ -167,7 +167,7 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage<< "* Kernel options --dslash-generic, --dslash-unroll, --dslash-asm" <<std::endl; std::cout << GridLogMessage<< "* Kernel options --dslash-generic, --dslash-unroll, --dslash-asm" <<std::endl;
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl; std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl; std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::Dhop "<<std::endl; std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionD::Dhop "<<std::endl;
std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl; std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl;
std::cout << GridLogMessage<< "* VComplex size is "<<sizeof(vComplex)<< " B"<<std::endl; std::cout << GridLogMessage<< "* VComplex size is "<<sizeof(vComplex)<< " B"<<std::endl;
if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl; if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
@ -181,7 +181,7 @@ int main (int argc, char ** argv)
if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3 WilsonKernels" <<std::endl; if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3 WilsonKernels" <<std::endl;
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl; std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
int ncall =1000; int ncall =1000;
if (1) { if (1) {
@ -291,7 +291,7 @@ int main (int argc, char ** argv)
// S-direction is INNERMOST and takes no part in the parity. // S-direction is INNERMOST and takes no part in the parity.
std::cout << GridLogMessage<< "*********************************************************" <<std::endl; std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::DhopEO "<<std::endl; std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionD::DhopEO "<<std::endl;
std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl; std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl;
if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl; if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
if ( sizeof(Real)==8 ) std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl; if ( sizeof(Real)==8 ) std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;

2
benchmarks/Benchmark_dwf_fp32.cc
View File

@ -125,7 +125,7 @@ void Benchmark(int Ls, Coordinate Dirichlet)
std::vector<int> seeds4({1,2,3,4}); std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8}); std::vector<int> seeds5({5,6,7,8});
#define DOUBLE #define SINGLE
#ifdef SINGLE #ifdef SINGLE
typedef vComplexF Simd; typedef vComplexF Simd;
typedef LatticeFermionF FermionField; typedef LatticeFermionF FermionField;

2
benchmarks/Benchmark_dwf_sweep.cc
View File

@ -168,7 +168,7 @@ void benchDw(std::vector<int> & latt4, int Ls, int threads,int report )
RealD M5 =1.8; RealD M5 =1.8;
RealD NP = UGrid->_Nprocessors; RealD NP = UGrid->_Nprocessors;
DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
double t0=usecond(); double t0=usecond();
Dw.Dhop(src,result,0); Dw.Dhop(src,result,0);

6
benchmarks/Benchmark_mooee.cc
View File

@ -67,17 +67,17 @@ int main (int argc, char ** argv)
const int ncall=1000; const int ncall=1000;
std::cout << GridLogMessage<< "*********************************************************" <<std::endl; std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::Dhop "<<std::endl; std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionD::Dhop "<<std::endl;
std::cout << GridLogMessage<< "*********************************************************" <<std::endl; std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeFermion src(FGrid); random(RNG5,src); LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); LatticeFermion result(FGrid);
DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
double t0,t1; double t0,t1;
typedef typename DomainWallFermionR::Coeff_t Coeff_t; typedef typename DomainWallFermionD::Coeff_t Coeff_t;
Vector<Coeff_t> diag = Dw.bs; Vector<Coeff_t> diag = Dw.bs;
Vector<Coeff_t> upper= Dw.cs; Vector<Coeff_t> upper= Dw.cs;
Vector<Coeff_t> lower= Dw.cs; Vector<Coeff_t> lower= Dw.cs;

6
benchmarks/Benchmark_staggered.cc
View File

@ -53,8 +53,8 @@ int main (int argc, char ** argv)
pRNG.SeedFixedIntegers(seeds); pRNG.SeedFixedIntegers(seeds);
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}); // pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
typedef typename ImprovedStaggeredFermionR::FermionField FermionField; typedef typename ImprovedStaggeredFermionD::FermionField FermionField;
typename ImprovedStaggeredFermionR::ImplParams params; typename ImprovedStaggeredFermionD::ImplParams params;
FermionField src (&Grid); random(pRNG,src); FermionField src (&Grid); random(pRNG,src);
FermionField result(&Grid); result=Zero(); FermionField result(&Grid); result=Zero();
@ -93,7 +93,7 @@ int main (int argc, char ** argv)
RealD c1=9.0/8.0; RealD c1=9.0/8.0;
RealD c2=-1.0/24.0; RealD c2=-1.0/24.0;
RealD u0=1.0; RealD u0=1.0;
ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params); ImprovedStaggeredFermionD Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
std::cout<<GridLogMessage << "Calling Ds"<<std::endl; std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
int ncall=1000; int ncall=1000;

4
benchmarks/Benchmark_wilson.cc
View File

@ -146,9 +146,9 @@ int main (int argc, char ** argv)
ref = -0.5*ref; ref = -0.5*ref;
RealD mass=0.1; RealD mass=0.1;
typename WilsonFermionR::ImplParams params; typename WilsonFermionD::ImplParams params;
WilsonFermionR Dw(Umu,Grid,RBGrid,mass,params); WilsonFermionD Dw(Umu,Grid,RBGrid,mass,params);
std::cout<<GridLogMessage << "Calling Dw"<<std::endl; std::cout<<GridLogMessage << "Calling Dw"<<std::endl;
int ncall=1000; int ncall=1000;

14
benchmarks/Benchmark_wilson_sweep.cc
View File

@ -40,21 +40,21 @@ Gamma::Algebra Gmu [] = {
void bench_wilson ( void bench_wilson (
LatticeFermion & src, LatticeFermion & src,
LatticeFermion & result, LatticeFermion & result,
WilsonFermionR & Dw, WilsonFermionD & Dw,
double const volume, double const volume,
int const dag ); int const dag );
void bench_wilson_eo ( void bench_wilson_eo (
LatticeFermion & src, LatticeFermion & src,
LatticeFermion & result, LatticeFermion & result,
WilsonFermionR & Dw, WilsonFermionD & Dw,
double const volume, double const volume,
int const dag ); int const dag );
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
typename WilsonFermionR::ImplParams params; typename WilsonFermionD::ImplParams params;
Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi(); Coordinate mpi_layout = GridDefaultMpi();
@ -66,7 +66,7 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl; std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
std::cout << GridLogMessage<< "*****************************************************************" <<std::endl; std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
std::cout << GridLogMessage<< "* Number of colours "<< Nc <<std::endl; std::cout << GridLogMessage<< "* Number of colours "<< Nc <<std::endl;
std::cout << GridLogMessage<< "* Benchmarking WilsonFermionR::Dhop "<<std::endl; std::cout << GridLogMessage<< "* Benchmarking WilsonFermionD::Dhop "<<std::endl;
std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl; std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl;
if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl; if ( sizeof(Real)==4 ) std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
if ( sizeof(Real)==8 ) std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl; if ( sizeof(Real)==8 ) std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;
@ -110,7 +110,7 @@ int main (int argc, char ** argv)
double volume = std::accumulate(latt_size.begin(),latt_size.end(),1,std::multiplies<int>()); double volume = std::accumulate(latt_size.begin(),latt_size.end(),1,std::multiplies<int>());
WilsonFermionR Dw(Umu,Grid,RBGrid,mass,params); WilsonFermionD Dw(Umu,Grid,RBGrid,mass,params);
// Full operator // Full operator
bench_wilson(src,result,Dw,volume,DaggerNo); bench_wilson(src,result,Dw,volume,DaggerNo);
@ -130,7 +130,7 @@ int main (int argc, char ** argv)
void bench_wilson ( void bench_wilson (
LatticeFermion & src, LatticeFermion & src,
LatticeFermion & result, LatticeFermion & result,
WilsonFermionR & Dw, WilsonFermionD & Dw,
double const volume, double const volume,
int const dag ) int const dag )
{ {
@ -149,7 +149,7 @@ void bench_wilson (
void bench_wilson_eo ( void bench_wilson_eo (
LatticeFermion & src, LatticeFermion & src,
LatticeFermion & result, LatticeFermion & result,
WilsonFermionR & Dw, WilsonFermionD & Dw,
double const volume, double const volume,
int const dag ) int const dag )
{ {

4
examples/Example_Mobius_spectrum.cc
View File

@ -253,7 +253,7 @@ int main (int argc, char ** argv)
int nmass = masses.size(); int nmass = masses.size();
std::vector<MobiusFermionR *> FermActs; std::vector<MobiusFermionD *> FermActs;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl; std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
@ -265,7 +265,7 @@ int main (int argc, char ** argv)
RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5; RealD c=0.5;
FermActs.push_back(new MobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c)); FermActs.push_back(new MobiusFermionD(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c));
} }

8
examples/Example_christoph.cc
View File

@ -387,8 +387,8 @@ int main (int argc, char ** argv)
int nmass = masses.size(); int nmass = masses.size();
typedef MobiusFermionR FermionActionR; typedef MobiusFermionD FermionActionD;
std::vector<FermionActionR *> FermActs; std::vector<FermionActionD *> FermActs;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"DomainWallFermion action"<<std::endl; std::cout<<GridLogMessage <<"DomainWallFermion action"<<std::endl;
@ -396,10 +396,10 @@ int main (int argc, char ** argv)
for(auto mass: masses) { for(auto mass: masses) {
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionActionR::ImplParams Params(boundary); FermionActionD::ImplParams Params(boundary);
RealD b=1.5; RealD b=1.5;
RealD c=0.5; RealD c=0.5;
FermActs.push_back(new FermionActionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c)); FermActs.push_back(new FermionActionD(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c));
} }
LatticePropagator point_source(UGrid); LatticePropagator point_source(UGrid);

16
examples/Example_taku1.cc
View File

@ -416,10 +416,10 @@ int main (int argc, char ** argv)
int nmass = masses.size(); int nmass = masses.size();
typedef DomainWallFermionR FermionActionR; typedef DomainWallFermionD FermionActionD;
// typedef MobiusFermionR FermionActionR; // typedef MobiusFermionD FermionActionD;
std::vector<FermionActionR *> FermActs; std::vector<FermionActionD *> FermActs;
std::vector<DomainWallFermionR *> DWFActs; std::vector<DomainWallFermionD *> DWFActs;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"DomainWallFermion action"<<std::endl; std::cout<<GridLogMessage <<"DomainWallFermion action"<<std::endl;
@ -427,13 +427,13 @@ int main (int argc, char ** argv)
for(auto mass: masses) { for(auto mass: masses) {
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionActionR::ImplParams Params(boundary); FermionActionD::ImplParams Params(boundary);
RealD b=1.5; RealD b=1.5;
RealD c=0.5; RealD c=0.5;
std::cout<<GridLogMessage <<"Making DomainWallFermion action"<<std::endl; std::cout<<GridLogMessage <<"Making DomainWallFermion action"<<std::endl;
// DWFActs.push_back(new DomainWallFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5)); // DWFActs.push_back(new DomainWallFermionD(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5));
FermActs.push_back(new FermionActionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,Params)); FermActs.push_back(new FermionActionD(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,Params));
// FermActs.push_back(new FermionActionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass+0.001,M5,b,c)); // FermActs.push_back(new FermionActionD(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass+0.001,M5,b,c));
std::cout<<GridLogMessage <<"Made DomainWallFermion action"<<std::endl; std::cout<<GridLogMessage <<"Made DomainWallFermion action"<<std::endl;
} }

6
examples/Example_wall_wall_3pt.cc
View File

@ -389,14 +389,14 @@ int main (int argc, char ** argv)
int nmass = masses.size(); int nmass = masses.size();
int nmom = momenta.size(); int nmom = momenta.size();
std::vector<MobiusFermionR *> FermActs; std::vector<MobiusFermionD *> FermActs;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl; std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
typedef MobiusFermionR FermionAction; typedef MobiusFermionD FermionAction;
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
for(int m=0;m<masses.size();m++) { for(int m=0;m<masses.size();m++) {
@ -413,7 +413,7 @@ int main (int argc, char ** argv)
FGrids.push_back(SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid)); FGrids.push_back(SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid));
FrbGrids.push_back(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid)); FrbGrids.push_back(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid));
FermActs.push_back(new MobiusFermionR(Utmp,*FGrids[m],*FrbGrids[m],*UGrid,*UrbGrid,mass,M5,b,c,Params)); FermActs.push_back(new MobiusFermionD(Utmp,*FGrids[m],*FrbGrids[m],*UGrid,*UrbGrid,mass,M5,b,c,Params));
} }
LatticePropagator z2wall_source(UGrid); LatticePropagator z2wall_source(UGrid);

6
examples/Example_wall_wall_spectrum.cc
View File

@ -340,13 +340,13 @@ int main (int argc, char ** argv)
int nmass = masses.size(); int nmass = masses.size();
std::vector<MobiusFermionR *> FermActs; std::vector<MobiusFermionD *> FermActs;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl; std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
std::cout<<GridLogMessage <<"======================"<<std::endl; std::cout<<GridLogMessage <<"======================"<<std::endl;
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
typedef MobiusFermionR FermionAction; typedef MobiusFermionD FermionAction;
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
for(int m=0;m<masses.size();m++) { for(int m=0;m<masses.size();m++) {
@ -363,7 +363,7 @@ int main (int argc, char ** argv)
FGrids.push_back(SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid)); FGrids.push_back(SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid));
FrbGrids.push_back(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid)); FrbGrids.push_back(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid));
FermActs.push_back(new MobiusFermionR(Utmp,*FGrids[m],*FrbGrids[m],*UGrid,*UrbGrid,mass,M5,b,c,Params)); FermActs.push_back(new MobiusFermionD(Utmp,*FGrids[m],*FrbGrids[m],*UGrid,*UrbGrid,mass,M5,b,c,Params));
} }
LatticePropagator z2wall_source(UGrid); LatticePropagator z2wall_source(UGrid);

32
tests/Test_cayley_even_odd_vec.cc
View File

@ -88,8 +88,8 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"DomainWallFermion vectorised test"<<std::endl; std::cout<<GridLogMessage <<"DomainWallFermion vectorised test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
TestWhat<DomainWallFermionR>(Ddwf,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<DomainWallFermionD>(Ddwf,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5; RealD c=0.5;
@ -97,54 +97,54 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl; std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); MobiusFermionD Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestWhat<MobiusFermionR>(Dmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<MobiusFermionD>(Dmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"Z-MobiusFermion test"<<std::endl; std::cout<<GridLogMessage <<"Z-MobiusFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::vector<ComplexD> gamma(Ls,std::complex<double>(1.0,0.0)); std::vector<ComplexD> gamma(Ls,std::complex<double>(1.0,0.0));
ZMobiusFermionR zDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c); ZMobiusFermionD zDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c);
TestWhat<ZMobiusFermionR>(zDmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<ZMobiusFermionD>(zDmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl; std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0); MobiusZolotarevFermionD Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
TestWhat<MobiusZolotarevFermionR>(Dzolo,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<MobiusZolotarevFermionD>(Dzolo,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl; std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0); ScaledShamirFermionD Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
TestWhat<ScaledShamirFermionR>(Dsham,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<ScaledShamirFermionD>(Dsham,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl; std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
ShamirZolotarevFermionR Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); ShamirZolotarevFermionD Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestWhat<ShamirZolotarevFermionR>(Dshamz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<ShamirZolotarevFermionD>(Dshamz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl; std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
OverlapWilsonCayleyTanhFermionR Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); OverlapWilsonCayleyTanhFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestWhat<OverlapWilsonCayleyTanhFermionR>(Dov,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<OverlapWilsonCayleyTanhFermionD>(Dov,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl; std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************************"<<std::endl; std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
OverlapWilsonCayleyZolotarevFermionR Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); OverlapWilsonCayleyZolotarevFermionD Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestWhat<OverlapWilsonCayleyZolotarevFermionR>(Dovz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); TestWhat<OverlapWilsonCayleyZolotarevFermionD>(Dovz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage<<"=============================================================="<<std::endl; std::cout<<GridLogMessage<<"=============================================================="<<std::endl;

4
tests/Test_compressed_lanczos_hot_start.cc
View File

@ -209,8 +209,8 @@ int main (int argc, char ** argv) {
std::cout << GridLogMessage << "Lattice dimensions: " << latt << " Ls: " << Ls << std::endl; std::cout << GridLogMessage << "Lattice dimensions: " << latt << " Ls: " << Ls << std::endl;
// ZMobius EO Operator // ZMobius EO Operator
ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, Params.omega,1.,0.); ZMobiusFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, Params.omega,1.,0.);
SchurDiagTwoOperator<ZMobiusFermionR,LatticeFermion> HermOp(Ddwf); SchurDiagTwoOperator<ZMobiusFermionD,LatticeFermion> HermOp(Ddwf);
// Eigenvector storage // Eigenvector storage
LanczosParams fine =Params.FineParams; LanczosParams fine =Params.FineParams;

183
tests/Test_gfield_shift.cc Normal file
View File

@ -0,0 +1,183 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_gfield_shift.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
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 */
//Test the shifting of the gauge field that respects the boundary conditions
#include <Grid/Grid.h>
using namespace Grid;
;
typedef ConjugateGimplR Gimpl; //can choose periodic / charge conjugate directions at wil
typedef Gimpl::GaugeField GaugeField;
typedef Gimpl::GaugeLinkField GaugeLinkField;
typedef Gimpl::SiteGaugeField SiteGaugeField;
typedef Gimpl::SiteGaugeLink SiteGaugeLink;
GaugeField CshiftGaugeField(const GaugeField &U, const int dir, const int shift){
GridBase *Grid = U.Grid();
GaugeField out(Grid);
GaugeLinkField Umu(Grid);
for(int mu=0;mu<Grid->Nd();mu++){
Umu = PeekIndex<LorentzIndex>(U, mu);
Umu = Gimpl::CshiftLink(Umu,dir,shift);
PokeIndex<LorentzIndex>(out,Umu,mu);
}
return out;
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
auto latt_size = GridDefaultLatt();
auto simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
auto mpi_layout = GridDefaultMpi();
std::vector<int> conj_dirs = {1,1,0,0};
Gimpl::setDirections(conj_dirs);
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridParallelRNG FineRNG(&Fine); FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
GaugeField U(&Fine);
GaugeField ShiftU(&Fine);
GaugeLinkField link_field(&Fine), link_field_2(&Fine);
//Like Test_cshift we put the lex coordinate index on each site but make it imaginary
//so we can tell when it was complex conjugated
LatticeComplex lex(&Fine);
lex=Zero();
U = Zero();
{
LatticeComplex coor(&Fine);
Integer stride =1;
for(int d=0;d<4;d++){
LatticeCoordinate(coor,d);
lex = lex + coor*stride;
stride=stride*latt_size[d];
}
PokeIndex<ColourIndex>(link_field, lex, 0,0); //place on 0,0 element of link
for(int mu=0;mu<Nd;mu++){
link_field_2 = link_field + mu*stride; //add in lex-mapping of mu
link_field_2 = ComplexD(0,1) * link_field_2; //make imaginary
PokeIndex<LorentzIndex>(U, link_field_2, mu);
}
}
std::stringstream ss;
ss<<"error";
for(int d=0;d<Fine._ndimension;d++){
ss<<"."<<Fine._processor_coor[d];
}
ss<<"_wr_"<<Fine._processor;
std::string fname(ss.str());
std::ofstream ferr(fname);
Integer vol4d = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
bool fail = false;
typename SiteGaugeField::scalar_object um;
TComplex cm;
for(int dir=0;dir<4;dir++){
for(int shift=-latt_size[dir]+1;shift<latt_size[dir];shift++){
if ( Fine.IsBoss() )
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction "<<dir
<< " dir is conj ? " << conj_dirs[dir] << std::endl;
ShiftU = CshiftGaugeField(U,dir,shift);
Coordinate coor(4);
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
for(coor[0]=0;coor[0]<latt_size[0];coor[0]++){
peekSite(um,ShiftU,coor);
Coordinate scoor(coor);
scoor[dir] = (scoor[dir]+shift + latt_size[dir])%latt_size[dir];
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
for(int mu = 0 ; mu < 4; mu++){
Integer slex_mu = slex + vol4d*mu;
Complex scm(0,slex_mu); //imaginary
if(
( shift > 0 && coor[dir] >= latt_size[dir]-shift && conj_dirs[dir] )
||
( shift < 0 && coor[dir] <= -shift-1 && conj_dirs[dir] )
)
scm = conjugate(scm); //CC if pulled over boundary
cm = um(mu)()(0,0);
RealD nrm = abs(scm-cm()()());
//std::cout << cm << " " << scm << std::endl;
Coordinate peer(4);
Complex tmp =cm;
Integer index=real(tmp);
Integer cm_mu = index / vol4d;
index = index % vol4d;
Lexicographic::CoorFromIndex(peer,index,latt_size);
if (nrm > 0){
ferr<<"FAIL mu " << mu << " shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
ferr<<"Got mu "<< cm_mu << " site " <<index<<" : " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Integer scm_mu = index / vol4d;
index = index % vol4d;
Lexicographic::CoorFromIndex(peer,index,latt_size);
ferr<<"Expect mu " << scm_mu << " site " <<index<<": " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
fail = true;
}
}
}}}}
}
}
if(fail) std::cout << "Test FAILED : see " << fname << " for more details" << std::endl;
else std::cout << "Test Passed" << std::endl;
Grid_finalize();
}

4
tests/Test_meson_field.cc
View File

@ -33,8 +33,8 @@ using namespace Grid;
const int TSRC = 0; //timeslice where rho is nonzero const int TSRC = 0; //timeslice where rho is nonzero
const int VDIM = 5; //length of each vector const int VDIM = 5; //length of each vector
typedef typename DomainWallFermionR::ComplexField ComplexField; typedef typename DomainWallFermionD::ComplexField ComplexField;
typedef typename DomainWallFermionR::FermionField FermionField; typedef typename DomainWallFermionD::FermionField FermionField;
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {

153
tests/smearing/Test_WilsonFlow_adaptive.cc Normal file
View File

@ -0,0 +1,153 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/hmc/Test_WilsonFlow_adaptive.cc
Copyright (C) 2017
Author: Christopher Kelly <ckelly@bnl.gov>
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;
//Linearly interpolate between two nearest times
RealD interpolate(const RealD t_int, const std::vector<std::pair<RealD,RealD> > &data){
RealD tdiff1=1e32; int t1_idx=-1;
RealD tdiff2=1e32; int t2_idx=-1;
for(int i=0;i<data.size();i++){
RealD diff = fabs(data[i].first-t_int);
//std::cout << "targ " << t_int << " cur " << data[i].first << " diff " << diff << " best diff1 " << tdiff1 << " diff2 " << tdiff2 << std::endl;
if(diff < tdiff1){
if(tdiff1 < tdiff2){ //swap out tdiff2
tdiff2 = tdiff1; t2_idx = t1_idx;
}
tdiff1 = diff; t1_idx = i;
}
else if(diff < tdiff2){ tdiff2 = diff; t2_idx = i; }
}
assert(t1_idx != -1 && t2_idx != -1);
RealD t2 = data[t2_idx].first, v2 = data[t2_idx].second;
RealD t1 = data[t1_idx].first, v1 = data[t1_idx].second;
//v = a + bt
//v2-v1 = b(t2-t1)
RealD b = (v2-v1)/(t2-t1);
RealD a = v1 - b*t1;
RealD vout = a + b*t_int;
//std::cout << "Interpolate to " << t_int << " two closest points " << t1 << " " << t2
//<< " with values " << v1 << " "<< v2 << " : got " << vout << std::endl;
return vout;
}
int main(int argc, char **argv) {
Grid_init(&argc, &argv);
GridLogLayout();
auto latt_size = GridDefaultLatt();
auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
auto mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size, simd_layout, mpi_layout);
GridRedBlackCartesian RBGrid(&Grid);
std::vector<int> seeds({1, 2, 3, 4, 5});
GridSerialRNG sRNG;
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(seeds);
LatticeGaugeField U(&Grid);
SU<Nc>::HotConfiguration(pRNG, U);
int Nstep = 300;
RealD epsilon = 0.01;
RealD maxTau = Nstep*epsilon;
RealD tolerance = 1e-4;
for(int i=1;i<argc;i++){
std::string sarg(argv[i]);
if(sarg == "--tolerance"){
std::stringstream ss; ss << argv[i+1]; ss >> tolerance;
}
}
std::cout << "Adaptive smear tolerance " << tolerance << std::endl;
//Setup iterative Wilson flow
WilsonFlow<PeriodicGimplD> wflow(epsilon,Nstep);
wflow.resetActions();
std::vector<std::pair<RealD, RealD> > meas_orig;
wflow.addMeasurement(1, [&wflow,&meas_orig](int step, RealD t, const LatticeGaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl;
meas_orig.push_back( {t, wflow.energyDensityCloverleaf(t,U)} );
});
//Setup adaptive Wilson flow
WilsonFlowAdaptive<PeriodicGimplD> wflow_ad(epsilon,maxTau,tolerance);
wflow_ad.resetActions();
std::vector<std::pair<RealD, RealD> > meas_adaptive;
wflow_ad.addMeasurement(1, [&wflow_ad,&meas_adaptive](int step, RealD t, const LatticeGaugeField &U){
std::cout << GridLogMessage << "[WilsonFlow] Computing Cloverleaf energy density for step " << step << std::endl;
meas_adaptive.push_back( {t, wflow_ad.energyDensityCloverleaf(t,U)} );
});
//Run
LatticeGaugeFieldD Vtmp(U.Grid());
wflow.smear(Vtmp, U); //basic smear
Vtmp = Zero();
wflow_ad.smear(Vtmp, U);
//Output values for plotting
{
std::ofstream out("wflow_t2E_orig.dat");
out.precision(16);
for(auto const &e: meas_orig){
out << e.first << " " << e.second << std::endl;
}
}
{
std::ofstream out("wflow_t2E_adaptive.dat");
out.precision(16);
for(auto const &e: meas_adaptive){
out << e.first << " " << e.second << std::endl;
}
}
//Compare at times available with adaptive smearing
for(int i=0;i<meas_adaptive.size();i++){
RealD t = meas_adaptive[i].first;
RealD v_adaptive = meas_adaptive[i].second;
RealD v_orig = interpolate(t, meas_orig); //should be very precise due to fine timestep
std::cout << t << " orig: " << v_orig << " adaptive: " << v_adaptive << " reldiff: " << (v_adaptive-v_orig)/v_orig << std::endl;
}
std::cout << GridLogMessage << "Done" << std::endl;
Grid_finalize();
}