1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00

Adding a resource manager

This commit is contained in:
Guido Cossu
2016-12-22 12:41:56 +00:00
parent ce1a115e0b
commit 5214846341
10 changed files with 435 additions and 311 deletions

View File

@ -36,14 +36,13 @@ template <class GaugeField>
class Action {
public:
bool is_smeared = false;
// Boundary conditions? // Heatbath?
virtual void refresh(const GaugeField& U,
GridParallelRNG& pRNG) = 0; // refresh pseudofermions
virtual RealD S(const GaugeField& U) = 0; // evaluate the action
virtual void deriv(const GaugeField& U,
GaugeField& dSdU) = 0; // evaluate the action derivative
virtual std::string action_name() = 0; // return the action name
virtual ~Action(){};
// Heatbath?
virtual void refresh(const GaugeField& U, const GridParallelRNG& pRNG) = 0; // refresh pseudofermions
virtual RealD S(const GaugeField& U) = 0; // evaluate the action
virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0; // evaluate the action derivative
virtual std::string action_name() = 0; // return the action name
virtual std::string LogParameters() = 0; // prints action parameters
virtual ~Action(){}
};
// Indexing of tuple types
@ -60,32 +59,10 @@ struct Index<T, std::tuple<U, Types...>> {
static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
};
/*
template <class GaugeField>
struct ActionLevel {
public:
typedef Action<GaugeField>*
ActPtr; // now force the same colours as the rest of the code
//Add supported representations here
unsigned int multiplier;
std::vector<ActPtr> actions;
ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) {
assert(mul >= 1);
};
void push_back(ActPtr ptr) { actions.push_back(ptr); }
};
*/
template <class Field, class Repr = NoHirep >
struct ActionLevel {
public:
unsigned int multiplier;
unsigned int multiplier;
// Fundamental repr actions separated because of the smearing
typedef Action<Field>* ActPtr;
@ -98,15 +75,13 @@ struct ActionLevel {
std::vector<ActPtr>& actions;
ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
explicit ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
// initialize the hirep vectors to zero.
//apply(this->resize, actions_hirep, 0); //need a working resize
// apply(this->resize, actions_hirep, 0); //need a working resize
assert(mul >= 1);
};
//void push_back(ActPtr ptr) { actions.push_back(ptr); }
}
// void push_back(ActPtr ptr) { actions.push_back(ptr); }
template < class GenField >
void push_back(Action<GenField>* ptr) {

View File

@ -43,18 +43,22 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
public:
INHERIT_GIMPL_TYPES(Gimpl);
// typedef LorentzScalar<GaugeField> GaugeLinkField;
private:
RealD beta;
public:
WilsonGaugeAction(RealD b) : beta(b){};
explicit WilsonGaugeAction(RealD b) : beta(b){}
virtual std::string action_name(){return "WilsonGaugeAction";}
virtual std::string action_name() {return "WilsonGaugeAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "[WilsonGaugeAction] Beta: " << beta << std::endl;
return sstream.str();
}
virtual void refresh(const GaugeField &U,
GridParallelRNG &pRNG){}; // noop as no pseudoferms
const GridParallelRNG &pRNG){}; // noop as no pseudoferms
virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
@ -80,8 +84,7 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
}
};
}
};

View File

