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

added Production tests for MixedRep, Adj, 2S, 2AS. Still missing QObs. The HMC is not printing correctly all the actions and forces.

Browse Source
This commit is contained in:
pretidav
2017-11-04 18:16:54 +01:00
parent 915f610da0
commit a493429218
5 changed files with 491 additions and 106 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
lib/qcd/action/fermion/FermionCore.h
View File

@@ -70,7 +70,9 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
#define TwoIndexFermOpTemplateInstantiate(A) \
template class A<WilsonTwoIndexSymmetricImplF>; \
template class A<WilsonTwoIndexSymmetricImplD>;
template class A<WilsonTwoIndexSymmetricImplD>; \
template class A<WilsonTwoIndexAntiSymmetricImplF>; \
template class A<WilsonTwoIndexAntiSymmetricImplD>;
#define FermOp5dVecTemplateInstantiate(A) \
template class A<DomainWallVec5dImplF>; \

3
lib/qcd/action/fermion/WilsonKernelsHand.cc
View File

@@ -946,5 +946,6 @@ INSTANTIATE_THEM(DomainWallVec5dImplFH);
INSTANTIATE_THEM(DomainWallVec5dImplDF);
INSTANTIATE_THEM(ZDomainWallVec5dImplFH);
INSTANTIATE_THEM(ZDomainWallVec5dImplDF);
INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplF);
INSTANTIATE_THEM(WilsonTwoIndexAntiSymmetricImplD);
}}

162
tests/hmc/Test_hmc_WC2ASFG_Production.cc
View File

@@ -2,12 +2,11 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
Copyright (C) 2015
Copyright (C) 2017
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp>
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
@@ -27,103 +26,188 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "Grid/Grid.h"
#include <Grid/Grid.h>
int main(int argc, char **argv) {
namespace Grid{
struct FermionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
double, mass,
double, csw,
double, StoppingCondition,
int, MaxCGIterations,
bool, ApplySmearing);
};
struct WilsonCloverHMCParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
double, gauge_beta,
FermionParameters, WilsonClover)
template <class ReaderClass >
WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
struct SmearingParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
double, rho,
Integer, Nsmear)
template <class ReaderClass >
SmearingParameters(Reader<ReaderClass>& Reader){
read(Reader, "StoutSmearing", *this);
}
};
}
int main(int argc, char **argv)
{
using namespace Grid;
using namespace Grid::QCD;
// Here change the allowed (higher) representations
typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation > TheRepresentations;
typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation > TheRepresentations;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
// Typedefs to simplify notation
typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonTwoIndexAntiSymmetricImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions
typedef WilsonTwoIndexAntiSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
typedef WilsonCloverTwoIndexAntiSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
typedef typename FermionAction::FermionField FermionField;
typedef Grid::JSONReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Possibile to create the module by hand
// hardcoding parameters or using a Reader
TheHMC.ReadCommandLine(argc, argv);
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
WilsonCloverHMCParameters MyParams(Reader);
// Apply smearing to the fermionic action
bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams;
CheckpointerParameters CPparams(Reader);
/*
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
*/
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGModuleParameters RNGpar(Reader);
/*
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
*/
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
//typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
//TopologyObsParameters TopParams(Reader);
//TheHMC.Resources.AddObservable<QObs>(TopParams);
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 2.25 ;
WilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
//RealD beta = 5.6;
WilsonGaugeActionR Waction(MyParams.gauge_beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
TwoIndexSymmetricRepresentation::LatticeField U(GridPtr);
TwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr);
Real mass = -0.95;
//Real mass = 0.01;
//Real csw = 1.0;
// Can we define an overloaded operator that does not need U and initialises
// it with zeroes?
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass);
Real mass = MyParams.WilsonClover.mass;
Real csw = MyParams.WilsonClover.csw;
ConjugateGradient<FermionField> CG(1.0e-8, 2000, false);
std::cout << "mass and csw" << mass << " and " << csw << std::endl;
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
// Set smearing (true/false), default: false
Nf2.is_smeared = false;
Nf2.is_smeared = ApplySmearing;
// Collect actions
ActionLevel<LatticeGaugeField, TheRepresentations > Level1(1);
// Collect actions
ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<LatticeGaugeField, TheRepresentations > Level2(4);
ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run(); // no smearing
/*
double rho = 0.1; // smearing parameter
int Nsmear = 2; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
UGrid, Nsmear, Stout);
*/
// HMC parameters are serialisable
TheHMC.Parameters.initialize(Reader);
//TheHMC.Parameters.MD.MDsteps = 20;
//TheHMC.Parameters.MD.trajL = 1.0;
if (ApplySmearing){
SmearingParameters SmPar(Reader);
//double rho = 0.1; // smearing parameter
//int Nsmear = 3; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
TheHMC.Run(SmearingPolicy); // for smearing
} else {
TheHMC.Run(); // no smearing
}
//TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
//TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main

