portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

First cut on generalised HMC

Browse Source 
Backward compatibility OK
This commit is contained in:
Guido Cossu 2016-10-03 15:28:00 +01:00
parent 7ea4b959a4
commit e415260961
7 changed files with 376 additions and 345 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
lib/qcd/action/gauge/GaugeImpl.h
View File

@ -38,15 +38,19 @@ namespace QCD {
template <class Gimpl> class WilsonLoops;
#define INHERIT_GIMPL_TYPES(GImpl) \
typedef typename GImpl::Simd Simd; \
typedef typename GImpl::GaugeLinkField GaugeLinkField; \
typedef typename GImpl::GaugeField GaugeField; \
typedef typename GImpl::SiteGaugeField SiteGaugeField; \
typedef typename GImpl::SiteGaugeLink SiteGaugeLink;
//
#define INHERIT_GIMPL_TYPES(GImpl) \
typedef typename GImpl::Simd Simd; \
typedef typename GImpl::LinkField GaugeLinkField; \
typedef typename GImpl::Field GaugeField; \
typedef typename GImpl::SiteField SiteGaugeField; \
typedef typename GImpl::SiteLink SiteGaugeLink;
//
template <class S, int Nrepresentation = Nc> class GaugeImplTypes {
#define INHERIT_FIELD_TYPES(Impl) \
typedef typename Impl::Field Field;
template <class S, int Nrepresentation = Nc > class GaugeImplTypes {
public:
typedef S Simd;
@ -55,16 +59,14 @@ public:
template <typename vtype>
using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation>>, Nd>;
typedef iImplGaugeLink<Simd> SiteGaugeLink;
typedef iImplGaugeField<Simd> SiteGaugeField;
typedef iImplGaugeLink<Simd> SiteLink;
typedef iImplGaugeField<Simd> SiteField;
typedef Lattice<SiteGaugeLink> GaugeLinkField; // bit ugly naming; polarised
// gauge field, lorentz... all
// ugly
typedef Lattice<SiteGaugeField> GaugeField;
typedef Lattice<SiteLink> LinkField;
typedef Lattice<SiteField> Field;
// Move this elsewhere? FIXME
static inline void AddGaugeLink(GaugeField &U, GaugeLinkField &W,
static inline void AddLink(Field &U, LinkField &W,
int mu) { // U[mu] += W
PARALLEL_FOR_LOOP
for (auto ss = 0; ss < U._grid->oSites(); ss++) {
@ -72,6 +74,28 @@ public:
U._odata[ss]._internal[mu] + W._odata[ss]._internal;
}
}
static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
// specific for SU gauge fields
LinkField Pmu(P._grid);
Pmu = zero;
for (int mu = 0; mu < Nd; mu++) {
SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
PokeIndex<LorentzIndex>(P, Pmu, mu);
}
}
static inline void update_field(Field& P, Field& U, double ep, unsigned int Nexp){
for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
auto Pmu = PeekIndex<LorentzIndex>(P, mu);
Umu = expMat(Pmu, ep, Nexp) * Umu;
PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
}
}
};
// Composition with smeared link, bc's etc.. probably need multiple inheritance

38
lib/qcd/hmc/HMC.h
View File

@ -60,24 +60,25 @@ struct HMCparameters {
/////////////////////////////////
}
void print() const {
std::cout << GridLogMessage << "[HMC parameter] Trajectories : " << Trajectories << "\n";
std::cout << GridLogMessage << "[HMC parameter] Start trajectory : " << StartTrajectory << "\n";
std::cout << GridLogMessage << "[HMC parameter] Metropolis test (on/off): " << MetropolisTest << "\n";
std::cout << GridLogMessage << "[HMC parameter] Thermalization trajs : " << NoMetropolisUntil << "\n";
void print_parameters() const {
std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n";
std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << MetropolisTest << "\n";
std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
}
};
template <class GaugeField>
template <class Field>
class HmcObservable {
public:
virtual void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
virtual void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
GridParallelRNG &pRNG) = 0;
};
// this is only defined for a gauge theory
template <class Gimpl>
class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> {
class PlaquetteLogger : public HmcObservable<typename Gimpl::Field> {
private:
std::string Stem;
@ -117,19 +118,19 @@ class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> {
}
};
// template <class GaugeField, class Integrator, class Smearer, class
// Boundary>
template <class GaugeField, class IntegratorType>
template <class IntegratorType>
class HybridMonteCarlo {
private:
const HMCparameters Params;
typedef typename IntegratorType::Field Field;
GridSerialRNG &sRNG; // Fixme: need a RNG management strategy.
GridParallelRNG &pRNG; // Fixme: need a RNG management strategy.
GaugeField &Ucur;
Field &Ucur;
IntegratorType &TheIntegrator;
std::vector<HmcObservable<GaugeField> *> Observables;
std::vector<HmcObservable<Field> *> Observables;
/////////////////////////////////////////////////////////
// Metropolis step
@ -164,7 +165,7 @@ class HybridMonteCarlo {
/////////////////////////////////////////////////////////
// Evolution
/////////////////////////////////////////////////////////
RealD evolve_step(GaugeField &U) {
RealD evolve_step(Field &U) {
TheIntegrator.refresh(U, pRNG); // set U and initialize P and phi's
RealD H0 = TheIntegrator.S(U); // initial state action
@ -191,20 +192,21 @@ class HybridMonteCarlo {
// Constructor
/////////////////////////////////////////
HybridMonteCarlo(HMCparameters Pams, IntegratorType &_Int,
GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U)
GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, Field &_U)
: Params(Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Ucur(_U) {}
~HybridMonteCarlo(){};
void AddObservable(HmcObservable<GaugeField> *obs) {
void AddObservable(HmcObservable<Field> *obs) {
Observables.push_back(obs);
}
void evolve(void) {
Real DeltaH;
GaugeField Ucopy(Ucur._grid);
Field Ucopy(Ucur._grid);
Params.print();
Params.print_parameters();
TheIntegrator.print_parameters();
// Actual updates (evolve a copy Ucopy then copy back eventually)
for (int traj = Params.StartTrajectory;

8
lib/qcd/hmc/HmcRunner.h
View File

@ -32,6 +32,7 @@ directory
namespace Grid {
namespace QCD {
// Class for HMC specific for gauge theories
template <class Gimpl, class RepresentationsPolicy = NoHirep >
class NerscHmcRunnerTemplate {
public:
@ -114,7 +115,7 @@ class NerscHmcRunnerTemplate {
*/
//////////////
NoSmearing<Gimpl> SmearingPolicy;
typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >
typedef MinimumNorm2<Gimpl, NoSmearing<Gimpl>, RepresentationsPolicy >
IntegratorType; // change here to change the algorithm
IntegratorParameters MDpar(20, 1.0);
IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
@ -157,10 +158,9 @@ class NerscHmcRunnerTemplate {
// smeared set
// notice that the unit configuration is singular in this procedure
std::cout << GridLogMessage << "Filling the smeared set\n";
SmearingPolicy.set_GaugeField(U);
SmearingPolicy.set_Field(U);
HybridMonteCarlo<GaugeField, IntegratorType> HMC(HMCpar, MDynamics, sRNG,
pRNG, U);
HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
HMC.AddObservable(&Checkpoint);
HMC.AddObservable(&PlaqLog);

2
lib/qcd/hmc/NerscCheckpointer.h
View File

@ -36,7 +36,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
namespace Grid{
namespace QCD{
// Only for gauge fields
template<class Gimpl>
class NerscHmcCheckpointer : public HmcObservable<typename Gimpl::GaugeField> {
private:

107
lib/qcd/hmc/integrators/Integrator.h
View File

@ -49,60 +49,53 @@ namespace Grid {
namespace QCD {
struct IntegratorParameters {
int Nexp;
int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize;
unsigned int
Nexp; // number of terms in the Taylor expansion of the exponential
unsigned int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize; // trajectory stepsize
IntegratorParameters(int MDsteps_, RealD trajL_ = 1.0, int Nexp_ = 12)
IntegratorParameters(int MDsteps_, RealD trajL_ = 1.0,
unsigned int Nexp_ = 12)
: Nexp(Nexp_),
MDsteps(MDsteps_),
trajL(trajL_),
stepsize(trajL / MDsteps){
// empty body constructor
};
};
void print_parameters() {
std::cout << GridLogMessage << "[Integrator] Trajectory length : " << trajL << std::endl;
std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl;
std::cout << GridLogMessage << "[Integrator] Step size : " << stepsize << std::endl;
std::cout << GridLogMessage << "[Integrator] Exponential approx.: " << Nexp << std::endl;
}
};
/*! @brief Class for Molecular Dynamics management */
template <class GaugeField, class SmearingPolicy, class RepresentationPolicy>
template <class FieldImplementation, class SmearingPolicy, class RepresentationPolicy>
class Integrator {
protected:
typedef IntegratorParameters ParameterType;
typedef typename FieldImplementation::Field MomentaField; //for readability
typedef typename FieldImplementation::Field Field;
IntegratorParameters Params;
const ActionSet<GaugeField, RepresentationPolicy> as;
const ActionSet<Field, RepresentationPolicy> as;
int levels; //
double t_U; // Track time passing on each level and for U and for P
std::vector<double> t_P; //
GaugeField P;
MomentaField P;
SmearingPolicy& Smearer;
RepresentationPolicy Representations;
// Should match any legal (SU(n)) gauge field
// Need to use this template to match Ncol to pass to SU<N> class
template <int Ncol, class vec>
void generate_momenta(Lattice<iVector<iScalar<iMatrix<vec, Ncol> >, Nd> >& P,
GridParallelRNG& pRNG) {
typedef Lattice<iScalar<iScalar<iMatrix<vec, Ncol> > > > GaugeLinkField;
GaugeLinkField Pmu(P._grid);
Pmu = zero;
for (int mu = 0; mu < Nd; mu++) {
SU<Ncol>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
PokeIndex<LorentzIndex>(P, Pmu, mu);
}
}
// ObserverList observers; // not yet
// typedef std::vector<Observer*> ObserverList;
// void register_observers();
// void notify_observers();
void update_P(GaugeField& U, int level, double ep) {
void update_P(Field& U, int level, double ep) {
t_P[level] += ep;
update_P(P, U, level, ep);
@ -129,12 +122,12 @@ class Integrator {
}
} update_P_hireps{};
void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) {
void update_P(MomentaField& Mom, Field& U, int level, double ep) {
// input U actually not used in the fundamental case
// Fundamental updates, include smearing
for (int a = 0; a < as[level].actions.size(); ++a) {
GaugeField force(U._grid);
GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
Field force(U._grid);
Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta
std::cout << GridLogIntegrator
@ -152,7 +145,7 @@ class Integrator {
as[level].apply(update_P_hireps, Representations, Mom, U, ep);
}
void update_U(GaugeField& U, double ep) {
void update_U(Field& U, double ep) {
update_U(P, U, ep);
t_U += ep;
@ -161,26 +154,21 @@ class Integrator {
<< "[" << fl << "] U "
<< " dt " << ep << " : t_U " << t_U << std::endl;
}
void update_U(GaugeField& Mom, GaugeField& U, double ep) {
// rewrite exponential to deal internally with the lorentz index?
for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu);
auto Pmu = PeekIndex<LorentzIndex>(Mom, mu);
Umu = expMat(Pmu, ep, Params.Nexp) * Umu;
PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
}
void update_U(MomentaField& Mom, Field& U, double ep) {
FieldImplementation::update_field(Mom, U, ep, Params.Nexp);
// Update the smeared fields, can be implemented as observer
Smearer.set_GaugeField(U);
Smearer.set_Field(U);
// Update the higher representations fields
Representations.update(U); // void functions if fundamental representation
}
virtual void step(GaugeField& U, int level, int first, int last) = 0;
virtual void step(Field& U, int level, int first, int last) = 0;
public:
Integrator(GridBase* grid, IntegratorParameters Par,
ActionSet<GaugeField, RepresentationPolicy>& Aset,
ActionSet<Field, RepresentationPolicy>& Aset,
SmearingPolicy& Sm)
: Params(Par),
as(Aset),
@ -195,6 +183,13 @@ class Integrator {
virtual ~Integrator() {}
virtual std::string integrator_name() = 0;
void print_parameters(){
std::cout << GridLogMessage << "[Integrator] Name : "<< integrator_name() << std::endl;
Params.print_parameters();
}
// to be used by the actionlevel class to iterate
// over the representations
struct _refresh {
@ -210,14 +205,14 @@ class Integrator {
} refresh_hireps{};
// Initialization of momenta and actions
void refresh(GaugeField& U, GridParallelRNG& pRNG) {
void refresh(Field& U, GridParallelRNG& pRNG) {
std::cout << GridLogIntegrator << "Integrator refresh\n";
generate_momenta(P, pRNG);
FieldImplementation::generate_momenta(P, pRNG);
// Update the smeared fields, can be implemented as observer
// necessary to keep the fields updated even after a reject
// of the Metropolis
Smearer.set_GaugeField(U);
Smearer.set_Field(U);
// Set the (eventual) representations gauge fields
Representations.update(U);
@ -228,7 +223,7 @@ class Integrator {
for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
// get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID
GaugeField& Us =
Field& Us =
Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
as[level].actions.at(actionID)->refresh(Us, pRNG);
}
@ -256,11 +251,12 @@ class Integrator {
} S_hireps{};
// Calculate action
RealD S(GaugeField& U) { // here also U not used
RealD S(Field& U) { // here also U not used
LatticeComplex Hloc(U._grid);
Hloc = zero;
// Momenta
// Momenta --- modify this (take the trace of field)
// momenta_action()
for (int mu = 0; mu < Nd; mu++) {
auto Pmu = PeekIndex<LorentzIndex>(P, mu);
Hloc -= trace(Pmu * Pmu);
@ -276,7 +272,7 @@ class Integrator {
for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
// get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID
GaugeField& Us =
Field& Us =
Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
Hterm = as[level].actions.at(actionID)->S(Us);
std::cout << GridLogMessage << "S Level " << level << " term "
@ -289,7 +285,7 @@ class Integrator {
return H;
}
void integrate(GaugeField& U) {
void integrate(Field& U) {
// reset the clocks
t_U = 0;
for (int level = 0; level < as.size(); ++level) {
@ -312,7 +308,12 @@ class Integrator {
// and that we indeed got to the end of the trajectory
assert(fabs(t_U - Params.trajL) < 1.0e-6);
}
};
}
}
#endif // INTEGRATOR_INCLUDED

508
lib/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -1,287 +1,291 @@
/*************************************************************************************
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/integrators/Integrator_algorithm.h
Source file: ./lib/qcd/hmc/integrators/Integrator_algorithm.h
Copyright (C) 2015
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp>
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 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.
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.
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 */
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
//--------------------------------------------------------------------
/*! @file Integrator_algorithm.h
* @brief Declaration of classes for the Molecular Dynamics algorithms
*
* @author Guido Cossu
*/
//--------------------------------------------------------------------
#ifndef INTEGRATOR_ALG_INCLUDED
#define INTEGRATOR_ALG_INCLUDED
namespace Grid{
namespace QCD{
namespace Grid {
namespace QCD {
/* PAB:
*
* Recursive leapfrog; explanation of nested stepping
*
* Nested 1:4; units in dt for top level integrator
*
* CHROMA IroIro
* 0 1 0
* P 1/2 P 1/2
* P 1/16 P1/16
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/8
* P 1 P 1
* P 1/16 * skipped --- avoids revaluating force
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/8
* P 1 P 1
* P 1/16 * skipped
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 * skipped
* P 1 P 1
* P 1/16 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/16
* P 1/2 P 1/2
*/
/* PAB:
*
* Recursive leapfrog; explanation of nested stepping
*
* Nested 1:4; units in dt for top level integrator
*
* CHROMA IroIro
* 0 1 0
* P 1/2 P 1/2
* P 1/16 P1/16
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/8
* P 1 P 1
* P 1/16 * skipped --- avoids revaluating force
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/8
* P 1 P 1
* P 1/16 * skipped
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 * skipped
* P 1 P 1
* P 1/16 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/8 P1/8
* U 1/8 U1/8
* P 1/16 P1/16
* P 1/2 P 1/2
*/
template<class GaugeField,
class SmearingPolicy,
class RepresentationPolicy = Representations< FundamentalRepresentation > > class LeapFrog :
public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
public:
template <class FieldImplementation, class SmearingPolicy,
class RepresentationPolicy =
Representations<FundamentalRepresentation> >
class LeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
RepresentationPolicy> {
public:
typedef LeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy>
Algorithm;
INHERIT_FIELD_TYPES(FieldImplementation);
typedef LeapFrog<GaugeField, SmearingPolicy, RepresentationPolicy> Algorithm;
std::string integrator_name(){return "LeapFrog";}
LeapFrog(GridBase* grid,
IntegratorParameters Par,
ActionSet<GaugeField, RepresentationPolicy> & Aset,
SmearingPolicy & Sm):
Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset,Sm) {};
LeapFrog(GridBase* grid, IntegratorParameters Par,
ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
: Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
grid, Par, Aset, Sm){};
void step(Field& U, int level, int _first, int _last) {
int fl = this->as.size() - 1;
// level : current level
// fl : final level
// eps : current step size
void step (GaugeField& U, int level,int _first, int _last){
// Get current level step size
RealD eps = this->Params.stepsize;
for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier;
int fl = this->as.size() -1;
// level : current level
// fl : final level
// eps : current step size
// Get current level step size
RealD eps = this->Params.stepsize;
for(int l=0; l<=level; ++l) eps/= this->as[l].multiplier;
int multiplier = this->as[level].multiplier;
for(int e=0; e<multiplier; ++e){
int multiplier = this->as[level].multiplier;
for (int e = 0; e < multiplier; ++e) {
int first_step = _first && (e == 0);
int last_step = _last && (e == multiplier - 1);
int first_step = _first && (e==0);
int last_step = _last && (e==multiplier-1);
if(first_step){ // initial half step
this->update_P(U, level,eps/2.0);
}
if(level == fl){ // lowest level
this->update_U(U, eps);
}else{ // recursive function call
this->step(U, level+1,first_step,last_step);
}
int mm = last_step ? 1 : 2;
this->update_P(U, level,mm*eps/2.0);
}
}
};
template<class GaugeField,
class SmearingPolicy,
class RepresentationPolicy = Representations < FundamentalRepresentation > > class MinimumNorm2 :
public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
private:
const RealD lambda = 0.1931833275037836;
public:
MinimumNorm2(GridBase* grid,
IntegratorParameters Par,
ActionSet<GaugeField, RepresentationPolicy> & Aset,
SmearingPolicy& Sm):
Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset,Sm) {};
void step (GaugeField& U, int level, int _first,int _last){
// level : current level
// fl : final level
// eps : current step size
int fl = this->as.size() -1;
RealD eps = this->Params.stepsize*2.0;
for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;
// Nesting: 2xupdate_U of size eps/2
// Next level is eps/2/multiplier
int multiplier = this->as[level].multiplier;
for(int e=0; e<multiplier; ++e){ // steps per step
int first_step = _first && (e==0);
int last_step = _last && (e==multiplier-1);
if(first_step){ // initial half step
this->update_P(U,level,lambda*eps);
}
if(level == fl){ // lowest level
this->update_U(U,0.5*eps);
}else{ // recursive function call
this->step(U,level+1,first_step,0);
}
this->update_P(U,level,(1.0-2.0*lambda)*eps);
if(level == fl){ // lowest level
this->update_U(U,0.5*eps);
}else{ // recursive function call
this->step(U,level+1,0,last_step);
}
int mm = (last_step) ? 1 : 2;
this->update_P(U,level,lambda*eps*mm);
}
}
};
template<class GaugeField,
class SmearingPolicy,
class RepresentationPolicy = Representations< FundamentalRepresentation > > class ForceGradient :
public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
private:
const RealD lambda = 1.0/6.0;;
const RealD chi = 1.0/72.0;
const RealD xi = 0.0;
const RealD theta = 0.0;
public:
// Looks like dH scales as dt^4. tested wilson/wilson 2 level.
ForceGradient(GridBase* grid,
IntegratorParameters Par,
ActionSet<GaugeField, RepresentationPolicy> & Aset,
SmearingPolicy &Sm):
Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset, Sm) {};
void FG_update_P(GaugeField&U, int level,double fg_dt,double ep){
GaugeField Ufg(U._grid);
GaugeField Pfg(U._grid);
Ufg = U;
Pfg = zero;
std::cout << GridLogMessage << "FG update "<<fg_dt<<" "<<ep<<std::endl;
// prepare_fg; no prediction/result cache for now
// could relax CG stopping conditions for the
// derivatives in the small step since the force gets multiplied by
// a tiny dt^2 term relative to main force.
//
// Presently 4 force evals, and should have 3, so 1.33x too expensive.
// could reduce this with sloppy CG to perhaps 1.15x too expensive
// even without prediction.
this->update_P(Pfg,Ufg,level,1.0);
this->update_U(Pfg,Ufg,fg_dt);
this->update_P(Ufg,level,ep);
if (first_step) { // initial half step
this->update_P(U, level, eps / 2.0);
}
void step (GaugeField& U, int level, int _first,int _last){
RealD eps = this->Params.stepsize*2.0;
for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;
RealD Chi = chi*eps*eps*eps;
int fl = this->as.size() -1;
int multiplier = this->as[level].multiplier;
for(int e=0; e<multiplier; ++e){ // steps per step
int first_step = _first && (e==0);
int last_step = _last && (e==multiplier-1);
if(first_step){ // initial half step
this->update_P(U,level,lambda*eps);
}
if(level == fl){ // lowest level
this->update_U(U,0.5*eps);
}else{ // recursive function call
this->step(U,level+1,first_step,0);
}
this->FG_update_P(U,level,2*Chi/((1.0-2.0*lambda)*eps),(1.0-2.0*lambda)*eps);
if(level == fl){ // lowest level
this->update_U(U,0.5*eps);
}else{ // recursive function call
this->step(U,level+1,0,last_step);
}
int mm = (last_step) ? 1 : 2;
this->update_P(U,level,lambda*eps*mm);
}
if (level == fl) { // lowest level
this->update_U(U, eps);
} else { // recursive function call
this->step(U, level + 1, first_step, last_step);
}
};
int mm = last_step ? 1 : 2;
this->update_P(U, level, mm * eps / 2.0);
}
}
};
template <class FieldImplementation, class SmearingPolicy,
class RepresentationPolicy =
Representations<FundamentalRepresentation> >
class MinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy,
RepresentationPolicy> {
private:
const RealD lambda = 0.1931833275037836;
public:
INHERIT_FIELD_TYPES(FieldImplementation);
MinimumNorm2(GridBase* grid, IntegratorParameters Par,
ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
: Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
grid, Par, Aset, Sm){};
std::string integrator_name(){return "MininumNorm2";}
void step(Field& U, int level, int _first, int _last) {
// level : current level
// fl : final level
// eps : current step size
int fl = this->as.size() - 1;
RealD eps = this->Params.stepsize * 2.0;
for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
// Nesting: 2xupdate_U of size eps/2
// Next level is eps/2/multiplier
int multiplier = this->as[level].multiplier;
for (int e = 0; e < multiplier; ++e) { // steps per step
int first_step = _first && (e == 0);
int last_step = _last && (e == multiplier - 1);
if (first_step) { // initial half step
this->update_P(U, level, lambda * eps);
}
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, first_step, 0);
}
this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, 0, last_step);
}
int mm = (last_step) ? 1 : 2;
this->update_P(U, level, lambda * eps * mm);
}
}
};
template <class FieldImplementation, class SmearingPolicy,
class RepresentationPolicy =
Representations<FundamentalRepresentation> >
class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy,
RepresentationPolicy> {
private:
const RealD lambda = 1.0 / 6.0;
;
const RealD chi = 1.0 / 72.0;
const RealD xi = 0.0;
const RealD theta = 0.0;
public:
INHERIT_FIELD_TYPES(FieldImplementation);
// Looks like dH scales as dt^4. tested wilson/wilson 2 level.
ForceGradient(GridBase* grid, IntegratorParameters Par,
ActionSet<Field, RepresentationPolicy>& Aset,
SmearingPolicy& Sm)
: Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
grid, Par, Aset, Sm){};
std::string integrator_name(){return "ForceGradient";}
void FG_update_P(Field& U, int level, double fg_dt, double ep) {
Field Ufg(U._grid);
Field Pfg(U._grid);
Ufg = U;
Pfg = zero;
std::cout << GridLogMessage << "FG update " << fg_dt << " " << ep
<< std::endl;
// prepare_fg; no prediction/result cache for now
// could relax CG stopping conditions for the
// derivatives in the small step since the force gets multiplied by
// a tiny dt^2 term relative to main force.
//
// Presently 4 force evals, and should have 3, so 1.33x too expensive.
// could reduce this with sloppy CG to perhaps 1.15x too expensive
// even without prediction.
this->update_P(Pfg, Ufg, level, 1.0);
this->update_U(Pfg, Ufg, fg_dt);
this->update_P(Ufg, level, ep);
}
void step(Field& U, int level, int _first, int _last) {
RealD eps = this->Params.stepsize * 2.0;
for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
RealD Chi = chi * eps * eps * eps;
int fl = this->as.size() - 1;
int multiplier = this->as[level].multiplier;
for (int e = 0; e < multiplier; ++e) { // steps per step
int first_step = _first && (e == 0);
int last_step = _last && (e == multiplier - 1);
if (first_step) { // initial half step
this->update_P(U, level, lambda * eps);
}
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, first_step, 0);
}
this->FG_update_P(U, level, 2 * Chi / ((1.0 - 2.0 * lambda) * eps),
(1.0 - 2.0 * lambda) * eps);
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, 0, last_step);
}
int mm = (last_step) ? 1 : 2;
this->update_P(U, level, lambda * eps * mm);
}
}
};
}
}
#endif//INTEGRATOR_INCLUDED
#endif // INTEGRATOR_INCLUDED

4
lib/qcd/smearing/GaugeConfiguration.h
View File

@ -21,7 +21,7 @@ public:
NoSmearing(): ThinLinks(NULL) {}
void set_GaugeField(GaugeField& U) { ThinLinks = &U; }
void set_Field(GaugeField& U) { ThinLinks = &U; }
void smeared_force(GaugeField& SigmaTilde) const {}
@ -227,7 +227,7 @@ class SmearedConfiguration {
// attach the smeared routines to the thin links U and fill the smeared set
void set_GaugeField(GaugeField& U) { fill_smearedSet(U); }
void set_Field(GaugeField& U) { fill_smearedSet(U); }
//====================================================================
void smeared_force(GaugeField& SigmaTilde) const {