@ -30,176 +30,167 @@ with this program; if not, write to the Free Software Foundation, Inc.,
#ifndef GRID_GENERIC_HMC_RUNNER
#define GRID_GENERIC_HMC_RUNNER
#include <unordered_map>
namespace Grid {
namespace QCD {
// Virtual Class for HMC specific for gauge theories
// implement a specific theory by defining the BuildTheAction
template <class Implementation, class RepresentationsPolicy = NoHirep>
template <class Implementation,
template < typename, typename, typename > class Integrator,
class RepresentationsPolicy = NoHirep >
class BinaryHmcRunnerTemplate {
public:
INHERIT_FIELD_TYPES(Implementation);
typedef Implementation ImplPolicy;
INHERIT_FIELD_TYPES(Implementation);
typedef Implementation ImplPolicy; // visible from outside
template < typename S = NoSmearing<Implementation> >
using IntegratorType = Integrator<Implementation,S,RepresentationsPolicy>;
enum StartType_t { ColdStart,
HotStart,
TepidStart,
CheckpointStart };
enum StartType_t
{
ColdStart,
HotStart,
TepidStart,
CheckpointStart,
FilenameStart
};
ActionSet<Field, RepresentationsPolicy> TheAction;
struct HMCPayload
{
StartType_t StartType;
HMCparameters Parameters;
// A vector of HmcObservable
// that can be injected from outside
std::vector<HmcObservable<typename Implementation::Field> *>
ObservablesList;
HMCPayload() { StartType = HotStart; }
};
IntegratorParameters MDparameters;
// These can be rationalised, some private
HMCPayload Payload; // Parameters
HMCResourceManager Resources;
IntegratorParameters MDparameters;
GridCartesian * UGrid;
GridRedBlackCartesian *UrbGrid;
ActionSet<Field, RepresentationsPolicy> TheAction;
// A vector of HmcObservable that can be injected from outside
std::vector<HmcObservable<typename Implementation::Field> *> ObservablesList;
//GridCartesian * UGrid;
// These two are unnecessary, eliminate
GridCartesian * FGrid;
GridRedBlackCartesian *FrbGrid;
// GridRedBlackCartesian *UrbGrid;
// GridCartesian * FGrid;
// GridRedBlackCartesian *FrbGrid;
std::vector<int> SerialSeed;
std::vector<int> ParallelSeed;
void ReadCommandLine(int argc, char ** argv) {
std::string arg;
void RNGSeeds(std::vector<int> S, std::vector<int> P) {
SerialSeed = S;
ParallelSeed = P;
}
if (GridCmdOptionExists(argv, argv + argc, "--StartType")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--StartType");
if (arg == "HotStart") {
Payload.StartType = HotStart;
} else if (arg == "ColdStart") {
Payload.StartType = ColdStart;
} else if (arg == "TepidStart") {
Payload.StartType = TepidStart;
} else if (arg == "CheckpointStart") {
Payload.StartType = CheckpointStart;
} else {
std::cout << GridLogError << "Unrecognized option in --StartType\n";
std::cout
<< GridLogError
<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
assert(0);
}
}
virtual void BuildTheAction(int argc, char **argv) = 0; // necessary?
if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
Payload.Parameters.StartTrajectory = ivec[0];
}
// A couple of wrapper classes
template <class IOCheckpointer>
void Run(int argc, char **argv, IOCheckpointer &Checkpoint) {
NoSmearing<Implementation> S;
Runner(argc, argv, Checkpoint, S);
}
if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--Trajectories");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
Payload.Parameters.Trajectories = ivec[0];
}
template <class IOCheckpointer, class SmearingPolicy>
void Run(int argc, char **argv, IOCheckpointer &CP, SmearingPolicy &S) {
Runner(argc, argv, CP, S);
}
//////////////////////////////
if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
Payload.Parameters.NoMetropolisUntil = ivec[0];
}
}
// A couple of wrapper functions
template <class IOCheckpointer> void Run(IOCheckpointer &CP) {
NoSmearing<Implementation> S;
Runner(CP, S);
}
template <class SmearingPolicy, class IOCheckpointer>
void Runner(int argc,
char ** argv,
IOCheckpointer &Checkpoint,
SmearingPolicy &Smearing) {
StartType_t StartType = HotStart;
template <class IOCheckpointer, class SmearingPolicy> void Run(IOCheckpointer &CP, SmearingPolicy &S) {
Runner(CP, S);
}
std::string arg;
//////////////////////////////////////////////////////////////////
if (GridCmdOptionExists(argv, argv + argc, "--StartType")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--StartType");
if (arg == "HotStart") {
StartType = HotStart;
} else if (arg == "ColdStart") {
StartType = ColdStart;
} else if (arg == "TepidStart") {
StartType = TepidStart;
} else if (arg == "CheckpointStart") {
StartType = CheckpointStart;
} else {
std::cout << GridLogError << "Unrecognized option in --StartType\n";
std::cout
<< GridLogError
<< "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
assert(0);
}
}
private:
template <class SmearingPolicy, class IOCheckpointer>
void Runner(IOCheckpointer &Checkpoint, SmearingPolicy &Smearing) {
auto UGrid = Resources.GetCartesian();
Resources.AddRNGs();
Field U(UGrid);
int StartTraj = 0;
if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
StartTraj = ivec[0];
}
typedef IntegratorType<SmearingPolicy> TheIntegrator;
TheIntegrator MDynamics(UGrid, MDparameters, TheAction, Smearing);
int NumTraj = 1;
if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--Trajectories");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
NumTraj = ivec[0];
}
int NumThermalizations = 10;
if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
arg = GridCmdOptionPayload(argv, argv + argc, "--Thermalizations");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg, ivec);
NumThermalizations = ivec[0];
}
GridSerialRNG sRNG;
GridParallelRNG pRNG(UGrid);
Field U(UGrid);
typedef MinimumNorm2<Implementation, SmearingPolicy, RepresentationsPolicy> IntegratorType; // change here to change the algorithm
IntegratorType MDynamics(UGrid, MDparameters, TheAction, Smearing);
HMCparameters HMCpar;
HMCpar.StartTrajectory = StartTraj;
HMCpar.Trajectories = NumTraj;
HMCpar.NoMetropolisUntil = NumThermalizations;
if (StartType == HotStart) {
if (Payload.StartType == HotStart) {
// Hot start
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerialSeed);
pRNG.SeedFixedIntegers(ParallelSeed);
Implementation::HotConfiguration(pRNG, U);
} else if (StartType == ColdStart) {
Payload.Parameters.MetropolisTest = true;
Resources.SeedFixedIntegers();
Implementation::HotConfiguration(Resources.GetParallelRNG(), U);
} else if (Payload.StartType == ColdStart) {
// Cold start
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerialSeed);
pRNG.SeedFixedIntegers(ParallelSeed);
Implementation::ColdConfiguration(pRNG, U);
} else if (StartType == TepidStart) {
Payload.Parameters.MetropolisTest = true;
Resources.SeedFixedIntegers();
Implementation::ColdConfiguration(Resources.GetParallelRNG(), U);
} else if (Payload.StartType == TepidStart) {
// Tepid start
HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerialSeed);
pRNG.SeedFixedIntegers(ParallelSeed);
Implementation::TepidConfiguration(pRNG, U);
} else if (StartType == CheckpointStart) {
HMCpar.MetropolisTest = true;
Payload.Parameters.MetropolisTest = true;
Resources.SeedFixedIntegers();
Implementation::TepidConfiguration(Resources.GetParallelRNG(), U);
} else if (Payload.StartType == CheckpointStart) {
Payload.Parameters.MetropolisTest = true;
// CheckpointRestart
Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
}
Checkpoint.CheckpointRestore(Payload.Parameters.StartTrajectory, U, Resources.GetSerialRNG(), Resources.GetParallelRNG());
}
Smearing.set_Field(U);
Smearing.set_Field(U);
HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
HybridMonteCarlo<TheIntegrator> HMC(Payload.Parameters, MDynamics, Resources.GetSerialRNG(), Resources.GetParallelRNG(), U);
for (int obs = 0; obs < ObservablesList.size(); obs++)
HMC.AddObservable(ObservablesList[obs]);
for (int obs = 0; obs < ObservablesList.size(); obs++)
HMC.AddObservable(ObservablesList[obs]);
// Run it
HMC.evolve();
}
HMC.evolve();
}
};
// These are for gauge fields
typedef BinaryHmcRunnerTemplate<PeriodicGimplR> BinaryHmcRunner;
typedef BinaryHmcRunnerTemplate<PeriodicGimplF> BinaryHmcRunnerF;
typedef BinaryHmcRunnerTemplate<PeriodicGimplD> BinaryHmcRunnerD;
// These are for gauge fields, default integrator MinimumNorm2
template <template <typename, typename, typename> class Integrator > using BinaryHmcRunner = BinaryHmcRunnerTemplate<PeriodicGimplR, Integrator > ;
template <template <typename, typename, typename> class Integrator > using BinaryHmcRunnerF = BinaryHmcRunnerTemplate<PeriodicGimplF, Integrator > ;
template <template <typename, typename, typename> class Integrator > using BinaryHmcRunnerD = BinaryHmcRunnerTemplate<PeriodicGimplD, Integrator > ;
template <class RepresentationsPolicy>
using BinaryHmcRunnerTemplateHirep = BinaryHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
template <class RepresentationsPolicy, template <typename, typename, typename> class Integrator >
using BinaryHmcRunnerTemplateHirep = BinaryHmcRunnerTemplate<PeriodicGimplR, Integrator, RepresentationsPolicy>;
typedef BinaryHmcRunnerTemplate<ScalarImplR, ScalarFields>
ScalarBinaryHmcRunner;
typedef BinaryHmcRunnerTemplate<ScalarImplR, MinimumNorm2, ScalarFields> ScalarBinaryHmcRunner;
} // namespace QCD
} // namespace Grid
#endif
#endif // GRID_GENERIC_HMC_RUNNER