215
tests/hmc/Test_hmc_WCMixedRepFG_Production.cc
View File

@@ -32,6 +32,40 @@ directory
#include "Grid/Grid.h"
namespace Grid{
struct FermionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
double, mass,
double, csw,
double, StoppingCondition,
int, MaxCGIterations,
bool, ApplySmearing);
};
struct WilsonCloverHMCParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
double, gauge_beta,
FermionParameters, WilsonCloverFund,
FermionParameters, WilsonCloverAS)
template <class ReaderClass >
WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
struct SmearingParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
double, rho,
Integer, Nsmear)
template <class ReaderClass >
SmearingParameters(Reader<ReaderClass>& Reader){
read(Reader, "StoutSmearing", *this);
}
};
}
int main(int argc, char **argv) {
@@ -39,7 +73,7 @@ int main(int argc, char **argv) {
using namespace Grid::QCD;
// Here change the allowed (higher) representations
typedef Representations< FundamentalRepresentation, AdjointRepresentation , TwoIndexSymmetricRepresentation> TheRepresentations;
typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation> TheRepresentations;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
@@ -49,91 +83,142 @@ int main(int argc, char **argv) {
// Typedefs to simplify notation
typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
typedef WilsonAdjImplR AdjImplPolicy; // gauge field implemetation for the pseudofermions
typedef WilsonAdjFermionR AdjFermionAction; // type of lattice fermions (Wilson, DW, ...)
typedef WilsonTwoIndexSymmetricImplR SymmImplPolicy;
typedef WilsonTwoIndexSymmetricFermionR SymmFermionAction;
typedef WilsonImplR FundImplPolicy;
typedef WilsonCloverFermionR FundFermionAction;
typedef typename FundFermionAction::FermionField FundFermionField;
typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy;
typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction;
typedef typename ASymmFermionAction::FermionField ASymmFermionField;
typedef typename AdjFermionAction::FermionField AdjFermionField;
typedef typename SymmFermionAction::FermionField SymmFermionField;
typedef Grid::JSONReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Possibile to create the module by hand
// hardcoding parameters or using a Reader
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
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
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 2.25 ;
WilsonGaugeActionR Waction(beta);
// Grid from the command line
TheHMC.ReadCommandLine(argc, argv);
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
WilsonCloverHMCParameters MyParams(Reader);
// Apply smearing to the fermionic action
bool ApplySmearingFund = MyParams.WilsonCloverFund.ApplySmearing;
bool ApplySmearingAS = MyParams.WilsonCloverAS.ApplySmearing;
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
AdjointRepresentation::LatticeField UA(GridPtr);
TwoIndexSymmetricRepresentation::LatticeField US(GridPtr);
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams(Reader);
/*
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
*/
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar(Reader);
/*
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
*/
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
//typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
//TopologyObsParameters TopParams(Reader);
//TheHMC.Resources.AddObservable<QObs>(TopParams);
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
//RealD beta = 5.6;
WilsonGaugeActionR Waction(MyParams.gauge_beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
FundamentalRepresentation::LatticeField UF(GridPtr);
TwoIndexAntiSymmetricRepresentation::LatticeField UAS(GridPtr);
Real adjoint_mass = -0.1;
Real symm_mass = -0.5;
AdjFermionAction AdjFermOp(UA, *GridPtr, *GridRBPtr, adjoint_mass);
SymmFermionAction SymmFermOp(US, *GridPtr, *GridRBPtr, symm_mass);
ConjugateGradient<AdjFermionField> CG_adj(1.0e-8, 10000, false);
ConjugateGradient<SymmFermionField> CG_symm(1.0e-8, 10000, false);
Real Fundmass = MyParams.WilsonCloverFund.mass;
Real Fundcsw = MyParams.WilsonCloverFund.csw;
Real ASmass = MyParams.WilsonCloverAS.mass;
Real AScsw = MyParams.WilsonCloverAS.csw;
// Pass two solvers: one for the force computation and one for the action
TwoFlavourPseudoFermionAction<AdjImplPolicy> Nf2_Adj(AdjFermOp, CG_adj, CG_adj);
TwoFlavourPseudoFermionAction<SymmImplPolicy> Nf2_Symm(SymmFermOp, CG_symm, CG_symm);
std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl;
std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl;
FundFermionAction FundFermOp(UF, *GridPtr, *GridRBPtr, Fundmass, Fundcsw, Fundcsw);
ConjugateGradient<FundFermionField> CG_Fund(MyParams.WilsonCloverFund.StoppingCondition, MyParams.WilsonCloverFund.MaxCGIterations);
TwoFlavourPseudoFermionAction<FundImplPolicy> Nf2_Fund(FundFermOp, CG_Fund, CG_Fund);
ASymmFermionAction ASFermOp(UAS, *GridPtr, *GridRBPtr, ASmass, AScsw, AScsw);
ConjugateGradient<ASymmFermionField> CG_AS(MyParams.WilsonCloverAS.StoppingCondition, MyParams.WilsonCloverAS.MaxCGIterations);
TwoFlavourPseudoFermionAction<ASymmImplPolicy> Nf2_AS(ASFermOp, CG_AS, CG_AS);
Nf2_Fund.is_smeared = ApplySmearingFund;
Nf2_AS.is_smeared = ApplySmearingAS;
// Collect actions
ActionLevel<LatticeGaugeField, TheRepresentations > Level1(1);
Level1.push_back(&Nf2_Adj);
Level1.push_back(&Nf2_Symm);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1);
Level1.push_back(&Nf2_Fund);
Level1.push_back(&Nf2_AS);
ActionLevel<LatticeGaugeField, TheRepresentations > Level2(4);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
// HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.Parameters.initialize(Reader);
//TheHMC.Parameters.MD.MDsteps = 20;
//TheHMC.Parameters.MD.trajL = 1.0;
/*
if (ApplySmearingFund || ApplySmearingAS){
SmearingParameters SmPar(Reader);
//double rho = 0.1; // smearing parameter
//int Nsmear = 3; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
TheHMC.Run(SmearingPolicy); // for smearing
} else {
TheHMC.Run(); // no smearing
}
*/
TheHMC.Run();
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run(); // no smearing
//TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
//TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main

213
tests/hmc/Test_hmc_WCadjFG_Production.cc Normal file
View File

@@ -0,0 +1,213 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
namespace Grid{
struct FermionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
double, mass,
double, csw,
double, StoppingCondition,
int, MaxCGIterations,
bool, ApplySmearing);
};
struct WilsonCloverHMCParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
double, gauge_beta,
FermionParameters, WilsonClover)
template <class ReaderClass >
WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
struct SmearingParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
double, rho,
Integer, Nsmear)
template <class ReaderClass >
SmearingParameters(Reader<ReaderClass>& Reader){
read(Reader, "StoutSmearing", *this);
}
};
}
int main(int argc, char **argv)
{
using namespace Grid;
using namespace Grid::QCD;
typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonAdjImplR FermionImplPolicy; // gauge field implemetation for the pseudofermions
typedef WilsonCloverAdjFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
typedef typename FermionAction::FermionField FermionField;
typedef Grid::JSONReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.ReadCommandLine(argc, argv);
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
WilsonCloverHMCParameters MyParams(Reader);
// Apply smearing to the fermionic action
bool ApplySmearing = MyParams.WilsonClover.ApplySmearing;
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams(Reader);
/*
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
*/
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar(Reader);
/*
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
*/
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
//typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
//TopologyObsParameters TopParams(Reader);
//TheHMC.Resources.AddObservable<QObs>(TopParams);
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
//RealD beta = 5.6;
WilsonGaugeActionR Waction(MyParams.gauge_beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
AdjointRepresentation::LatticeField U(GridPtr);
//Real mass = 0.01;
//Real csw = 1.0;
Real mass = MyParams.WilsonClover.mass;
Real csw = MyParams.WilsonClover.csw;
std::cout << "mass and csw" << mass << " and " << csw << std::endl;
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw, csw);
ConjugateGradient<FermionField> CG(MyParams.WilsonClover.StoppingCondition, MyParams.WilsonClover.MaxCGIterations);
TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
// Set smearing (true/false), default: false
Nf2.is_smeared = ApplySmearing;
// Collect actions
ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
/*
double rho = 0.1; // smearing parameter
int Nsmear = 2; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
UGrid, Nsmear, Stout);
*/
// HMC parameters are serialisable
TheHMC.Parameters.initialize(Reader);
//TheHMC.Parameters.MD.MDsteps = 20;
//TheHMC.Parameters.MD.trajL = 1.0;
if (ApplySmearing){
SmearingParameters SmPar(Reader);
//double rho = 0.1; // smearing parameter
//int Nsmear = 3; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
TheHMC.Run(SmearingPolicy); // for smearing
} else {
TheHMC.Run(); // no smearing
}
//TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
//TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main