diff --git a/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h b/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h index 8a88cb07..de1cfe01 100644 --- a/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h +++ b/Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h @@ -154,6 +154,7 @@ public: // dynamic sized arrays on stack; 2d is a pain with vector RealD bs[nshift]; RealD rsq[nshift]; + RealD rsqf[nshift]; RealD z[nshift][2]; int converged[nshift]; @@ -190,6 +191,7 @@ public: for(int s=0;s &in, Lattice &out) template void precisionChangeFast(Lattice &out, const Lattice &in) { - typedef typename VobjOut::scalar_object SobjOut; - typedef typename VobjIn::scalar_object SobjIn; + typedef typename VobjOut::vector_type Vout; + typedef typename VobjIn::vector_type Vin; + const int N = sizeof(VobjOut)/sizeof(Vout); conformable(out.Grid(),in.Grid()); out.Checkerboard() = in.Checkerboard(); int nsimd = out.Grid()->Nsimd(); autoView( out_v , out, AcceleratorWrite); autoView( in_v , in, AcceleratorRead); - accelerator_for(idx,out.Grid()->oSites(),nsimd,{ - auto itmp = coalescedRead(in_v[idx]); - auto otmp = coalescedRead(out_v[idx]); -#ifdef GRID_SIMT - otmp=itmp; -#endif - coalescedWrite(out_v[idx],otmp); + accelerator_for(idx,out.Grid()->oSites(),1,{ + Vout *vout = (Vout *)&out_v[idx]; + Vin *vin = (Vin *)&in_v[idx]; + precisionChange(vout,vin,N); }); } //Convert a Lattice from one precision to another @@ -1116,7 +1114,7 @@ void precisionChange(Lattice &out, const Lattice &in) int ndim = out.Grid()->Nd(); int out_nsimd = out_grid->Nsimd(); - + int in_nsimd = in_grid->Nsimd(); std::vector out_icoor(out_nsimd); for(int lane=0; lane < out_nsimd; lane++){ diff --git a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h index d154fda5..cd8c4d64 100644 --- a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h +++ b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h @@ -36,14 +36,18 @@ NAMESPACE_BEGIN(Grid); // cf. GeneralEvenOddRational.h for details ///////////////////////////////////////////////////////////////////////////////////////////////////////////// - template + template class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction { private: + typedef typename ImplD2::FermionField FermionFieldD2; typedef typename ImplD::FermionField FermionFieldD; typedef typename ImplF::FermionField FermionFieldF; FermionOperator & NumOpD; FermionOperator & DenOpD; + + FermionOperator & NumOpD2; + FermionOperator & DenOpD2; FermionOperator & NumOpF; FermionOperator & DenOpF; @@ -53,37 +57,70 @@ NAMESPACE_BEGIN(Grid); //Allow derived classes to override the multishift CG virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){ - SchurDifferentiableOperator schurOpD(numerator ? NumOpD : DenOpD); +#if 0 + SchurDifferentiableOperator schurOp(numerator ? NumOp : DenOp); + ConjugateGradientMultiShift msCG(MaxIter, approx); + msCG(schurOp,in, out); +#else + SchurDifferentiableOperator schurOpD2(numerator ? NumOpD2 : DenOpD2); SchurDifferentiableOperator schurOpF(numerator ? NumOpF : DenOpF); - - ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); - msCG(schurOpD, in, out); + FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid()); + FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid()); + + ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); + precisionChange(inD2,in); + std::cout << "msCG single solve "< &out_elems, FermionFieldD &out){ - SchurDifferentiableOperator schurOpD(numerator ? NumOpD : DenOpD); + SchurDifferentiableOperator schurOpD2(numerator ? NumOpD2 : DenOpD2); SchurDifferentiableOperator schurOpF(numerator ? NumOpF : DenOpF); - ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); - msCG(schurOpD, in, out_elems, out); + FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid()); + FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid()); + std::vector out_elemsD2(out_elems.size(),NumOpD2.FermionRedBlackGrid()); + ConjugateGradientMultiShiftMixedPrec msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq); + precisionChange(inD2,in); + std::cout << "msCG in "< &_NumOpD, FermionOperator &_DenOpD, FermionOperator &_NumOpF, FermionOperator &_DenOpF, + FermionOperator &_NumOpD2, FermionOperator &_DenOpD2, const RationalActionParams & p, Integer _ReliableUpdateFreq ) : GeneralEvenOddRatioRationalPseudoFermionAction(_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";} }; diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index 078bf845..302c3a5d 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -67,8 +67,9 @@ NAMESPACE_BEGIN(Grid); virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} }; - template - class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction { + template + class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction + : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction { public: typedef OneFlavourRationalParams Params; private: @@ -90,9 +91,11 @@ NAMESPACE_BEGIN(Grid); FermionOperator &_DenOp, FermionOperator &_NumOpF, FermionOperator &_DenOpF, + FermionOperator &_NumOpD2, + FermionOperator &_DenOpD2, const Params & p, Integer ReliableUpdateFreq ) : - GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(_NumOp, _DenOp,_NumOpF, _DenOpF, transcribe(p),ReliableUpdateFreq){} + GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(_NumOp, _DenOp,_NumOpF, _DenOpF,_NumOpD2, _DenOpD2, transcribe(p),ReliableUpdateFreq){} virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";} }; diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index 86796f05..a9fe6f46 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -53,6 +53,7 @@ struct HMCparameters: Serializable { Integer, Trajectories, /* @brief Number of sweeps in this run */ bool, MetropolisTest, Integer, NoMetropolisUntil, + bool, PerformRandomShift, /* @brief Randomly shift the gauge configuration at the start of a trajectory */ std::string, StartingType, IntegratorParameters, MD) @@ -63,6 +64,7 @@ struct HMCparameters: Serializable { StartTrajectory = 0; Trajectories = 10; 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] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\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"; MD.print_parameters(); } @@ -95,6 +98,7 @@ private: const HMCparameters Params; typedef typename IntegratorType::Field Field; + typedef typename IntegratorType::FieldImplementation FieldImplementation; typedef std::vector< HmcObservable * > ObsListType; //pass these from the resource manager @@ -138,26 +142,37 @@ private: GridBase *Grid = U.Grid(); - ////////////////////////////////////////////////////////////////////////////////////////////////////// - // Mainly for DDHMC perform a random translation of U modulo volume - ////////////////////////////////////////////////////////////////////////////////////////////////////// - std::cout << GridLogMessage << "--------------------------------------------------\n"; - std::cout << GridLogMessage << "Random shifting gauge field by ["; - for(int d=0;dNd();d++) { + if(Params.PerformRandomShift){ + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // Mainly for DDHMC perform a random translation of U modulo volume + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Random shifting gauge field by ["; - int L = Grid->GlobalDimensions()[d]; + std::vector Umu(Grid->Nd(), U.Grid()); + for(int mu=0;muNd();mu++) Umu[mu] = PeekIndex(U, mu); - RealD rn_uniform; random(sRNG, rn_uniform); + for(int d=0;dNd();d++) { - int shift = (int) (rn_uniform*L); + int L = Grid->GlobalDimensions()[d]; - std::cout << shift; - if(dNd()-1) std::cout <<","; - else std::cout <<"]\n"; + RealD rn_uniform; random(sRNG, rn_uniform); + + int shift = (int) (rn_uniform*L); + + std::cout << shift; + if(dNd()-1) std::cout <<","; + 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;muNd();mu++) PokeIndex(U,Umu[mu],mu); + + std::cout << GridLogMessage << "--------------------------------------------------\n"; } - std::cout << GridLogMessage << "--------------------------------------------------\n"; TheIntegrator.reset_timer(); diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 48219c04..016b1a84 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -63,10 +63,10 @@ public: }; /*! @brief Class for Molecular Dynamics management */ -template +template class Integrator { protected: - + typedef FieldImplementation_ FieldImplementation; typedef typename FieldImplementation::Field MomentaField; //for readability typedef typename FieldImplementation::Field Field; diff --git a/Grid/qcd/hmc/integrators/Integrator_algorithm.h b/Grid/qcd/hmc/integrators/Integrator_algorithm.h index ad556659..9c70fd1f 100644 --- a/Grid/qcd/hmc/integrators/Integrator_algorithm.h +++ b/Grid/qcd/hmc/integrators/Integrator_algorithm.h @@ -92,10 +92,11 @@ NAMESPACE_BEGIN(Grid); * P 1/2 P 1/2 */ -template > -class LeapFrog : public Integrator +template > +class LeapFrog : public Integrator { public: + typedef FieldImplementation_ FieldImplementation; typedef LeapFrog Algorithm; INHERIT_FIELD_TYPES(FieldImplementation); @@ -135,13 +136,14 @@ public: } }; -template > -class MinimumNorm2 : public Integrator +template > +class MinimumNorm2 : public Integrator { private: const RealD lambda = 0.1931833275037836; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) @@ -192,8 +194,8 @@ public: } }; -template > -class ForceGradient : public Integrator +template > +class ForceGradient : public Integrator { private: const RealD lambda = 1.0 / 6.0; @@ -202,6 +204,7 @@ private: const RealD theta = 0.0; public: + typedef FieldImplementation_ FieldImplementation; INHERIT_FIELD_TYPES(FieldImplementation); // Looks like dH scales as dt^4. tested wilson/wilson 2 level. diff --git a/Grid/qcd/observables/topological_charge.h b/Grid/qcd/observables/topological_charge.h index 7c09a180..220ed738 100644 --- a/Grid/qcd/observables/topological_charge.h +++ b/Grid/qcd/observables/topological_charge.h @@ -31,15 +31,16 @@ directory NAMESPACE_BEGIN(Grid); + struct TopologySmearingParameters : Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(TopologySmearingParameters, - int, steps, - float, step_size, 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): - steps(s), step_size(ss), meas_interval(mi), maxTau(mT){} + TopologySmearingParameters(float ss = 0.0f, int mi = 0, float mT = 0.0f, float tol = 1e-4): + init_step_size(ss), meas_interval(mi), maxTau(mT), tolerance(tol){} template < class ReaderClass > TopologySmearingParameters(Reader& Reader){ @@ -97,8 +98,8 @@ public: if (Pars.do_smearing){ // using wilson flow by default here - WilsonFlow WF(Pars.Smearing.steps, Pars.Smearing.step_size, Pars.Smearing.meas_interval); - WF.smear_adaptive(Usmear, U, Pars.Smearing.maxTau); + WilsonFlowAdaptive WF(Pars.Smearing.init_step_size, Pars.Smearing.maxTau, Pars.Smearing.tolerance, Pars.Smearing.meas_interval); + WF.smear(Usmear, U); Real T0 = WF.energyDensityPlaquette(Pars.Smearing.maxTau, Usmear); std::cout << GridLogMessage << std::setprecision(std::numeric_limits::digits10 + 1) << "T0 : [ " << traj << " ] "<< T0 << std::endl; diff --git a/Grid/qcd/smearing/WilsonFlow.h b/Grid/qcd/smearing/WilsonFlow.h index 0d1ee5d2..f169d02b 100644 --- a/Grid/qcd/smearing/WilsonFlow.h +++ b/Grid/qcd/smearing/WilsonFlow.h @@ -33,27 +33,25 @@ directory NAMESPACE_BEGIN(Grid); template -class WilsonFlow: public Smear{ +class WilsonFlowBase: public Smear{ public: //Store generic measurements to take during smearing process using std::function typedef std::function FunctionType; //int: step, RealD: flow time, GaugeField : the gauge field - -private: - unsigned int Nstep; - RealD epsilon; //for regular smearing this is the time step, for adaptive it is the initial time step - + +protected: std::vector< std::pair > functions; //The int maps to the measurement frequency mutable WilsonGaugeAction 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: INHERIT_GIMPL_TYPES(Gimpl) + explicit WilsonFlowBase(unsigned int meas_interval =1): + SG(WilsonGaugeAction(3.0)) { + // WilsonGaugeAction with beta 3.0 + setDefaultMeasurements(meas_interval); + } + void resetActions(){ functions.clear(); } void addMeasurement(int meas_interval, FunctionType meas){ functions.push_back({meas_interval, meas}); } @@ -64,34 +62,11 @@ public: //and output to stdout void setDefaultMeasurements(int topq_meas_interval = 1); - explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1): - Nstep(Nstep), - epsilon(epsilon), - SG(WilsonGaugeAction(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 { + void derivative(GaugeField&, const GaugeField&, const GaugeField&) const override{ assert(0); // undefined for WilsonFlow } - void smear_adaptive(GaugeField&, const GaugeField&, RealD maxTau) const; - //Compute t^2 for time t from the plaquette static RealD energyDensityPlaquette(const RealD t, const GaugeField& U); @@ -115,82 +90,63 @@ public: std::vector flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval = 1); }; +//Basic iterative Wilson flow +template +class WilsonFlow: public WilsonFlowBase{ +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(meas_interval), Nstep(Nstep), epsilon(epsilon){} + + void smear(GaugeField& out, const GaugeField& in) const override; +}; + +//Wilson flow with adaptive step size +template +class WilsonFlowAdaptive: public WilsonFlowBase{ +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(meas_interval), init_epsilon(init_epsilon), maxTau(maxTau), tolerance(tolerance){} + + void smear(GaugeField& out, const GaugeField& in) const override; +}; //////////////////////////////////////////////////////////////////////////////// // Implementations //////////////////////////////////////////////////////////////////////////////// template -void WilsonFlow::evolve_step(typename Gimpl::GaugeField &U, RealD &tau) const{ - 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 -void WilsonFlow::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 -RealD WilsonFlow::energyDensityPlaquette(const RealD t, const GaugeField& U){ +RealD WilsonFlowBase::energyDensityPlaquette(const RealD t, const GaugeField& U){ static WilsonGaugeAction SG(3.0); return 2.0 * t * t * SG.S(U)/U.Grid()->gSites(); } //Compute t^2 for time from the 1x1 cloverleaf form template -RealD WilsonFlow::energyDensityCloverleaf(const RealD t, const GaugeField& U){ +RealD WilsonFlowBase::energyDensityCloverleaf(const RealD t, const GaugeField& U){ typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeField GaugeLorentz; @@ -215,7 +171,7 @@ RealD WilsonFlow::energyDensityCloverleaf(const RealD t, const GaugeField template -std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityPlaquette(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector out; resetActions(); addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ @@ -227,13 +183,13 @@ std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(GaugeFie } template -std::vector WilsonFlow::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityPlaquette(const GaugeField& U, int measure_interval){ GaugeField V(U); return flowMeasureEnergyDensityPlaquette(V,U, measure_interval); } template -std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityCloverleaf(GaugeField &V, const GaugeField& U, int measure_interval){ std::vector out; resetActions(); addMeasurement(measure_interval, [&out](int step, RealD t, const typename Gimpl::GaugeField &U){ @@ -245,16 +201,52 @@ std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(GaugeFi } template -std::vector WilsonFlow::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ +std::vector WilsonFlowBase::flowMeasureEnergyDensityCloverleaf(const GaugeField& U, int measure_interval){ GaugeField V(U); return flowMeasureEnergyDensityCloverleaf(V,U, measure_interval); } +template +void WilsonFlowBase::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::TopologicalCharge(U) << std::endl; + }); +} -//#define WF_TIMING + +template +void WilsonFlow::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 void WilsonFlow::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; RealD taus = 0.; for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement @@ -266,37 +258,93 @@ void WilsonFlow::smear(GaugeField& out, const GaugeField& in) const{ std::cout << "Time to evolve " << diff.count() << " s\n"; #endif //Perform measurements - for(auto const &meas : functions) + for(auto const &meas : this->functions) if( step % meas.first == 0 ) meas.second(step,taus,out); } } + + template -void WilsonFlow::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau) const{ - out = in; - RealD taus = 0.; - RealD eps = epsilon; - unsigned int step = 0; - do{ - step++; - //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl; - evolve_step_adaptive(out, taus, eps, maxTau); - //Perform measurements - for(auto const &meas : functions) - if( step % meas.first == 0 ) meas.second(step,taus,out); - } while (taus < maxTau); +int WilsonFlowAdaptive::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(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 -void WilsonFlow::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::TopologicalCharge(U) << std::endl; - }); +void WilsonFlowAdaptive::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; + RealD taus = 0.; + RealD eps = init_epsilon; + unsigned int step = 0; + do{ + int step_success = evolve_step_adaptive(out, taus, eps); + step += step_success; //step will not be incremented if the integration step fails + + //Perform measurements + if(step_success) + for(auto const &meas : this->functions) + if( step % meas.first == 0 ) meas.second(step,taus,out); + } while (taus < maxTau); } + NAMESPACE_END(Grid); diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h index 79cf8e0f..23984145 100644 --- a/Grid/qcd/utils/CovariantCshift.h +++ b/Grid/qcd/utils/CovariantCshift.h @@ -227,26 +227,38 @@ namespace ConjugateBC { //shift = -1 //Out(x) = U_\mu(x-mu) | 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 Lattice CshiftLink(const Lattice &Link, int mu, int shift) { GridBase *grid = Link.Grid(); - int Lmu = grid->GlobalDimensions()[mu] - 1; + int Lmu = grid->GlobalDimensions()[mu]; + assert(abs(shift) < Lmu && "Invalid shift value"); Lattice> coor(grid); LatticeCoordinate(coor, mu); Lattice tmp(grid); - if(shift == 1){ - tmp = Cshift(Link, mu, 1); - tmp = where(coor == Lmu, conjugate(tmp), tmp); + if(shift > 0){ + tmp = Cshift(Link, mu, shift); + tmp = where(coor >= Lmu-shift, conjugate(tmp), tmp); return tmp; - }else if(shift == -1){ + }else if(shift < 0){ tmp = Link; - tmp = where(coor == Lmu, conjugate(tmp), tmp); - return Cshift(tmp, mu, -1); - }else assert(0 && "Invalid shift value"); - return tmp; //shuts up the compiler fussing about the return type + tmp = where(coor >= Lmu+shift, conjugate(tmp), tmp); + return Cshift(tmp, mu, shift); + } + + //shift == 0 + return Link; } } diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index d9d03c54..fc723fe3 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -72,12 +72,12 @@ public: //Fix the gauge field Umu //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(); GaugeMat xform(grid); 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 GridBase *grid = Umu.Grid(); diff --git a/Grid/simd/Grid_sse4.h b/Grid/simd/Grid_sse4.h index eb76427e..2b2a80fd 100644 --- a/Grid/simd/Grid_sse4.h +++ b/Grid/simd/Grid_sse4.h @@ -35,7 +35,7 @@ Author: neo */ // Time-stamp: <2015-06-16 23:27:54 neo> //---------------------------------------------------------------------- - +#include #include NAMESPACE_BEGIN(Grid); diff --git a/Grid/simd/Simd.h b/Grid/simd/Simd.h index ddee6a36..cacb567f 100644 --- a/Grid/simd/Simd.h +++ b/Grid/simd/Simd.h @@ -258,6 +258,13 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexD &o){ stream<<">"; return stream; } +inline std::ostream& operator<< (std::ostream& stream, const vComplexD2 &o){ + stream<<"<"; + stream<"; + return stream; +} inline std::ostream& operator<< (std::ostream& stream, const vRealF &o){ int nn=vRealF::Nsimd(); diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc index 6de636bd..35ec2246 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc @@ -492,6 +492,7 @@ NAMESPACE_END(Grid); int main(int argc, char **argv) { +#if 0 Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); // 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; Grid_finalize(); return 0; +#endif } // main diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc index 42f54edd..004a0953 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc @@ -491,6 +491,7 @@ NAMESPACE_END(Grid); int main(int argc, char **argv) { +#if 0 Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); // 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::FermionField FermionFieldF; - typedef GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction MixedPrecRHMC; + typedef GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction MixedPrecRHMC; typedef GeneralEvenOddRatioRationalPseudoFermionAction DoublePrecRHMC; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -870,4 +871,5 @@ int main(int argc, char **argv) { std::cout << GridLogMessage << " Done" << std::endl; Grid_finalize(); return 0; +#endif } // main diff --git a/HMC/Mobius2p1f_DD_RHMC_96I.cc b/HMC/Mobius2p1f_DD_RHMC_96I.cc index 5158aed9..7ba5f5b5 100644 --- a/HMC/Mobius2p1f_DD_RHMC_96I.cc +++ b/HMC/Mobius2p1f_DD_RHMC_96I.cc @@ -286,8 +286,8 @@ int main(int argc, char **argv) { // ii) Break low bound, how rapidly? // iii) Run lanczos // iv) Have CG return spectral range estimate - FermionField vec(StrangeOp.FermionDedBlackGrid()); - FermionField res(StrangeOp.FermionDedBlackGrid()); + FermionField vec(StrangeOp.FermionRedBlackGrid()); + FermionField res(StrangeOp.FermionRedBlackGrid()); vec = 1; // Fill with any old junk std::cout << "Bounds check on strange operator mass "<< StrangeOp.Mass()< -Author: Guido Cossu - -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 - -NAMESPACE_BEGIN(Grid); - -template - class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction { - public: - typedef typename FermionOperatorD::FermionField FieldD; - typedef typename FermionOperatorF::FermionField FieldF; - - using OperatorFunction::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 " < &LinOpU, const FieldD &src, FieldD &psi) { - - std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<(&LinOpU); - - // std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <_Mat)<_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 "<Nd()<Nd()<(FermOpD.Umu, mu); - precisionChange(Umu_f,Umu_d); - PokeIndex(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" < 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 HMCWrapper; - // MD.name = std::string("Leap Frog"); - // typedef GenericHMCRunner HMCWrapper; - // MD.name = std::string("Force Gradient"); - typedef GenericHMCRunner 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 PlaqObs; - TheHMC.Resources.AddObservable(); - ////////////////////////////////////////////// - - 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 hasenbusch({ light_mass, 3.8e-3, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); - std::vector 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 LinearOperatorF; - typedef SchurDiagMooeeOperator LinearOperatorD; - typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; - - //////////////////////////////////////////////////////////////// - // Domain decomposed - //////////////////////////////////////////////////////////////// - Coordinate latt4 = GridPtr->GlobalDimensions(); - Coordinate mpi = GridPtr->ProcessorGrid(); - Coordinate shm; - - GlobalSharedMemory::GetShmDims(mpi,shm); - - Coordinate CommDim(Nd); - for(int d=0;d1 ? 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(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 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 CG(StoppingCondition,MaxCGIterations); - ConjugateGradient MDCG(MDStoppingCondition,MaxCGIterations); - - //////////////////////////////////// - // Collect actions - //////////////////////////////////// - ActionLevel Level1(1); - ActionLevel Level2(4); - ActionLevel 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 StrangePseudoFermionBdy(StrangeOpDir,StrangeOp, - StrangeOpDirF,StrangeOpF, - SFRp,500); - OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction StrangePseudoFermionLocal(StrangePauliVillarsOpDir,StrangeOpDir, - StrangePauliVillarsOpDirF,StrangeOpDirF, - SFRp,500); - OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction StrangePseudoFermionPVBdy(StrangePauliVillarsOp,StrangePauliVillarsOpDir, - StrangePauliVillarsOpF,StrangePauliVillarsOpDirF, - SFRp,500); -#else - OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionBdy(StrangeOpDir,StrangeOp,SFRp); - OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionLocal(StrangePauliVillarsOpDir,StrangeOpDir,SFRp); - OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermionPVBdy(StrangePauliVillarsOp,StrangePauliVillarsOpDir,SFRp); -#endif - Level1.push_back(&StrangePseudoFermionBdy); // ok - Level2.push_back(&StrangePseudoFermionLocal); - Level1.push_back(&StrangePseudoFermionPVBdy); //ok - - //////////////////////////////////// - // up down action - //////////////////////////////////// - std::vector light_den; - std::vector light_num; - std::vector dirichlet_den; - std::vector dirichlet_num; - - int n_hasenbusch = hasenbusch.size(); - light_den.push_back(light_mass); dirichlet_den.push_back(0); - for(int h=0;h Numerators; - std::vector NumeratorsF; - std::vector Denominators; - std::vector DenominatorsF; - std::vector *> Quotients; - -#define MIXED_PRECISION -#ifdef MIXED_PRECISION - std::vector *> Bdys; -#else - std::vector *> Bdys; -#endif - std::vector ActionMPCG; - std::vector MPCG; - - typedef SchurDiagMooeeOperator LinearOperatorF; - typedef SchurDiagMooeeOperator LinearOperatorD; - std::vector LinOpD; - std::vector LinOpF; - - for(int h=0;h(*Numerators[h],*Denominators[h],MDCG,CG)); - Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],CG)); - } else { -#ifdef MIXED_PRECISION - Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction( - *Numerators[h],*Denominators[h], - *NumeratorsF[h],*DenominatorsF[h], - OFRp, 500) ); - Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction( - *Numerators[h],*Denominators[h], - *NumeratorsF[h],*DenominatorsF[h], - OFRp, 500) ); -#else - Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction(*Numerators[h],*Denominators[h],OFRp)); - Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction(*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 seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); -#define DOUBLE +#define SINGLE #ifdef SINGLE typedef vComplexF Simd; typedef LatticeFermionF FermionField; diff --git a/benchmarks/Benchmark_dwf_sweep.cc b/benchmarks/Benchmark_dwf_sweep.cc index 2f13b667..2f11eb22 100644 --- a/benchmarks/Benchmark_dwf_sweep.cc +++ b/benchmarks/Benchmark_dwf_sweep.cc @@ -168,7 +168,7 @@ void benchDw(std::vector & latt4, int Ls, int threads,int report ) RealD M5 =1.8; 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(); Dw.Dhop(src,result,0); diff --git a/benchmarks/Benchmark_mooee.cc b/benchmarks/Benchmark_mooee.cc index 400199e2..22c1edbf 100644 --- a/benchmarks/Benchmark_mooee.cc +++ b/benchmarks/Benchmark_mooee.cc @@ -67,17 +67,17 @@ int main (int argc, char ** argv) const int ncall=1000; std::cout << GridLogMessage<< "*********************************************************" < diag = Dw.bs; Vector upper= Dw.cs; Vector lower= Dw.cs; diff --git a/benchmarks/Benchmark_staggered.cc b/benchmarks/Benchmark_staggered.cc index 34e1e470..a2be7f62 100644 --- a/benchmarks/Benchmark_staggered.cc +++ b/benchmarks/Benchmark_staggered.cc @@ -53,8 +53,8 @@ int main (int argc, char ** argv) pRNG.SeedFixedIntegers(seeds); // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); - typedef typename ImprovedStaggeredFermionR::FermionField FermionField; - typename ImprovedStaggeredFermionR::ImplParams params; + typedef typename ImprovedStaggeredFermionD::FermionField FermionField; + typename ImprovedStaggeredFermionD::ImplParams params; FermionField src (&Grid); random(pRNG,src); FermionField result(&Grid); result=Zero(); @@ -93,7 +93,7 @@ int main (int argc, char ** argv) RealD c1=9.0/8.0; RealD c2=-1.0/24.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<()); - WilsonFermionR Dw(Umu,Grid,RBGrid,mass,params); + WilsonFermionD Dw(Umu,Grid,RBGrid,mass,params); // Full operator bench_wilson(src,result,Dw,volume,DaggerNo); @@ -130,7 +130,7 @@ int main (int argc, char ** argv) void bench_wilson ( LatticeFermion & src, LatticeFermion & result, - WilsonFermionR & Dw, + WilsonFermionD & Dw, double const volume, int const dag ) { @@ -149,7 +149,7 @@ void bench_wilson ( void bench_wilson_eo ( LatticeFermion & src, LatticeFermion & result, - WilsonFermionR & Dw, + WilsonFermionD & Dw, double const volume, int const dag ) { diff --git a/examples/Example_Mobius_spectrum.cc b/examples/Example_Mobius_spectrum.cc index f4cd3335..b604eec4 100644 --- a/examples/Example_Mobius_spectrum.cc +++ b/examples/Example_Mobius_spectrum.cc @@ -253,7 +253,7 @@ int main (int argc, char ** argv) int nmass = masses.size(); - std::vector FermActs; + std::vector FermActs; std::cout< FermActs; + typedef MobiusFermionD FermionActionD; + std::vector FermActs; std::cout< boundary = {1,1,1,-1}; - FermionActionR::ImplParams Params(boundary); + FermionActionD::ImplParams Params(boundary); RealD b=1.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); diff --git a/examples/Example_taku1.cc b/examples/Example_taku1.cc index b3b01d6f..7c1b3526 100644 --- a/examples/Example_taku1.cc +++ b/examples/Example_taku1.cc @@ -416,10 +416,10 @@ int main (int argc, char ** argv) int nmass = masses.size(); - typedef DomainWallFermionR FermionActionR; - // typedef MobiusFermionR FermionActionR; - std::vector FermActs; - std::vector DWFActs; + typedef DomainWallFermionD FermionActionD; + // typedef MobiusFermionD FermionActionD; + std::vector FermActs; + std::vector DWFActs; std::cout< boundary = {1,1,1,-1}; - FermionActionR::ImplParams Params(boundary); + FermionActionD::ImplParams Params(boundary); RealD b=1.5; RealD c=0.5; std::cout< FermActs; + std::vector FermActs; std::cout< boundary = {1,1,1,-1}; - typedef MobiusFermionR FermionAction; + typedef MobiusFermionD FermionAction; FermionAction::ImplParams Params(boundary); for(int m=0;m FermActs; + std::vector FermActs; std::cout< boundary = {1,1,1,-1}; - typedef MobiusFermionR FermionAction; + typedef MobiusFermionD FermionAction; FermionAction::ImplParams Params(boundary); for(int m=0;m(Ddwf,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + TestWhat(Ddwf,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; @@ -97,54 +97,54 @@ int main (int argc, char ** argv) std::cout<(Dmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + MobiusFermionD Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + TestWhat(Dmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout< gamma(Ls,std::complex(1.0,0.0)); - ZMobiusFermionR zDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c); - TestWhat(zDmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + ZMobiusFermionD zDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c); + TestWhat(zDmob,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dzolo,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + TestWhat(Dzolo,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + TestWhat(Dsham,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dshamz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + TestWhat(Dshamz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dov,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + TestWhat(Dov,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dovz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); + TestWhat(Dovz,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5); std::cout< HermOp(Ddwf); + ZMobiusFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, Params.omega,1.,0.); + SchurDiagTwoOperator HermOp(Ddwf); // Eigenvector storage LanczosParams fine =Params.FineParams; diff --git a/tests/Test_gfield_shift.cc b/tests/Test_gfield_shift.cc new file mode 100644 index 00000000..cf450d3c --- /dev/null +++ b/tests/Test_gfield_shift.cc @@ -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 +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +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;muNd();mu++){ + Umu = PeekIndex(U, mu); + Umu = Gimpl::CshiftLink(Umu,dir,shift); + PokeIndex(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 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({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(link_field, lex, 0,0); //place on 0,0 element of link + + for(int mu=0;mu(U, link_field_2, mu); + } + } + + std::stringstream ss; + ss<<"error"; + for(int d=0;d 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<<" ["< + +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 + +using namespace Grid; + +//Linearly interpolate between two nearest times +RealD interpolate(const RealD t_int, const std::vector > &data){ + RealD tdiff1=1e32; int t1_idx=-1; + RealD tdiff2=1e32; int t2_idx=-1; + + for(int i=0;i seeds({1, 2, 3, 4, 5}); + GridSerialRNG sRNG; + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + + LatticeGaugeField U(&Grid); + SU::HotConfiguration(pRNG, U); + + int Nstep = 300; + RealD epsilon = 0.01; + RealD maxTau = Nstep*epsilon; + RealD tolerance = 1e-4; + + for(int i=1;i> tolerance; + } + } + std::cout << "Adaptive smear tolerance " << tolerance << std::endl; + + //Setup iterative Wilson flow + WilsonFlow wflow(epsilon,Nstep); + wflow.resetActions(); + + std::vector > 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 wflow_ad(epsilon,maxTau,tolerance); + wflow_ad.resetActions(); + + std::vector > 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