View File

@ -55,7 +55,7 @@ struct HMCparameters {
MetropolisTest = true;
NoMetropolisUntil = 10;
StartTrajectory = 0;
Trajectories = 200;
Trajectories = 10;
/////////////////////////////////
}
@ -206,10 +206,10 @@ class HybridMonteCarlo {
Params.print_parameters();
TheIntegrator.print_parameters();
TheIntegrator.print_actions();
// Actual updates (evolve a copy Ucopy then copy back eventually)
for (int traj = Params.StartTrajectory;
traj < Params.Trajectories + Params.StartTrajectory; ++traj) {
for (int traj = Params.StartTrajectory; traj < Params.Trajectories + Params.StartTrajectory; ++traj) {
std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
Ucopy = Ucur;

96
lib/qcd/hmc/HMCModules.h Normal file
View File

@ -0,0 +1,96 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/GenericHmcRunner.h
Copyright (C) 2015
Copyright (C) 2016
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 */
#ifndef GRID_HMC_MODULES
#define GRID_HMC_MODULES
namespace Grid {
namespace QCD {
///////////////////////////////////////////////////
// Modules
class GridModule {
public:
GridCartesian* get_full() { return grid_.get(); }
GridRedBlackCartesian* get_rb() { return rbgrid_.get(); }
void set_full(GridCartesian* grid) { grid_.reset(grid); }
void set_rb(GridRedBlackCartesian* rbgrid) { rbgrid_.reset(rbgrid); }
protected:
std::unique_ptr<GridCartesian> grid_;
std::unique_ptr<GridRedBlackCartesian> rbgrid_;
};
// helpers
class GridFourDimModule : public GridModule {
public:
GridFourDimModule() {
set_full(SpaceTimeGrid::makeFourDimGrid(
GridDefaultLatt(), GridDefaultSimd(4, vComplex::Nsimd()),
GridDefaultMpi()));
set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
}
};
class RNGModule{
// Random number generators
GridSerialRNG sRNG_;
std::unique_ptr<GridParallelRNG> pRNG_;
std::vector<int> SerialSeed_;
std::vector<int> ParallelSeed_;
public:
void set_pRNG(GridParallelRNG* pRNG){
pRNG_.reset(pRNG);
}
void set_RNGSeeds(const std::vector<int> S, const std::vector<int> P) {
SerialSeed_ = S;
ParallelSeed_ = P;
}
GridSerialRNG& get_sRNG(){return sRNG_;}
GridParallelRNG& get_pRNG(){return *pRNG_.get();}
void seed(){
sRNG_.SeedFixedIntegers(SerialSeed_);
pRNG_->SeedFixedIntegers(ParallelSeed_);
}
};
} // namespace QCD
} // namespace Grid
#endif // GRID_HMC_MODULES

View File

@ -0,0 +1,110 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/GenericHmcRunner.h
Copyright (C) 2015
Author: paboyle <paboyle@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 */
#ifndef HMC_RESOURCE_MANAGER_H
#define HMC_RESOURCE_MANAGER_H
#include <unordered_map>
namespace Grid {
namespace QCD {
// HMC Resource manager
class HMCResourceManager {
// Storage for grid pairs (std + red-black)
std::unordered_map<std::string, GridModule> Grids;
RNGModule RNGs;
bool have_RNG;
public:
HMCResourceManager():have_RNG(false){}
void AddGrid(std::string s, GridModule& M){
// Check for name clashes
auto search = Grids.find(s);
if(search != Grids.end()) {
std::cout << GridLogError << "Grid with name \"" << search->first << "\" already present. Terminating\n" ;
exit(1);
}
Grids[s] = std::move(M);
}
// Add a named grid set
void AddFourDimGrid(std::string s) {
GridFourDimModule Mod;
AddGrid(s,Mod);
}
GridCartesian* GetCartesian(std::string s="") {
if (s.empty()) s = Grids.begin()->first;
std::cout << GridLogDebug << "Getting cartesian grid from: "<< s << std::endl;
return Grids[s].get_full();
}
GridRedBlackCartesian* GetRBCartesian(std::string s="") {
if (s.empty()) s = Grids.begin()->first;
std::cout << GridLogDebug << "Getting rb-cartesian grid from: "<< s << std::endl;
return Grids[s].get_rb();
}
void AddRNGs(std::string s="") {
// Couple the RNGs to the GridModule tagged by s
// default is the first grid set
assert(Grids.size()>0 && !have_RNG );
if (s.empty()) s = Grids.begin()->first;
std::cout << GridLogDebug << "Adding RNG to grid: "<< s << std::endl;
RNGs.set_pRNG(new GridParallelRNG(GetCartesian(s)));
//pRNG.reset(new GridParallelRNG(GetCartesian(s)));
have_RNG = true;
}
void AddRNGSeeds(const std::vector<int> S, const std::vector<int> P) {
RNGs.set_RNGSeeds(S,P);
//SerialSeed = S;
//ParallelSeed = P;
}
GridSerialRNG& GetSerialRNG() {return RNGs.get_sRNG();}
GridParallelRNG& GetParallelRNG() {
assert(have_RNG);
return RNGs.get_pRNG();
}
void SeedFixedIntegers() {
assert(have_RNG);
RNGs.seed();
//sRNG.SeedFixedIntegers(SerialSeed);
//pRNG->SeedFixedIntegers(ParallelSeed);
}
};
}
}
#endif // HMC_RESOURCE_MANAGER_H

View File

@ -57,7 +57,7 @@ struct IntegratorParameters {
void print_parameters() {
std::cout << GridLogMessage << "[Integrator] Trajectory length : " << trajL << std::endl;
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;
}
@ -173,6 +173,20 @@ class Integrator {
Params.print_parameters();
}
void print_actions(){
std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
for (int level = 0; level < as.size(); ++level) {
std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
std::cout << as[level].actions.at(actionID)->LogParameters();
}
}
std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
}
// to be used by the actionlevel class to iterate
// over the representations
struct _refresh {
@ -189,8 +203,7 @@ class Integrator {
// Initialization of momenta and actions
void refresh(Field& U, GridParallelRNG& pRNG) {
//assert(P._grid == U._grid);
P.reset(U._grid);
assert(P._grid == U._grid);
std::cout << GridLogIntegrator << "Integrator refresh\n";
FieldImplementation::generate_momenta(P, pRNG);