1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-13 04:37:05 +01:00

removing Hadrons

This commit is contained in:
2020-03-27 12:03:32 +00:00
parent 05ebc458e2
commit 7e13724882
240 changed files with 1 additions and 31183 deletions

View File

@ -1,3 +0,0 @@
AM_LDFLAGS += -L../../Hadrons
include Make.inc

View File

@ -1,265 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_QED.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <v.m.guelpers@soton.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 <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<std::string> flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"};
std::vector<double> mass = {.2}; //{.01, .04, .2 , .25 , .3 };
unsigned int nt = GridDefaultLatt()[Tp];
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// pt source
MSource::Point::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point>("pt", ptPar);
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
std::string twist = "0. 0. 0. 0.";
//stochastic photon field
MGauge::StochEm::Par photonPar;
photonPar.gauge = PhotonR::Gauge::feynman;
photonPar.zmScheme = PhotonR::ZmScheme::qedL;
application.createModule<MGauge::StochEm>("ph_field", photonPar);
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 8;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.twist = "0. 0. 0. 0.";
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
// propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = "pt";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i],
quarkPar);
//seq sources with tadpole insertion
MSource::SeqConserved::Par seqPar_T;
seqPar_T.q = "Qpt_" + flavour[i] + "_5d";
seqPar_T.action = "DWF_" + flavour[i];
seqPar_T.tA = 0;
seqPar_T.tB = nt-1;
seqPar_T.curr_type = Current::Tadpole;
seqPar_T.mu_min = 0;
seqPar_T.mu_max = 3;
seqPar_T.mom = "0. 0. 0. 0.";
application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i]
+ "_seq_T", seqPar_T);
// seq propagator with tadpole insertion
MFermion::GaugeProp::Par quarkPar_seq_T;
quarkPar_seq_T.solver = "CG_" + flavour[i];
quarkPar_seq_T.source = "Qpt_" + flavour[i] + "_seq_T";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i]
+ "_seq_T" + flavour[i],
quarkPar_seq_T);
//seq sources with conserved vector and photon insertion
MSource::SeqConserved::Par seqPar_V;
seqPar_V.q = "Qpt_" + flavour[i] + "_5d";
seqPar_V.action = "DWF_" + flavour[i];
seqPar_V.tA = 0;
seqPar_V.tB = nt-1;
seqPar_V.curr_type = Current::Vector;
seqPar_V.mu_min = 0;
seqPar_V.mu_max = 3;
seqPar_V.mom = "0. 0. 0. 0.";
seqPar_V.photon = "ph_field";
application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i]
+ "_seq_V_ph", seqPar_V);
// seq propagator with conserved vector and photon insertion
MFermion::GaugeProp::Par quarkPar_seq_V;
quarkPar_seq_V.solver = "CG_" + flavour[i];
quarkPar_seq_V.source = "Qpt_" + flavour[i] + "_seq_V_ph";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i]
+ "_seq_V_ph_" + flavour[i],
quarkPar_seq_V);
//double seq sources with conserved vector and photon insertion
//(for self energy)
MSource::SeqConserved::Par seqPar_VV;
seqPar_VV.q = "Qpt_" + flavour[i] + "_seq_V_ph_"
+ flavour[i] + "_5d";
seqPar_VV.action = "DWF_" + flavour[i];
seqPar_VV.tA = 0;
seqPar_VV.tB = nt-1;
seqPar_VV.curr_type = Current::Vector;
seqPar_VV.mu_min = 0;
seqPar_VV.mu_max = 3;
seqPar_VV.mom = "0. 0. 0. 0.";
seqPar_VV.photon = "ph_field";
application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i]
+ "_seq_V_ph" + flavour[i]
+ "_seq_V_ph", seqPar_VV);
//double seq propagator with conserved vector and photon insertion
MFermion::GaugeProp::Par quarkPar_seq_VV;
quarkPar_seq_VV.solver = "CG_" + flavour[i];
quarkPar_seq_VV.source = "Qpt_" + flavour[i] + "_seq_V_ph"
+ flavour[i] + "_seq_V_ph";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i]
+ "_seq_V_ph_" + flavour[i]
+ "_seq_V_ph_" + flavour[i],
quarkPar_seq_VV);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
//2pt function contraction
MContraction::Meson::Par mesPar;
mesPar.output = "QED/pt_" + flavour[i] + flavour[j];
mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "(Gamma5 Gamma5)";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_pt_"
+ flavour[i] + flavour[j],
mesPar);
//tadpole contraction
MContraction::Meson::Par mesPar_seq_T;
mesPar_seq_T.output = "QED/tadpole_pt_" + flavour[i] + "_T_"
+ flavour[i] + "__" + flavour[j];
mesPar_seq_T.q1 = "Qpt_" + flavour[i] + "_seq_T" + flavour[i];
mesPar_seq_T.q2 = "Qpt_" + flavour[j];
mesPar_seq_T.gammas = "(Gamma5 Gamma5)";
mesPar_seq_T.sink = "sink";
application.createModule<MContraction::Meson>("meson_tadpole_pt_" +
flavour[i] + "_seq_T"
+ flavour[i] + flavour[j],
mesPar_seq_T);
//photon exchange contraction
MContraction::Meson::Par mesPar_seq_E;
mesPar_seq_E.output = "QED/exchange_pt_" + flavour[i] + "_V_ph_"
+ flavour[i] + "__" + flavour[j] + "_V_ph_"
+ flavour[j];
mesPar_seq_E.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i];
mesPar_seq_E.q2 = "Qpt_" + flavour[j] + "_seq_V_ph_" + flavour[j];
mesPar_seq_E.gammas = "(Gamma5 Gamma5)";
mesPar_seq_E.sink = "sink";
application.createModule<MContraction::Meson>("meson_exchange_pt_"
+ flavour[i] + "_seq_V_ph_" + flavour[i]
+ flavour[j] + "_seq_V_ph_" + flavour[j],
mesPar_seq_E);
//self energy contraction
MContraction::Meson::Par mesPar_seq_S;
mesPar_seq_S.output = "QED/selfenergy_pt_" + flavour[i] + "_V_ph_"
+ flavour[i] + "_V_ph_" + flavour[i] + "__"
+ flavour[j];
mesPar_seq_S.q1 = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i]
+ "_seq_V_ph_" + flavour[i];
mesPar_seq_S.q2 = "Qpt_" + flavour[j];
mesPar_seq_S.gammas = "(Gamma5 Gamma5)";
mesPar_seq_S.sink = "sink";
application.createModule<MContraction::Meson>("meson_selfenergy_pt_"
+ flavour[i] + "_seq_V_ph_"
+ flavour[i] + "_seq_V_ph_"
+ flavour[i] + flavour[j],
mesPar_seq_S);
}
// execution
application.saveParameterFile("QED.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,113 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_diskvector.cc
Copyright (C) 2015-2018
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 */
#define DV_DEBUG
#include <Hadrons/DiskVector.hpp>
using namespace Grid;
using namespace Grid::Hadrons;
GRID_SERIALIZABLE_ENUM(Enum, undef, red, 1, blue, 2, green, 3);
class Object: Serializable {
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Object,
Enum, e,
SpinColourMatrix, scm);
};
#ifdef HAVE_HDF5
typedef Hdf5Reader TestReader;
typedef Hdf5Writer TestWriter;
#else
typedef BinaryReader TestReader;
typedef BinaryWriter TestWriter;
#endif
int main(int argc, char *argv[])
{
Grid_init(&argc, &argv);
GridSerialRNG rng;
Object obj, v2w, v2r, v13w, v13r;
SerializableDiskVector<Object, TestReader, TestWriter> v("diskvector_test", 1000, 4);
obj.e = Enum::red;
random(rng, obj.scm);
v[32] = obj;
random(rng, obj.scm);
v[2] = obj;
v2w = obj;
random(rng, obj.scm);
v[6] = obj;
random(rng, obj.scm);
v[7] = obj;
random(rng, obj.scm);
v[8] = obj;
random(rng, obj.scm);
v[9] = obj;
random(rng, obj.scm);
v[10] = obj;
random(rng, obj.scm);
v[11] = obj;
random(rng, obj.scm);
v[12] = obj;
random(rng, obj.scm);
v[13] = obj;
v13w = obj;
random(rng, obj.scm);
v[14] = obj;
random(rng, obj.scm);
v[15] = obj;
v2r = v[2];
LOG(Message) << "v[2] correct? "
<< ((v2r == v2w) ? "yes" : "no" ) << std::endl;
v13r = v[13];
LOG(Message) << "v[13] correct? "
<< ((v13r == v13w) ? "yes" : "no" ) << std::endl;
LOG(Message) << "hit ratio " << v.hitRatio() << std::endl;
EigenDiskVector<ComplexD> w("eigendiskvector_test", 1000, 4);
EigenDiskVector<ComplexD>::Matrix m,n;
w[2] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
m = w[2];
w[3] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
w[4] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
w[5] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
w[6] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
w[7] = EigenDiskVectorMat<ComplexD>::Random(2000, 2000);
n = w[2];
LOG(Message) << "w[2] correct? "
<< ((m == n) ? "yes" : "no" ) << std::endl;
LOG(Message) << "hit ratio " << w.hitRatio() << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,382 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_distil.cc
Copyright (C) 2015-2019
Author: Felix Erben <ferben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@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 <typeinfo>
#include <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
// Very simple iterators for Eigen tensors
// The only way I could get these iterators to work is to put the begin() and end() functions in the Eigen namespace
// So if Eigen ever defines these, we'll have a conflict and have to change this
namespace Eigen {
template <typename ET>
inline typename std::enable_if<EigenIO::is_tensor<ET>::value, typename EigenIO::Traits<ET>::scalar_type *>::type
begin( ET & et ) { return reinterpret_cast<typename Grid::EigenIO::Traits<ET>::scalar_type *>(et.data()); }
template <typename ET>
inline typename std::enable_if<EigenIO::is_tensor<ET>::value, typename EigenIO::Traits<ET>::scalar_type *>::type
end( ET & et ) { return begin(et) + et.size() * EigenIO::Traits<ET>::count; }
}
/////////////////////////////////////////////////////////////
// Test creation of laplacian eigenvectors
/////////////////////////////////////////////////////////////
void test_Global(Application &application)
{
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1100;
globalPar.trajCounter.end = 1120;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
globalPar.graphFile = "";
globalPar.scheduleFile = "";
globalPar.saveSchedule = "false";
globalPar.parallelWriteMaxRetry = -1;
application.setPar(globalPar);
}
/////////////////////////////////////////////////////////////
// Create a random gauge with the correct name
/////////////////////////////////////////////////////////////
std::string test_Gauge(Application &application )
{
std::string sGaugeName{ "gauge" };
application.createModule<MGauge::Random>( sGaugeName );
return sGaugeName;
}
/////////////////////////////////////////////////////////////
// Test creation of laplacian eigenvectors
/////////////////////////////////////////////////////////////
void test_LapEvec(Application &application)
{
const char szModuleName[] = "LapEvec";
test_Gauge( application );
MDistil::LapEvecPar p;
p.gauge = "gauge";
p.Stout.steps = 3;
p.Stout.rho = 0.2;
p.Cheby.PolyOrder = 11;
p.Cheby.alpha = 0.55;
p.Cheby.beta = 35.5;
p.Lanczos.Nvec = 5;
p.Lanczos.Nk = 6;
p.Lanczos.Np = 2;
p.Lanczos.MaxIt = 1000;
p.Lanczos.resid = 1e-2;
p.Lanczos.IRLLog = 0;
application.createModule<MDistil::LapEvec>(szModuleName,p);
}
/////////////////////////////////////////////////////////////
// Test creation Solver
/////////////////////////////////////////////////////////////
std::string SolverName( const char * pSuffix = nullptr ) {
std::string sSolverName{ "CG" };
if( pSuffix && pSuffix[0] ) {
sSolverName.append( "_" );
sSolverName.append( pSuffix );
}
return sSolverName;
}
std::string test_Solver(Application &application, const char * pSuffix = nullptr )
{
std::string sActionName{ "DWF" };
if( pSuffix && pSuffix[0] ) {
sActionName.append( "_" );
sActionName.append( pSuffix );
}
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 16;
actionPar.M5 = 1.8;
actionPar.mass = 0.005;
actionPar.boundary = "1 1 1 -1";
actionPar.twist = "0. 0. 0. 0.";
application.createModule<MAction::DWF>( sActionName, actionPar );
MSolver::RBPrecCG::Par solverPar;
solverPar.action = sActionName;
solverPar.residual = 1.0e-2;
solverPar.maxIteration = 10000;
std::string sSolverName{ SolverName( pSuffix ) };
application.createModule<MSolver::RBPrecCG>( sSolverName, solverPar );
return sSolverName;
}
/////////////////////////////////////////////////////////////
// DistilParameters
/////////////////////////////////////////////////////////////
std::string test_DPar(Application &application) {
MDistil::DistilParameters DPar;
DPar.nvec = 5;
DPar.nnoise = 1;
DPar.tsrc = 0;
DPar.LI = 5;
DPar.TI = 8;
DPar.SI = 4;
std::string sDParName{"DPar_l"};
application.createModule<MDistil::DistilPar>(sDParName,DPar);
return sDParName;
}
/////////////////////////////////////////////////////////////
// Noises
/////////////////////////////////////////////////////////////
std::string test_Noises(Application &application, const std::string &sNoiseBaseName ) {
MDistil::NoisesPar NoisePar;
NoisePar.DistilParams = "DPar_l";
NoisePar.NoiseFileName = "noise";
std::string sNoiseName{"noise"};
application.createModule<MDistil::Noises>(sNoiseName,NoisePar);
return sNoiseName;
}
/////////////////////////////////////////////////////////////
// Perambulators
/////////////////////////////////////////////////////////////
std::string PerambulatorName( const char * pszSuffix = nullptr )
{
std::string sPerambulatorName{ "Peramb" };
if( pszSuffix && pszSuffix[0] )
sPerambulatorName.append( pszSuffix );
return sPerambulatorName;
}
void test_LoadPerambulators( Application &application, const char * pszSuffix = nullptr )
{
std::string sModuleName{ "Peramb_load" };
MIO::LoadPerambulator::Par PerambPar;
PerambPar.PerambFileName = "Peramb";
PerambPar.DistilParams = "DPar_l";
test_Noises(application, sModuleName); // I want these written after solver stuff
application.createModule<MIO::LoadPerambulator>( sModuleName, PerambPar );
}
void test_Perambulators( Application &application, const char * pszSuffix = nullptr )
{
std::string sModuleName{ PerambulatorName( pszSuffix ) };
// Perambulator parameters
MDistil::Perambulator::Par PerambPar;
PerambPar.lapevec = "LapEvec";
PerambPar.PerambFileName = sModuleName;
PerambPar.solver = test_Solver( application, pszSuffix );
PerambPar.DistilParams = "DPar_l";
PerambPar.noise = "noise";
test_Noises(application, sModuleName); // I want these written after solver stuff
application.createModule<MDistil::Perambulator>( sModuleName, PerambPar );
}
/////////////////////////////////////////////////////////////
// DistilVectors
/////////////////////////////////////////////////////////////
void test_DistilVectors(Application &application, const char * pszSuffix = nullptr, const char * pszNvec = nullptr )
{
std::string sModuleName{"DistilVecs"};
if( pszSuffix )
sModuleName.append( pszSuffix );
std::string sPerambName{"Peramb"};
if( pszSuffix )
sPerambName.append( pszSuffix );
MDistil::DistilVectors::Par DistilVecPar;
DistilVecPar.noise = "noise";
DistilVecPar.rho = "rho";
DistilVecPar.phi = "phi";
DistilVecPar.perambulator = sPerambName;
DistilVecPar.lapevec = "LapEvec";
DistilVecPar.DistilParams = "DPar_l";
application.createModule<MDistil::DistilVectors>(sModuleName,DistilVecPar);
}
/////////////////////////////////////////////////////////////
// MesonSink
/////////////////////////////////////////////////////////////
void test_MesonSink(Application &application)
{
// DistilVectors parameters
MContraction::A2AMesonField::Par A2AMesonFieldPar;
//A2AMesonFieldPar.left="Peramb_unsmeared_sink";
A2AMesonFieldPar.left="g5phi";
A2AMesonFieldPar.right="Peramb_unsmeared_sink";
A2AMesonFieldPar.output="DistilFields";
A2AMesonFieldPar.gammas="Identity";
A2AMesonFieldPar.mom={"0 0 0"};
A2AMesonFieldPar.cacheBlock=2;
A2AMesonFieldPar.block=4;
application.createModule<MContraction::A2AMesonField>("DistilMesonSink",A2AMesonFieldPar);
}
/////////////////////////////////////////////////////////////
// MesonFields
/////////////////////////////////////////////////////////////
void test_MesonField(Application &application, const char * pszFileSuffix,
const char * pszObjectLeft = nullptr, const char * pszObjectRight = nullptr )
{
// DistilVectors parameters
if( pszObjectLeft == nullptr )
pszObjectLeft = pszFileSuffix;
if( pszObjectRight == nullptr )
pszObjectRight = pszObjectLeft;
MContraction::A2AMesonField::Par A2AMesonFieldPar;
A2AMesonFieldPar.left="";
A2AMesonFieldPar.right=A2AMesonFieldPar.left;
A2AMesonFieldPar.left.append( pszObjectLeft );
A2AMesonFieldPar.right.append( pszObjectRight );
A2AMesonFieldPar.output="MesonSinks";
A2AMesonFieldPar.output.append( pszFileSuffix );
A2AMesonFieldPar.gammas="Identity";
A2AMesonFieldPar.mom={"0 0 0"};
A2AMesonFieldPar.cacheBlock=2;
A2AMesonFieldPar.block=4;
std::string sObjectName{"DistilMesonField"};
sObjectName.append( pszFileSuffix );
application.createModule<MContraction::A2AMesonField>(sObjectName, A2AMesonFieldPar);
}
bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
{
if( bGobbleWhiteSpace )
while( std::isspace(static_cast<unsigned char>(*pstr)) )
pstr++;
const char * p = pstr;
bool bMinus = false;
char c = * p++;
if( c == '+' )
c = * p++;
else if( c == '-' ) {
bMinus = true;
c = * p++;
}
int n = c - '0';
if( n < 0 || n > 9 )
return false;
while( * p >= '0' && * p <= '9' ) {
n = n * 10 + ( * p ) - '0';
p++;
}
if( bMinus )
n *= -1;
ri = n;
pstr = p;
return true;
}
int main(int argc, char *argv[])
{
// Decode command-line parameters. 1st one is which test to run
int iTestNum = -1;
for(int i = 1 ; i < argc ; i++ ) {
std::cout << "argv[" << i << "]=\"" << argv[i] << "\"" << std::endl;
const char * p = argv[i];
if( * p == '/' || * p == '-' ) {
p++;
char c = * p++;
switch(toupper(c)) {
case 'T':
if( bNumber( iTestNum, p ) ) {
std::cout << "Test " << iTestNum << " requested";
if( * p )
std::cout << " (ignoring trailer \"" << p << "\")";
std::cout << std::endl;
}
else
std::cout << "Invalid test \"" << &argv[i][2] << "\"" << std::endl;
break;
default:
std::cout << "Ignoring switch \"" << &argv[i][1] << "\"" << std::endl;
break;
}
}
}
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
// For now perform free propagator test - replace this with distillation test(s)
LOG(Message) << "====== Creating xml for test " << iTestNum << " ======" << std::endl;
switch(iTestNum) {
default: // 0
LOG(Message) << "Computing Meson 2pt-function" << std::endl;
test_Global( application );
test_LapEvec( application );
test_DPar( application );
test_Perambulators( application );
test_DistilVectors( application );
test_MesonField( application, "Phi", "phi" );
test_MesonField( application, "Rho", "rho" );
break;
case 1:
LOG(Message) << "Computing Meson 2pt-function by loading perambulators" << std::endl;
test_Global( application );
test_LapEvec( application );
test_DPar( application );
test_LoadPerambulators( application );
test_DistilVectors( application, "_load" );
test_MesonField( application, "Phi", "phi" );
test_MesonField( application, "Rho", "rho" );
break;
}
// execution
static const char XmlFileName[] = "test_distil.xml";
application.saveParameterFile( XmlFileName );
const Grid::Coordinate &lat{GridDefaultLatt()};
if( lat.size() == 4 && lat[0] == 4 && lat[1] == 4 && lat[2] == 4 && lat[3] == 8 )
application.run();
else
LOG(Warning) << "The parameters in " << XmlFileName << " are designed to run on a laptop usid --grid 4.4.4.8" << std::endl;
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,246 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_free_prop.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <v.m.guelpers@soton.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 <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<std::string> flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"};
std::vector<double> mass = {.2}; //{.01, .04, .2 , .25 , .3 };
std::vector<std::string> lepton_flavour = {"mu"};
std::vector<double> lepton_mass = {.2};
unsigned int nt = GridDefaultLatt()[Tp];
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// unit gauge field for lepton
application.createModule<MGauge::Unit>("free_gauge");
// pt source
MSource::Point::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point>("pt", ptPar);
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
//Propagators from FFT and Feynman rules
for (unsigned int i = 0; i < lepton_mass.size(); ++i)
{
//DWF actions
MAction::DWF::Par actionPar_lep;
actionPar_lep.gauge = "free_gauge";
actionPar_lep.Ls = 8;
actionPar_lep.M5 = 1.8;
actionPar_lep.mass = lepton_mass[i];
actionPar_lep.boundary = boundary;
application.createModule<MAction::DWF>("free_DWF_" + lepton_flavour[i], actionPar_lep);
//DWF free propagators
MFermion::FreeProp::Par freePar;
freePar.source = "pt";
freePar.action = "free_DWF_" + lepton_flavour[i];
freePar.twist = "0 0 0 0.5";
freePar.mass = lepton_mass[i];
application.createModule<MFermion::FreeProp>("Lpt_" + lepton_flavour[i],
freePar);
//Wilson actions
MAction::Wilson::Par actionPar_lep_W;
actionPar_lep_W.gauge = "free_gauge";
actionPar_lep_W.mass = lepton_mass[i];
actionPar_lep_W.boundary = boundary;
application.createModule<MAction::Wilson>("free_W_" + lepton_flavour[i], actionPar_lep_W);
//Wilson free propagators
MFermion::FreeProp::Par freePar_W;
freePar_W.source = "pt";
freePar_W.action = "free_W_" + lepton_flavour[i];
freePar_W.twist = "0 0 0 0.5";
freePar_W.mass = lepton_mass[i];
application.createModule<MFermion::FreeProp>("W_Lpt_" + lepton_flavour[i],
freePar_W);
}
//Propagators from inversion
for (unsigned int i = 0; i < flavour.size(); ++i)
{
//DWF actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 8;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
//DWF propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = "pt";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i],
quarkPar);
//Wilson actions
MAction::Wilson::Par actionPar_W;
actionPar_W.gauge = "gauge";
actionPar_W.mass = mass[i];
actionPar_W.boundary = boundary;
application.createModule<MAction::Wilson>("W_" + flavour[i], actionPar_W);
// solvers
MSolver::RBPrecCG::Par solverPar_W;
solverPar_W.action = "W_" + flavour[i];
solverPar_W.residual = 1.0e-8;
solverPar_W.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("W_CG_" + flavour[i],
solverPar_W);
//Wilson propagators
MFermion::GaugeProp::Par quarkPar_W;
quarkPar_W.solver = "W_CG_" + flavour[i];
quarkPar_W.source = "pt";
application.createModule<MFermion::GaugeProp>("W_Qpt_" + flavour[i],
quarkPar_W);
}
//2pt contraction for Propagators from FFT and Feynman rules
for (unsigned int i = 0; i < lepton_flavour.size(); ++i)
for (unsigned int j = i; j < lepton_flavour.size(); ++j)
{
//2pt function contraction DWF
MContraction::Meson::Par freemesPar;
freemesPar.output = "2pt_free/DWF_L_pt_" + lepton_flavour[i] + lepton_flavour[j];
freemesPar.q1 = "Lpt_" + lepton_flavour[i];
freemesPar.q2 = "Lpt_" + lepton_flavour[j];
freemesPar.gammas = "(Gamma5 Gamma5)";
freemesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_L_pt_"
+ lepton_flavour[i] + lepton_flavour[j],
freemesPar);
//2pt function contraction Wilson
MContraction::Meson::Par freemesPar_W;
freemesPar_W.output = "2pt_free/W_L_pt_" + lepton_flavour[i] + lepton_flavour[j];
freemesPar_W.q1 = "W_Lpt_" + lepton_flavour[i];
freemesPar_W.q2 = "W_Lpt_" + lepton_flavour[j];
freemesPar_W.gammas = "(Gamma5 Gamma5)";
freemesPar_W.sink = "sink";
application.createModule<MContraction::Meson>("W_meson_L_pt_"
+ lepton_flavour[i] + lepton_flavour[j],
freemesPar_W);
}
//2pt contraction for Propagators from inverion
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
//2pt function contraction DWF
MContraction::Meson::Par mesPar;
mesPar.output = "2pt_free/DWF_pt_" + flavour[i] + flavour[j];
mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "(Gamma5 Gamma5)";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_pt_"
+ flavour[i] + flavour[j],
mesPar);
//2pt function contraction Wilson
MContraction::Meson::Par mesPar_W;
mesPar_W.output = "2pt_free/W_pt_" + flavour[i] + flavour[j];
mesPar_W.q1 = "W_Qpt_" + flavour[i];
mesPar_W.q2 = "W_Qpt_" + flavour[j];
mesPar_W.gammas = "(Gamma5 Gamma5)";
mesPar_W.sink = "sink";
application.createModule<MContraction::Meson>("W_meson_pt_"
+ flavour[i] + flavour[j],
mesPar_W);
}
// execution
application.saveParameterFile("free_prop.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,188 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_meson_3pt.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<std::string> flavour = {"l", "s", "c1", "c2", "c3"};
std::vector<double> mass = {.01, .04, .2 , .25 , .3 };
unsigned int nt = GridDefaultLatt()[Tp];
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
globalPar.genetic.maxGen = 1000;
globalPar.genetic.maxCstGen = 200;
globalPar.genetic.popSize = 20;
globalPar.genetic.mutationRate = .1;
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
std::string twist = "0. 0. 0. 0.";
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 12;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.twist = twist;
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
}
for (unsigned int t = 0; t < nt; t += 1)
{
std::string srcName;
std::vector<std::string> qName;
std::vector<std::vector<std::string>> seqName;
// Z2 source
MSource::Z2::Par z2Par;
z2Par.tA = t;
z2Par.tB = t;
srcName = "z2_" + std::to_string(t);
application.createModule<MSource::Z2>(srcName, z2Par);
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// sequential sources
MSource::SeqGamma::Par seqPar;
qName.push_back("QZ2_" + flavour[i] + "_" + std::to_string(t));
seqPar.q = qName[i];
seqPar.tA = (t + nt/4) % nt;
seqPar.tB = (t + nt/4) % nt;
seqPar.mom = "1. 0. 0. 0.";
seqName.push_back(std::vector<std::string>(Nd));
for (unsigned int mu = 0; mu < Nd; ++mu)
{
seqPar.gamma = 0x1 << mu;
seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
+ "_" + std::to_string(seqPar.tA) + "-"
+ qName[i];
application.createModule<MSource::SeqGamma>(seqName[i][mu], seqPar);
}
// propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = srcName;
application.createModule<MFermion::GaugeProp>(qName[i], quarkPar);
for (unsigned int mu = 0; mu < Nd; ++mu)
{
quarkPar.source = seqName[i][mu];
seqName[i][mu] = "Q_" + flavour[i] + "-" + seqName[i][mu];
application.createModule<MFermion::GaugeProp>(seqName[i][mu], quarkPar);
}
}
// contractions
MContraction::Meson::Par mesPar;
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
mesPar.q1 = qName[i];
mesPar.q2 = qName[j];
mesPar.gammas = "all";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_Z2_"
+ std::to_string(t)
+ "_"
+ flavour[i]
+ flavour[j],
mesPar);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = 0; j < flavour.size(); ++j)
for (unsigned int mu = 0; mu < Nd; ++mu)
{
MContraction::Meson::Par mesPar;
mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_"
+ std::to_string(mu);
mesPar.q1 = qName[i];
mesPar.q2 = seqName[j][mu];
mesPar.gammas = "all";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("3pt_Z2_"
+ std::to_string(t)
+ "_"
+ flavour[i]
+ flavour[j]
+ "_"
+ std::to_string(mu),
mesPar);
}
}
// execution
application.saveParameterFile("meson3pt.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,155 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_spectrum.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<std::string> flavour = {"l", "s", "c1", "c2", "c3"};
std::vector<std::string> flavour_baryon = {"l", "s", "a", "b", "c"}; //needs to be a single character
std::vector<double> mass = {.01, .04, .2 , .25 , .3 };
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// sources
MSource::Z2::Par z2Par;
z2Par.tA = 0;
z2Par.tB = 0;
application.createModule<MSource::Z2>("z2", z2Par);
MSource::Point::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point>("pt", ptPar);
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
std::string twist = "0. 0. 0. 0.";
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 12;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.twist = twist;
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
// propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = "pt";
application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i], quarkPar);
quarkPar.source = "z2";
application.createModule<MFermion::GaugeProp>("QZ2_" + flavour[i], quarkPar);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
MContraction::Meson::Par mesPar;
mesPar.output = "mesons/pt_" + flavour[i] + flavour[j];
mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "all";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_pt_"
+ flavour[i] + flavour[j],
mesPar);
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
mesPar.q1 = "QZ2_" + flavour[i];
mesPar.q2 = "QZ2_" + flavour[j];
mesPar.gammas = "all";
mesPar.sink = "sink";
application.createModule<MContraction::Meson>("meson_Z2_"
+ flavour[i] + flavour[j],
mesPar);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
for (unsigned int k = j; k < flavour.size(); ++k)
{
MContraction::Baryon::Par barPar;
barPar.output = "baryons/pt_" + flavour[i] + flavour[j] + flavour[k];
barPar.q1 = "Qpt_" + flavour[i];
barPar.q2 = "Qpt_" + flavour[j];
barPar.q3 = "Qpt_" + flavour[k];
barPar.gammas = "(j12 j12) (j32X j32Y)";
barPar.quarks = flavour_baryon[i] + flavour_baryon[j] + flavour_baryon[k];
barPar.prefactors = "1.0";
barPar.sink = "sink";
application.createModule<MContraction::Baryon>(
"baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
}
// execution
application.saveParameterFile("spectrum.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

View File

@ -1,154 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_spectrum.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
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 <Hadrons/Application.hpp>
#include <Hadrons/Modules.hpp>
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<std::string> flavour = {"l", "s", "c"};
std::vector<double> mass = {.01, .04, .2 };
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
// sources
MSource::Point::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point>("pt_0", ptPar);
ptPar.position = "0 0 0 4";
application.createModule<MSource::Point>("pt_4", ptPar);
// sink
MSink::Point::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint>("sink", sinkPar);
application.createModule<MSink::Point>("sink_spec", sinkPar);
// set fermion boundary conditions to be periodic space, antiperiodic time.
std::string boundary = "1 1 1 -1";
std::string twist = "0. 0. 0. 0.";
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 12;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.twist = twist;
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
solverPar);
}
// propagators
MFermion::GaugeProp::Par quarkPar;
quarkPar.solver = "CG_l";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_l_0", quarkPar);
quarkPar.source = "pt_4";
application.createModule<MFermion::GaugeProp>("Qpt_l_4", quarkPar);
quarkPar.solver = "CG_s";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_s_0", quarkPar);
//This should be a loop - how do I make this?
quarkPar.solver = "CG_c";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_c_loop", quarkPar);
quarkPar.solver = "CG_l";
quarkPar.source = "pt_0";
application.createModule<MFermion::GaugeProp>("Qpt_l_loop", quarkPar);
MSink::Smear::Par smearPar;
smearPar.q="Qpt_l_0";
smearPar.sink = "sink_spec";
application.createModule<MSink::Smear>("Qpt_u_spec",smearPar);
MContraction::SigmaToNucleonEye::Par EyePar;
EyePar.output = "SigmaToNucleon/Eye_u";
EyePar.qqLoop = "Qpt_l_loop";
EyePar.quSpec = "Qpt_u_spec";
EyePar.qdTf = "Qpt_l_4";
EyePar.qsTi = "Qpt_s_0";
EyePar.tf = 4;
EyePar.sink = "sink";
application.createModule<MContraction::SigmaToNucleonEye>("SigmaToNucleonEye_u", EyePar);
EyePar.output = "SigmaToNucleon/Eye_c";
EyePar.qqLoop = "Qpt_c_loop";
application.createModule<MContraction::SigmaToNucleonEye>("SigmaToNucleonEye_c", EyePar);
MContraction::SigmaToNucleonNonEye::Par NonEyePar;
NonEyePar.output = "SigmaToNucleon/NonEye";
NonEyePar.quTi = "Qpt_l_0";
NonEyePar.quTf = "Qpt_l_4";
NonEyePar.quSpec = "Qpt_u_spec";
NonEyePar.qdTf = "Qpt_l_4";
NonEyePar.qsTi = "Qpt_s_0";
NonEyePar.tf = 4;
NonEyePar.sink = "sink";
application.createModule<MContraction::SigmaToNucleonNonEye>("SigmaToNucleonNonEye", NonEyePar);
// execution
application.saveParameterFile("stn.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}