portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-14 13:57:07 +01:00
Code Releases Activity

Merge branch 'develop' into feature/gpu-port

Browse Source
This commit is contained in:
Peter Boyle
2019-07-16 11:55:17 +01:00
parent 7c11525d1a c3d0c176ab
commit fa9cd50c5b
274 changed files with 7120 additions and 4663 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
tests/IO/Test_ildg_io.cc
View File

@ -53,7 +53,6 @@ int main (int argc, char ** argv)
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridCartesian Coarse(clatt_size,simd_layout,mpi_layout);
GridParallelRNG pRNGa(&Fine);
GridParallelRNG pRNGb(&Fine);
GridSerialRNG sRNGa;
@ -94,6 +93,27 @@ int main (int argc, char ** argv)
_IldgReader.close();
Umu_diff = Umu - Umu_saved;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
std::cout <<GridLogMessage<<"** Writing out ILDG conf *********"<<std::endl;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
file = std::string("./ckpoint_scidac.4000");
emptyUserRecord record;
ScidacWriter _ScidacWriter(Fine.IsBoss());
_ScidacWriter.open(file);
_ScidacWriter.writeScidacFieldRecord(Umu,record);
_ScidacWriter.close();
Umu_saved = Umu;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
std::cout <<GridLogMessage<<"** Reading back ILDG conf *********"<<std::endl;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
ScidacReader _ScidacReader;
_ScidacReader.open(file);
_ScidacReader.readScidacFieldRecord(Umu,record);
_ScidacReader.close();
Umu_diff = Umu - Umu_saved;
std::cout <<GridLogMessage<< "norm2 Gauge Diff = "<<norm2(Umu_diff)<<std::endl;
Grid_finalize();

193
tests/IO/Test_serialisation.cc
View File

@ -4,11 +4,12 @@
Source file: ./tests/Test_serialisation.cc
Copyright (C) 2015-2016
Copyright (C) 2015-2019
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.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
@ -28,6 +29,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <typeinfo>
using namespace Grid;
@ -80,26 +82,159 @@ double d = 2*M_PI;
bool b = false;
template <typename W, typename R, typename O>
void ioTest(const std::string &filename, const O &object, const std::string &name)
void ioTest(const std::string &filename, const O &object, const std::string &name,
const char * tag = "testobject", unsigned short Precision = 0 )
{
std::cout << "IO test: " << name << " -> " << filename << " ...";
// writer needs to be destroyed so that writing physically happens
{
W writer(filename);
write(writer, "testobject", object);
if( Precision )
writer.setPrecision(Precision);
write(writer, tag , object);
}
std::cout << " done. reading ...";
R reader(filename);
O buf;
bool good;
std::unique_ptr<O> buf( new O ); // In case object too big for stack
read(reader, "testobject", buf);
good = (object == buf);
std::cout << name << " IO test: " << (good ? "success" : "failure");
std::cout << std::endl;
if (!good) exit(EXIT_FAILURE);
read(reader, tag, *buf);
bool good = Serializable::CompareMember(object, *buf);
if (!good) {
std::cout << " failure!" << std::endl;
exit(EXIT_FAILURE);
}
std::cout << " done." << std::endl;
}
// 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; }
}
// Perform I/O tests on a range of tensor types
// Test coverage: scalars, complex and GridVectors in single, double and default precision
class TensorIO : public Serializable {
using TestScalar = ComplexD;
using SR3 = Eigen::Sizes<9,4,2>;
using SR5 = Eigen::Sizes<5,4,3,2,1>;
using ESO = Eigen::StorageOptions;
using TensorRank3 = Eigen::Tensor<ComplexF, 3, ESO::RowMajor>;
using TensorR5 = Eigen::TensorFixedSize<Real, SR5>;
using TensorR5Alt = Eigen::TensorFixedSize<Real, SR5, ESO::RowMajor>;
using Tensor942 = Eigen::TensorFixedSize<TestScalar, SR3, ESO::RowMajor>;
using aTensor942 = std::vector<Tensor942>;
using Perambulator = Eigen::Tensor<SpinColourVector, 6, ESO::RowMajor>;
using LSCTensor = Eigen::TensorFixedSize<SpinColourMatrix, Eigen::Sizes<6,5>>;
static const Real FlagR;
static const Complex Flag;
static const ComplexF FlagF;
static const TestScalar FlagTS;
static const char * const pszFilePrefix;
void Init(unsigned short Precision)
{
for( auto &s : Perambulator1 ) s = Flag;
for( auto &s : Perambulator2 ) s = Flag;
for( auto &s : tensorR5 ) s = FlagR;
for( auto &s : tensorRank3 ) s = FlagF;
for( auto &s : tensor_9_4_2 ) s = FlagTS;
for( auto &t : atensor_9_4_2 )
for( auto &s : t ) s = FlagTS;
for( auto &s : MyLSCTensor ) s = Flag;
}
// Perform an I/O test for a single Eigen tensor (of any type)
template <typename W, typename R, typename T, typename... IndexTypes>
static void TestOne(const char * MyTypeName, unsigned short Precision, std::string &filename,
const char * pszExtension, unsigned int &TestNum,
typename EigenIO::Traits<T>::scalar_type Flag, IndexTypes... otherDims)
{
using Traits = EigenIO::Traits<T>;
using scalar_type = typename Traits::scalar_type;
std::unique_ptr<T> pTensor{new T(otherDims...)};
for( auto &s : * pTensor ) s = Flag;
filename = pszFilePrefix + std::to_string(++TestNum) + "_" + MyTypeName + pszExtension;
ioTest<W, R, T>(filename, * pTensor, MyTypeName, MyTypeName);
}
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(TensorIO
, SpinColourVector, spinColourVector
, SpinColourMatrix, spinColourMatrix
, std::vector<std::string>, DistilParameterNames
, std::vector<int>, DistilParameterValues
, Perambulator, Perambulator1
, Perambulator, Perambulator2
, TensorR5, tensorR5
, TensorRank3, tensorRank3
, Tensor942, tensor_9_4_2
, aTensor942, atensor_9_4_2
, LSCTensor, MyLSCTensor
);
TensorIO()
: DistilParameterNames {"do", "androids", "dream", "of", "electric", "sheep?"}
, DistilParameterValues{2,3,1,4,5,1}
, Perambulator1(2,3,1,4,5,1)
, Perambulator2(7,1,6,1,5,1)
, tensorRank3(7,3,2)
, atensor_9_4_2(3) {}
#define TEST_PARAMS( T ) #T, Precision, filename, pszExtension, TestNum
// Perform a series of I/O tests for Eigen tensors, including a serialisable object
template <typename WTR_, typename RDR_>
static void Test(const char * pszExtension, unsigned short Precision = 0)
{
// Perform a series of tests on progressively more complex tensors
unsigned int TestNum = 0;
std::string filename;
// Rank 1 tensor containing a single integer
using TensorSingle = Eigen::TensorFixedSize<Integer, Eigen::Sizes<1>>;
TestOne<WTR_, RDR_, TensorSingle>( TEST_PARAMS( TensorSingle ), 7 ); // lucky!
// Rather convoluted way of defining four complex numbers
using TensorSimple = Eigen::Tensor<iMatrix<TestScalar,2>, 6>;
using I = typename TensorSimple::Index; // NB: Never specified, so same for all my test tensors
// Try progressively more complicated tensors
TestOne<WTR_, RDR_, TensorSimple, I,I,I,I,I,I>( TEST_PARAMS( TensorSimple ), FlagTS, 1,1,1,1,1,1 );
TestOne<WTR_, RDR_, TensorRank3, I, I, I>( TEST_PARAMS( TensorRank3 ), FlagF, 6, 3, 2 );
TestOne<WTR_, RDR_, Tensor942>(TEST_PARAMS( Tensor942 ), FlagTS);
TestOne<WTR_, RDR_, LSCTensor>(TEST_PARAMS( LSCTensor ), Flag );
TestOne<WTR_, RDR_, TensorR5>(TEST_PARAMS( TensorR5 ), FlagR);
// Now test a serialisable object containing a number of tensors
{
static const char MyTypeName[] = "TensorIO";
filename = pszFilePrefix + std::to_string(++TestNum) + "_" + MyTypeName + pszExtension;
std::unique_ptr<TensorIO> pObj{new TensorIO()};
pObj->Init(Precision);
ioTest<WTR_, RDR_, TensorIO>(filename, * pObj, MyTypeName, MyTypeName, Precision);
}
// Stress test. Too large for the XML or text readers and writers!
#ifdef STRESS_TEST
const std::type_info &tw = typeid( WTR_ );
if( tw == typeid( Hdf5Writer ) || tw == typeid( BinaryWriter ) ) {
using LCMTensor=Eigen::TensorFixedSize<iMatrix<iVector<iMatrix<iVector<LorentzColourMatrix,5>,2>,7>,3>,
Eigen::Sizes<2,4,11,10,9>, Eigen::StorageOptions::RowMajor>;
std::cout << "sizeof( LCMTensor ) = " << sizeof( LCMTensor ) / 1024 / 1024 << " MB" << std::endl;
TestOne<WTR_, RDR_, LCMTensor>(TEST_PARAMS( LCMTensor ), Flag);
}
#endif
}
};
const Real TensorIO::FlagR {1};
const Complex TensorIO::Flag {1,-1};
const ComplexF TensorIO::FlagF {1,-1};
const TensorIO::TestScalar TensorIO::FlagTS{1,-1};
const char * const TensorIO::pszFilePrefix = "tensor_";
template <typename T>
void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
{
@ -121,12 +256,12 @@ void tensorConvTestFn(GridSerialRNG &rng, const std::string label)
int main(int argc,char **argv)
{
Grid_init(&argc,&argv);
std::cout << std::boolalpha << "==== basic IO" << std::endl; // display true / false for boolean
GridSerialRNG rng;
rng.SeedFixedIntegers(std::vector<int>({42,10,81,9}));
std::cout << "==== basic IO" << std::endl;
XmlWriter WR("bother.xml");
// test basic type writing
@ -146,7 +281,6 @@ int main(int argc,char **argv)
// test serializable class writing
myclass obj(1234); // non-trivial constructor
std::vector<myclass> vec;
std::pair<myenum, myenum> pair;
std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
write(WR,"obj",obj);
@ -154,15 +288,15 @@ int main(int argc,char **argv)
vec.push_back(obj);
vec.push_back(myclass(5678));
vec.push_back(myclass(3838));
pair = std::make_pair(myenum::red, myenum::blue);
write(WR, "objvec", vec);
std::cout << "-- serialisable class writing to std::cout:" << std::endl;
std::cout << obj << std::endl;
std::cout << "-- serialisable class comparison:" << std::endl;
std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
std::cout << "vec[0] == obj: " << (vec[0] == obj) << std::endl;
std::cout << "vec[1] == obj: " << (vec[1] == obj) << std::endl;
std::cout << "-- pair writing to std::cout:" << std::endl;
std::pair<myenum, myenum> pair = std::make_pair(myenum::red, myenum::blue);
std::cout << pair << std::endl;
// read tests
@ -184,7 +318,15 @@ int main(int argc,char **argv)
#ifdef HAVE_HDF5
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5 (object) ");
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5 (vector of objects)");
std::cout << "\n==== detailed Hdf5 tensor tests (Grid::EigenIO)" << std::endl;
TensorIO::Test<Hdf5Writer, Hdf5Reader>(".h5");
#endif
std::cout << "\n==== detailed binary tensor tests (Grid::EigenIO)" << std::endl;
TensorIO::Test<BinaryWriter, BinaryReader>(".bin");
std::cout << "\n==== detailed xml tensor tests (Grid::EigenIO)" << std::endl;
TensorIO::Test<XmlWriter, XmlReader>(".xml", 6);
std::cout << "\n==== detailed text tensor tests (Grid::EigenIO)" << std::endl;
TensorIO::Test<TextWriter, TextReader>(".dat", 5);
std::cout << "\n==== vector flattening/reconstruction" << std::endl;
typedef std::vector<std::vector<std::vector<double>>> vec3d;
@ -227,4 +369,19 @@ int main(int argc,char **argv)
tensorConvTest(rng, ColourVector);
tensorConvTest(rng, SpinMatrix);
tensorConvTest(rng, SpinVector);
{
HMCparameters HMCparams;
HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.StartTrajectory =7;
HMCparams.Trajectories =1000;
HMCparams.NoMetropolisUntil=0;
HMCparams.MD.name =std::string("Force Gradient");
HMCparams.MD.MDsteps = 10;
HMCparams.MD.trajL = 1.0;
XmlWriter HMCwr("HMCparameters.xml");
write(HMCwr,"HMCparameters",HMCparams);
}
Grid_finalize();
}

40
tests/core/Test_fft_gfix.cc
View File

@ -60,11 +60,18 @@ int main (int argc, char ** argv)
std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
// int coulomb_dir = -1;
int coulomb_dir = Nd-1;
LatticeGaugeField Umu(&GRID);
LatticeGaugeField Urnd(&GRID);
LatticeGaugeField Uorg(&GRID);
LatticeGaugeField Utmp(&GRID);
LatticeColourMatrix g(&GRID); // Gauge xform
LatticeColourMatrix xform1(&GRID); // Gauge xform
LatticeColourMatrix xform2(&GRID); // Gauge xform
LatticeColourMatrix xform3(&GRID); // Gauge xform
SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
Uorg=Umu;
@ -78,7 +85,14 @@ int main (int argc, char ** argv)
Real alpha=0.1;
Umu = Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
// Check the gauge xform matrices
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform1);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
@ -92,7 +106,13 @@ int main (int argc, char ** argv)
std::cout<< "* Testing Fourier accelerated fixing *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
Umu=Urnd;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
Utmp=Urnd;
SU<Nc>::GaugeTransform(Utmp,xform2);
Utmp = Utmp - Umu;
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
@ -111,6 +131,22 @@ int main (int argc, char ** argv)
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
std::cout<< "*****************************************************************" <<std::endl;
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
Umu=Urnd;
SU3::HotConfiguration(pRNG,Umu); // Unit gauge
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
std::cout << Umu<<std::endl;
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Final plaquette "<<plaq << std::endl;
Grid_finalize();
}

12
tests/forces/Test_rect_force.cc
View File

@ -57,9 +57,10 @@ int main (int argc, char ** argv)
SU3::HotConfiguration(pRNG,U);
double beta = 1.0;
double c1 = 0.331;
double c1 = -0.331;
PlaqPlusRectangleActionR Action(beta,c1);
IwasakiGaugeActionR Action(beta);
// PlaqPlusRectangleActionR Action(beta,c1);
// WilsonGaugeActionR Action(beta);
ComplexD S = Action.S(U);
@ -91,6 +92,12 @@ int main (int argc, char ** argv)
auto mom_v = mom.View();
thread_foreach(i,mom_v,{ // exp(pmu dt) * Umu
Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt ;
Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt
+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
});
}
@ -117,6 +124,7 @@ int main (int argc, char ** argv)
}
ComplexD dSpred = sum(dS);
std::cout << std::setprecision(15)<<std::endl;
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;

7
tests/forces/Test_wilson_force.cc
View File

@ -170,6 +170,13 @@ int main (int argc, char ** argv)
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
// Smom = - P^2 ;
// dSmom = trace ( (mom+f/2dt)(mom+f/2dt) ) - trace mom*mom
// = trace(mom*f) dt + 0.25*dt*dt * trace(f*f).
//
// can we improve on this in HMC???
//
//
dSmom = dSmom - trace(mommu*forcemu) * dt;
dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt);

665
tests/hadrons/Test_hadrons.hpp
View File

@ -1,665 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons.hpp
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.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;
/*******************************************************************************
* Macros to reduce code duplication.
******************************************************************************/
// Common initialisation
#define HADRONS_DEFAULT_INIT \
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;
#define HADRONS_DEFAULT_GLOBALS(application) \
{ \
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); \
}
// Useful definitions
#define ZERO_MOM "0. 0. 0. 0."
#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n))
#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2)
#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3)
#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn)
#define LABEL_5D(s) s + "_5d";
// Wall source/sink macros
#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom)
#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM)
#define NAME_POINT_SOURCE(pos) ("point_" + pos)
// Meson module "gammas" special values
#define ALL_GAMMAS "all"
#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\
{\
std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\
makeWallSource(application, srcName, tW, mom);\
makePropagator(application, propName, srcName, solver);\
}
#define MAKE_WALL_PROP(tW, propName, solver)\
MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
// Sequential source macros
#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, seqPropName, solver, gamma)\
{\
std::string srcName = seqPropName + "_src";\
makeSequentialSource(application, srcName, qSrc, tS, gamma, mom);\
makePropagator(application, seqPropName, srcName, solver);\
}
// Point source macros
#define MAKE_POINT_PROP(pos, propName, solver)\
{\
std::string srcName = NAME_POINT_SOURCE(pos);\
makePointSource(application, srcName, pos);\
makePropagator(application, propName, srcName, solver);\
}
/*******************************************************************************
* Action setups.
******************************************************************************/
/*******************************************************************************
* Name: makeWilsonAction
* Parameters: application - main application that stores modules.
* actionName - name of action module to create.
* gaugeField - gauge field module.
* mass - quark mass.
* boundary - fermion boundary conditions (default to periodic
* space, antiperiodic time).
* Returns: None.
******************************************************************************/
inline void makeWilsonAction(Application &application, std::string actionName,
std::string &gaugeField, double mass,
std::string boundary = "1 1 1 -1")
{
if (!(VirtualMachine::getInstance().hasModule(actionName)))
{
MAction::Wilson::Par actionPar;
actionPar.gauge = gaugeField;
actionPar.mass = mass;
actionPar.boundary = boundary;
actionPar.twist = "0. 0. 0. 0.";
application.createModule<MAction::Wilson>(actionName, actionPar);
}
}
/*******************************************************************************
* Name: makeDWFAction
* Parameters: application - main application that stores modules.
* actionName - name of action module to create.
* gaugeField - gauge field module.
* mass - quark mass.
* M5 - domain wall height.
* Ls - fifth dimension extent.
* boundary - fermion boundary conditions (default to periodic
* space, antiperiodic time).
* Returns: None.
******************************************************************************/
inline void makeDWFAction(Application &application, std::string actionName,
std::string &gaugeField, double mass, double M5,
unsigned int Ls, std::string boundary = "1 1 1 -1")
{
if (!(VirtualMachine::getInstance().hasModule(actionName)))
{
MAction::DWF::Par actionPar;
actionPar.gauge = gaugeField;
actionPar.Ls = Ls;
actionPar.M5 = M5;
actionPar.mass = mass;
actionPar.boundary = boundary;
actionPar.twist = "0. 0. 0. 0.";
application.createModule<MAction::DWF>(actionName, actionPar);
}
}
/*******************************************************************************
* Functions for propagator construction.
******************************************************************************/
/*******************************************************************************
* Name: makeRBPrecCGSolver
* Purpose: Make RBPrecCG solver module for specified action.
* Parameters: application - main application that stores modules.
* solverName - name of solver module to create.
* actionName - action module corresponding to propagators to be
* computed.
* residual - CG target residual.
* Returns: None.
******************************************************************************/
inline void makeRBPrecCGSolver(Application &application, std::string &solverName,
std::string &actionName, double residual = 1e-8)
{
if (!(VirtualMachine::getInstance().hasModule(solverName)))
{
MSolver::RBPrecCG::Par solverPar;
solverPar.action = actionName;
solverPar.residual = residual;
solverPar.maxIteration = 10000;
application.createModule<MSolver::RBPrecCG>(solverName,
solverPar);
}
}
/*******************************************************************************
* Name: makePointSource
* Purpose: Construct point source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* pos - Position of point source.
* Returns: None.
******************************************************************************/
inline void makePointSource(Application &application, std::string srcName,
std::string pos)
{
// If the source already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Point::Par pointPar;
pointPar.position = pos;
application.createModule<MSource::Point>(srcName, pointPar);
}
}
/*******************************************************************************
* Name: makeSequentialSource
* Purpose: Construct sequential source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* qSrc - Input quark for sequential inversion.
* tS - sequential source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeSequentialSource(Application &application, std::string srcName,
std::string qSrc, unsigned int tS,
Gamma::Algebra gamma = Gamma::Algebra::GammaT,
std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::SeqGamma::Par seqPar;
seqPar.q = qSrc;
seqPar.tA = tS;
seqPar.tB = tS;
seqPar.mom = mom;
seqPar.gamma = gamma;
application.createModule<MSource::SeqGamma>(srcName, seqPar);
}
}
/*******************************************************************************
* Name: makeConservedSequentialSource
* Purpose: Construct sequential source with conserved current insertion and
* add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* qSrc - Input quark for sequential inversion.
* actionName - action corresponding to quark.
* tS - sequential source timeslice.
* curr - conserved current type to insert.
* mu - Lorentz index of current to insert.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeConservedSequentialSource(Application &application,
std::string &srcName,
std::string &qSrc,
std::string &actionName,
unsigned int tS,
Current curr,
unsigned int mu,
std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::SeqConserved::Par seqPar;
seqPar.q = qSrc;
seqPar.action = actionName;
seqPar.tA = tS;
seqPar.tB = tS;
seqPar.curr_type = curr;
seqPar.mu_min = mu;
seqPar.mu_min = mu;
seqPar.mom = mom;
application.createModule<MSource::SeqConserved>(srcName, seqPar);
}
}
/*******************************************************************************
* Name: makeNoiseSource
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* tA - lower source timeslice limit.
* tB - upper source timeslice limit.
* Returns: None.
******************************************************************************/
inline void makeNoiseSource(Application &application, std::string &srcName,
unsigned int tA, unsigned int tB)
{
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Z2::Par noisePar;
noisePar.tA = tA;
noisePar.tB = tB;
application.createModule<MSource::Z2>(srcName, noisePar);
}
}
/*******************************************************************************
* Name: makeWallSource
* Purpose: Construct wall source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* tW - wall source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeWallSource(Application &application, std::string &srcName,
unsigned int tW, std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(srcName)))
{
MSource::Wall::Par wallPar;
wallPar.tW = tW;
wallPar.mom = mom;
application.createModule<MSource::Wall>(srcName, wallPar);
}
}
/*******************************************************************************
* Name: makePointSink
* Purpose: Create function for point sink smearing of a propagator.
* Parameters: application - main application that stores modules.
* propName - name of input propagator.
* sinkFnct - name of output sink smearing module.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makePointSink(Application &application, std::string &sinkFnct,
std::string mom = ZERO_MOM)
{
// If the sink function already exists, don't make it again.
if (!(VirtualMachine::getInstance().hasModule(sinkFnct)))
{
MSink::Point::Par pointPar;
pointPar.mom = mom;
application.createModule<MSink::Point>(sinkFnct, pointPar);
}
}
/*******************************************************************************
* Name: sinkSmear
* Purpose: Perform sink smearing of a propagator.
* Parameters: application - main application that stores modules.
* sinkFnct - sink smearing module.
* propName - propagator to smear.
* smearedProp - name of output smeared propagator.
* Returns: None.
******************************************************************************/
inline void sinkSmear(Application &application, std::string &sinkFnct,
std::string &propName, std::string &smearedProp)
{
// If the propagator has already been smeared, don't smear it again.
if (!(VirtualMachine::getInstance().hasModule(smearedProp)))
{
MSink::Smear::Par smearPar;
smearPar.q = propName;
smearPar.sink = sinkFnct;
application.createModule<MSink::Smear>(smearedProp, smearPar);
}
}
/*******************************************************************************
* Name: makePropagator
* Purpose: Construct source and propagator then add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of source module to use.
* solver - solver to use (default is CG).
* Returns: None.
******************************************************************************/
inline void makePropagator(Application &application, std::string &propName,
std::string &srcName, std::string &solver)
{
// If the propagator already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(propName)))
{
MFermion::GaugeProp::Par quarkPar;
quarkPar.source = srcName;
quarkPar.solver = solver;
application.createModule<MFermion::GaugeProp>(propName, quarkPar);
}
}
/*******************************************************************************
* Name: makeLoop
* Purpose: Use noise source and inversion result to make loop propagator, then
* add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of noise source module to use.
* resName - name of inversion result on given noise source.
* Returns: None.
******************************************************************************/
inline void makeLoop(Application &application, std::string &propName,
std::string &srcName, std::string &resName)
{
// If the loop propagator already exists, don't make the module again.
if (!(VirtualMachine::getInstance().hasModule(propName)))
{
MLoop::NoiseLoop::Par loopPar;
loopPar.q = resName;
loopPar.eta = srcName;
application.createModule<MLoop::NoiseLoop>(propName, loopPar);
}
}
/*******************************************************************************
* Contraction module creation.
******************************************************************************/
/*******************************************************************************
* Name: mesonContraction
* Purpose: Create meson contraction module and add to application module.
* Parameters: application - main application that stores modules.
* modName - unique module name.
* output - name of output files.
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* sink - sink smearing module.
* gammas - gamma insertions at source and sink.
* Returns: None.
******************************************************************************/
inline void mesonContraction(Application &application,
std::string &modName, std::string &output,
std::string &q1, std::string &q2,
std::string &sink,
std::string gammas = "<Gamma5 Gamma5>")
{
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::Meson::Par mesPar;
mesPar.output = output;
mesPar.q1 = q1;
mesPar.q2 = q2;
mesPar.sink = sink;
mesPar.gammas = gammas;
application.createModule<MContraction::Meson>(modName, mesPar);
}
}
/*******************************************************************************
* Name: gamma3ptContraction
* Purpose: Create gamma3pt contraction module and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1, sink smeared.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* label - unique label to construct module name.
* tSnk - sink position of sink for q1.
* gamma - gamma insertions between q2 and q3.
* Returns: None.
******************************************************************************/
inline void gamma3ptContraction(Application &application, unsigned int npt,
std::string &q1, std::string &q2,
std::string &q3, std::string &label,
unsigned int tSnk = 0,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
std::string modName = std::to_string(npt) + "pt_" + label;
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::Gamma3pt::Par gamma3ptPar;
gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
gamma3ptPar.q1 = q1;
gamma3ptPar.q2 = q2;
gamma3ptPar.q3 = q3;
gamma3ptPar.tSnk = tSnk;
gamma3ptPar.gamma = gamma;
application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
}
}
/*******************************************************************************
* Name: weakContraction[Eye,NonEye]
* Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
* and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* q4 - quark propagator 4.
* label - unique label to construct module name.
* tSnk - time position of sink (for sink smearing).
* Returns: None.
******************************************************************************/
#define HW_CONTRACTION(top) \
inline void weakContraction##top(Application &application, unsigned int npt,\
std::string &q1, std::string &q2, \
std::string &q3, std::string &q4, \
std::string &label, unsigned int tSnk = 0)\
{\
std::string modName = std::to_string(npt) + "pt_" + label;\
if (!(VirtualMachine::getInstance().hasModule(modName)))\
{\
MContraction::WeakHamiltonian##top::Par weakPar;\
weakPar.output = std::to_string(npt) + "pt/" + label;\
weakPar.q1 = q1;\
weakPar.q2 = q2;\
weakPar.q3 = q3;\
weakPar.q4 = q4;\
weakPar.tSnk = tSnk;\
application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
}\
}
HW_CONTRACTION(Eye) // weakContractionEye
HW_CONTRACTION(NonEye) // weakContractionNonEye
/*******************************************************************************
* Name: disc0Contraction
* Purpose: Create contraction module for 4pt Weak Hamiltonian + current
* disconnected topology for neutral mesons and add to application
* module.
* Parameters: application - main application that stores modules.
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* q4 - quark propagator 4.
* label - unique label to construct module name.
* Returns: None.
******************************************************************************/
inline void disc0Contraction(Application &application,
std::string &q1, std::string &q2,
std::string &q3, std::string &q4,
std::string &label)
{
std::string modName = "4pt_" + label;
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::WeakNeutral4ptDisc::Par disc0Par;
disc0Par.output = "4pt/" + label;
disc0Par.q1 = q1;
disc0Par.q2 = q2;
disc0Par.q3 = q3;
disc0Par.q4 = q4;
application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
}
}
/*******************************************************************************
* Name: discLoopContraction
* Purpose: Create contraction module for disconnected loop and add to
* application module.
* Parameters: application - main application that stores modules.
* q_loop - loop quark propagator.
* modName - unique module name.
* gamma - gamma matrix to use in contraction.
* Returns: None.
******************************************************************************/
inline void discLoopContraction(Application &application,
std::string &q_loop, std::string &modName,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::DiscLoop::Par discPar;
discPar.output = "disc/" + modName;
discPar.q_loop = q_loop;
discPar.gamma = gamma;
application.createModule<MContraction::DiscLoop>(modName, discPar);
}
}
/*******************************************************************************
* Name: makeWITest
* Purpose: Create module to test Ward Identities for conserved current
* contractions and add to application module.
* Parameters: application - main application that stores modules.
* modName - name of module to create.
* propName - 4D quark propagator.
* actionName - action used to compute quark propagator.
* mass - mass of quark.
* Ls - length of 5th dimension (default = 1).
* test_axial - whether or not to check PCAC relation.
* Returns: None.
******************************************************************************/
inline void makeWITest(Application &application, std::string &modName,
std::string &propName, std::string &actionName,
double mass, unsigned int Ls = 1, bool test_axial = false)
{
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MContraction::WardIdentity::Par wiPar;
if (Ls > 1)
{
wiPar.q = LABEL_5D(propName);
}
else
{
wiPar.q = propName;
}
wiPar.action = actionName;
wiPar.mass = mass;
wiPar.test_axial = test_axial;
application.createModule<MContraction::WardIdentity>(modName, wiPar);
}
}
/*******************************************************************************
* Name: makeSeqCurrComparison
* Purpose: Create module to compare sequential insertion of conserved current
* against sink contraction and add to application module.
* Parameters: application - main application that stores modules.
* modName - name of module to create.
* propName - quark propagator (point source), 5D if available.
* seqName - 4D quark propagator with sequential insertion of
* conserved current.
* actionName - action used to compute quark propagators.
* origin - origin of point source propagator.
* t_J - time at which sequential current is inserted.
* mu - Lorentz index of sequential current.
* curr - type of conserved current inserted.
* Returns: None.
******************************************************************************/
inline void makeSeqCurrComparison(Application &application, std::string &modName,
std::string &propName, std::string &seqName,
std::string &actionName, std::string &origin,
unsigned int t_J, unsigned int mu, Current curr)
{
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MUtilities::TestSeqConserved::Par seqPar;
seqPar.q = propName;
seqPar.qSeq = seqName;
seqPar.action = actionName;
seqPar.origin = origin;
seqPar.t_J = t_J;
seqPar.mu = mu;
seqPar.curr = curr;
application.createModule<MUtilities::TestSeqConserved>(modName, seqPar);
}
}
/*******************************************************************************
* Name: makeSeqGamComparison
* Purpose: Create module to compare sequential insertion of gamma matrix
* against sink contraction and add to application module.
* Parameters: application - main application that stores modules.
* modName - name of module to create.
* propName - 4D quark propagator.
* seqProp - 4D quark propagator with sequential insertion of
* gamma matrix.
* gamma - Inserted gamma matrix.
* t_g - time at which gamma matrix is inserted
* sequentially.
* Returns: None.
******************************************************************************/
inline void makeSeqGamComparison(Application &application, std::string &modName,
std::string &propName, std::string &seqProp,
std::string &origin, Gamma::Algebra gamma,
unsigned int t_g)
{
if (!(VirtualMachine::getInstance().hasModule(modName)))
{
MUtilities::TestSeqGamma::Par seqPar;
seqPar.q = propName;
seqPar.qSeq = seqProp;
seqPar.origin = origin;
seqPar.t_g = t_g;
seqPar.gamma = gamma;
application.createModule<MUtilities::TestSeqGamma>(modName, seqPar);
}
}

123
tests/hadrons/Test_hadrons_3pt_contractions.cc
View File

@ -1,123 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_3pt_contractions.cc
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.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 "Test_hadrons.hpp"
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
HADRONS_DEFAULT_INIT;
// run setup ///////////////////////////////////////////////////////////////
Application application;
double mass = 0.04;
double M5 = 1.8;
unsigned int Ls = 12;
unsigned int nt = GridDefaultLatt()[Tp];
unsigned int t_i = 0;
unsigned int t_f = nt / 2;
std::string mom = "1. 0. 0. 0.";
// global parameters
HADRONS_DEFAULT_GLOBALS(application);
// gauge field
std::string gaugeField = "gauge";
application.createModule<MGauge::Unit>(gaugeField);
// Action & solver setup.
std::string action = "DWF";
std::string solver = "CG";
makeDWFAction(application, action, gaugeField, mass, M5, Ls);
makeRBPrecCGSolver(application, solver, action);
/***************************************************************************
* Weak Contraction test: Non-Eye class.
**************************************************************************/
// Make wall source propagators for each leg of 4-quark vertex.
std::string q_i_0 = "q_i_0";
std::string q_i_p = "q_i_p";
std::string q_f_0 = "q_f_0";
std::string q_f_p = "q_f_p";
MAKE_WALL_PROP(t_i, q_i_0, solver);
MAKE_WALL_PROP(t_f, q_f_0, solver);
MAKE_3MOM_WALL_PROP(t_i, mom, q_i_p, solver);
MAKE_3MOM_WALL_PROP(t_f, mom, q_f_p, solver);
// Perform contractions, zero and non-zero momentum.
std::string HW_CW_0 = LABEL_3PT("HW_CW_0", t_i, t_f);
std::string HW_CW_p = LABEL_3PT("HW_CW_p", t_i, t_f);
weakContractionNonEye(application, 3, q_i_0, q_i_0, q_f_0, q_f_0, HW_CW_0);
weakContractionNonEye(application, 3, q_i_0, q_i_p, q_f_p, q_f_0, HW_CW_p);
/***************************************************************************
* Weak Contraction test: Eye-class.
**************************************************************************/
// Create random propagator for loop.
std::string eta = "noise_source";
makeNoiseSource(application, eta, 0, nt - 1);
std::string loopProp = "loop";
std::string loopRes = loopProp + "_res";
makePropagator(application, loopRes, eta, solver);
makeLoop(application, loopProp, eta, loopRes);
// Wall sink smear the propagator directly connecting the source & sink.
// (i.e. make point sink but smear before the contraction)
std::string wallSink = "wall_sink";
std::string qWall = "q_wall";
makePointSink(application, wallSink);
sinkSmear(application, wallSink, q_i_0, qWall);
// Perform contractions, zero and non-zero momentum.
std::string HW_SE_0 = LABEL_3PT("HW_SE_0", t_i, t_f);
std::string HW_SE_p = LABEL_3PT("HW_SE_p", t_i, t_f);
weakContractionEye(application, 3, qWall, q_i_0, q_f_p, loopProp, HW_SE_0, t_f);
weakContractionEye(application, 3, qWall, q_i_p, q_f_p, loopProp, HW_SE_p, t_f);
/***************************************************************************
* Gamma insertion test.
**************************************************************************/
Gamma::Algebra gamma = Gamma::Algebra::GammaT;
std::string sd_0 = LABEL_3PT("sd_0", t_i, t_f);
std::string sd_p = LABEL_3PT("sd_p", t_i, t_f);
gamma3ptContraction(application, 3, qWall, q_i_0, q_f_0, sd_0, t_f, gamma);
gamma3ptContraction(application, 3, qWall, q_i_p, q_f_p, sd_p, t_f, gamma);
// execution
application.saveParameterFile("ContractionTest3pt.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

151
tests/hadrons/Test_hadrons_conserved_current.cc
View File

@ -1,151 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_conserved_current.cc
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.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 "Test_hadrons.hpp"
using namespace Grid;
using namespace Hadrons;
inline void setupSeqCurrTests(Application &application, std::string modStem,
std::string &pointProp, std::string &seqStem,
std::string &actionName, std::string &solverName,
std::string &origin, Current curr,
unsigned int t_J, unsigned int mu,
unsigned int Ls = 1)
{
std::string modName = ADD_INDEX(modStem, mu);
std::string seqProp = ADD_INDEX(seqStem, mu);
std::string seqSrc = seqProp + "_src";
// 5D actions require 5D propagator as input for conserved current
// insertions.
std::string propIn;
if (Ls > 1)
{
propIn = LABEL_5D(pointProp);
}
else
{
propIn = pointProp;
}
makeConservedSequentialSource(application, seqSrc, propIn,
actionName, t_J, curr, mu);
makePropagator(application, seqProp, seqSrc, solverName);
makeSeqCurrComparison(application, modName, propIn, seqProp,
actionName, origin, t_J, mu, curr);
}
inline void setupWardIdentityTests(Application &application,
std::string &actionName,
double mass,
unsigned int Ls = 1,
bool perform_axial_tests = false)
{
// solver
std::string solverName = actionName + "_CG";
makeRBPrecCGSolver(application, solverName, actionName);
unsigned int nt = GridDefaultLatt()[Tp];
unsigned int t_J = nt/2;
/***************************************************************************
* Conserved current sink contractions: use a single point propagator for
* the Ward Identity test.
**************************************************************************/
std::string pointProp = actionName + "_q_0";
std::string origin = "0 0 0 0";
std::string modName = actionName + " Ward Identity Test";
MAKE_POINT_PROP(origin, pointProp, solverName);
makeWITest(application, modName, pointProp, actionName, mass, Ls,
perform_axial_tests);
/***************************************************************************
* Conserved current tests with sequential insertion of vector/axial
* current. If above Ward Identity passes, sufficient to test sequential
* insertion of conserved current agrees with contracted version.
**************************************************************************/
// Compare sequential insertion to contraction. Should be enough to perform
// for time and one space component.
std::string seqStem = ADD_INDEX(pointProp + "seq_V", t_J);
std::string modStem = actionName + " Vector Sequential Test mu";
setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName,
solverName, origin, Current::Vector, t_J, Tp, Ls);
setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName,
solverName, origin, Current::Vector, t_J, Xp, Ls);
// Perform axial tests only if partially-conserved axial current exists for
// the action.
if (perform_axial_tests)
{
seqStem = ADD_INDEX(pointProp + "seq_A", t_J);
modStem = actionName + " Axial Sequential Test mu";
setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName,
solverName, origin, Current::Axial, t_J, Tp, Ls);
setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName,
solverName, origin, Current::Axial, t_J, Xp, Ls);
}
}
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
HADRONS_DEFAULT_INIT;
// run setup ///////////////////////////////////////////////////////////////
Application application;
double mass = 0.04;
double M5 = 1.8;
unsigned int Ls = 12;
// global parameters
HADRONS_DEFAULT_GLOBALS(application);
// gauge field
std::string gaugeField = "gauge";
application.createModule<MGauge::Unit>(gaugeField);
// Setup each action and the conserved current tests relevant to it.
std::string actionName = "DWF";
makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
setupWardIdentityTests(application, actionName, mass, Ls, true);
actionName = "Wilson";
makeWilsonAction(application, actionName, gaugeField, mass);
setupWardIdentityTests(application, actionName, mass);
// execution
application.saveParameterFile("ConservedCurrentTest.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

116
tests/hadrons/Test_hadrons_meson_conserved_3pt.cc
View File

@ -1,116 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_meson_conserved_3pt.cc
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.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 "Test_hadrons.hpp"
using namespace Grid;
using namespace Hadrons;
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
HADRONS_DEFAULT_INIT;
// run setup ///////////////////////////////////////////////////////////////
Application application;
// actions parameters
double mass = 0.04;
unsigned int Ls = 16;
double M5 = 1.8;
// kinematics
unsigned int nt = GridDefaultLatt()[Tp];
unsigned int tSrc = 0;
unsigned int tJ = nt / 4;
std::string kmom = "0. 0. 0. 0.";
std::string pmom = "1. 0. 0. 0.";
// Global parameters.
HADRONS_DEFAULT_GLOBALS(application);
// Unit gauge field.
std::string gaugeField = "Unit gauge";
application.createModule<MGauge::Unit>(gaugeField);
// DWF action
std::string actionName = "DWF";
makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
// Solver
std::string solver = "CG";
makeRBPrecCGSolver(application, solver, actionName);
// main test body //////////////////////////////////////////////////////////
// Point sink modules.
std::string sink_0 = "sink_0";
std::string sink_p = "sink_p";
MSink::Point::Par sinkPar;
sinkPar.mom = kmom;
application.createModule<MSink::ScalarPoint>(sink_0, sinkPar);
sinkPar.mom = pmom;
application.createModule<MSink::ScalarPoint>(sink_p, sinkPar);
// 2pt pion contraction, zero momentum.
std::string q_0 = "Q_0";
MAKE_WALL_PROP(tSrc, q_0, solver);
std::string modName = INIT_INDEX("2pt_pion_WP", tSrc);
std::string output = "2pt/pion_WP_0";
mesonContraction(application, modName, output, q_0, q_0, sink_0);
// 2pt pion contraction, with momentum p.
std::string q_p = "Q_p";
MAKE_3MOM_WALL_PROP(tSrc, pmom, q_p, solver);
modName = INIT_INDEX("2pt_pion_WP_p", tSrc);
output = "2pt/pion_WP_p";
mesonContraction(application, modName, output, q_0, q_p, sink_p);
// 3pt pion(0) -> pion(p), with sequentially inserted vector current in
// time direction.
std::string qSeq = q_0 + INIT_INDEX("_seq_Vc3", tJ);
std::string q5d = LABEL_5D(q_0); // Need 5D prop for DWF conserved current.
std::string srcName = qSeq + "_src";
modName = LABEL_3PT("3pt_pion_Vc3", tSrc, tJ);
output = "3pt/pion_Vc3_p";
makeConservedSequentialSource(application, srcName, q5d, actionName,
tJ, Current::Vector, Tp, pmom);
makePropagator(application, qSeq, srcName, solver);
mesonContraction(application, modName, output, q_0, qSeq, sink_p);
std::string par_file_name = "conserved_3pt.xml";
application.saveParameterFile(par_file_name);
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

90
tests/hadrons/Test_hadrons_seq_gamma.cc
View File

@ -1,90 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_seq_gamma.cc
Copyright (C) 2015-2018
Author: Andrew Lawson <andrew.lawson1991@gmail.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 "Test_hadrons.hpp"
using namespace Grid;
using namespace QCD;
using namespace Hadrons;
/*******************************************************************************
* Consistency test for sequential gamma insertion.
******************************************************************************/
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
HADRONS_DEFAULT_INIT;
// run setup ///////////////////////////////////////////////////////////////
Application application;
unsigned int nt = GridDefaultLatt()[Tp];
unsigned int tS = nt / 2;
unsigned int Ls = 12;
double mass = 0.04;
double M5 = 1.8;
// global parameters
HADRONS_DEFAULT_GLOBALS(application);
// gauge field
std::string gaugeField = "gauge";
application.createModule<MGauge::Unit>(gaugeField);
// action
std::string actionName = "DWF";
makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
// solver
std::string solverName = "CG";
makeRBPrecCGSolver(application, solverName, actionName);
// test sequential propagator, with g5 insertion.
Gamma::Algebra g = Gamma::Algebra::Gamma5;
std::string pointProp = "q_0";
std::string point5d = LABEL_5D(pointProp);
std::string origin = "0 0 0 0";
MAKE_POINT_PROP(origin, pointProp, solverName);
std::string seqProp = ADD_INDEX(pointProp + "_seqg5", tS);
std::string seqSrc = seqProp + "_src";
MAKE_SEQUENTIAL_PROP(tS, pointProp, ZERO_MOM, seqProp, solverName, g);
std::string modName = "Test g5 sequential insertion";
makeSeqGamComparison(application, modName, pointProp, seqProp, origin, g, tS);
// execution
application.saveParameterFile("SeqGamma5Test.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

159
tests/hadrons/Test_hadrons_wilsonFund.cc
View File

@ -1,159 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Tests/Hadrons/Test_hadrons_wilsonFund.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"};
std::vector<double> mass = {-0.1};
double csw = 0.0;
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 309;
globalPar.trajCounter.end = 310;
globalPar.trajCounter.step = 1;
globalPar.runId = "test";
application.setPar(globalPar);
// gauge field
application.createModule<MIO::LoadNersc>("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";
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::WilsonClover::Par actionPar;
actionPar.gauge = "gauge";
actionPar.mass = mass[i];
actionPar.boundary = boundary;
actionPar.csw_r = csw;
actionPar.csw_t = csw;
// !!!!! Check if Anisotropy works !!!!!
actionPar.clover_anisotropy.isAnisotropic= false;
actionPar.clover_anisotropy.t_direction = 3 ; // Explicit for D=4
actionPar.clover_anisotropy.xi_0 = 1.0 ;
actionPar.clover_anisotropy.nu = 1.0 ;
application.createModule<MAction::WilsonClover>("WilsonClover_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "WilsonClover_" + 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 = "Fund_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 = "Fund_baryons/pt_" + flavour[i] + flavour[j] + flavour[k];
barPar.q1 = "Qpt_" + flavour[i];
barPar.q2 = "Qpt_" + flavour[j];
barPar.q3 = "Qpt_" + flavour[k];
application.createModule<MContraction::Baryon>(
"baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
}
// execution
application.saveParameterFile("WilsonClover_spectrum.xml");
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

217
tests/hmc/Test_hmc_Mobius2p1f.cc Normal file
View File

@ -0,0 +1,217 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_EODWFRatio.cc
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@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 */
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
using namespace Grid::QCD;
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 WilsonImplR FermionImplPolicy;
typedef MobiusFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
typedef Grid::XmlReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD;
typedef GenericHMCRunner<LeapFrog> HMCWrapper;
MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient");
//typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
// MD.name = std::string("MinimumNorm2");
MD.MDsteps = 40;
MD.trajL = 1.0;
HMCparameters HMCparams;
HMCparams.StartTrajectory = 0;
HMCparams.Trajectories = 1;
HMCparams.NoMetropolisUntil= 20;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
HMCparams.StartingType =std::string("ColdStart");
HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams);
// Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 10;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
const int Ls = 4;
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.04;
Real pv_mass = 1.0;
RealD M5 = 1.8;
RealD b = 1.5; // Scale factor two
RealD c = 0.5;
// RHMC
// OneFlavourRationalParams OFRp;
// OFRp.lo = 1.0e-2;
// OFRp.hi = 25;
// OFRp.MaxIter = 10000;
// OFRp.tolerance= 1.0e-7;
// OFRp.degree = 10;
// OFRp.precision= 40;
// EOFA
OneFlavourRationalParams OFRp;
OFRp.lo = 0.98;
OFRp.hi = 25.0;
OFRp.MaxIter = 10000;
OFRp.tolerance= 1.0e-7;
OFRp.degree = 10;
OFRp.precision= 40;
std::vector<Real> hasenbusch({ 0.1 });
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
IwasakiGaugeActionR GaugeAction(beta);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
double StoppingCondition = 1e-10;
double MaxCGIterations = 30000;
ConjugateGradient<LatticeFermion> CG(StoppingCondition,MaxCGIterations);
////////////////////////////////////
// Collect actions
////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(4);
////////////////////////////////////
// Strange action
////////////////////////////////////
// Setup for RHMC
// FermionAction StrangeOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_mass,M5,b,c, Params);
// OneFlavourRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermion(StrangeOp,OFRp);
// Level1.push_back(&StrangePseudoFermion);
// DJM: setup for EOFA ratio (Shamir)
// DomainWallEOFAFermionR Strange_Op_L(U, *FGrid, *FrbGrid, *GridPtr, *GridRBPtr, strange_mass, strange_mass, pv_mass, 0.0, -1, M5);
// DomainWallEOFAFermionR Strange_Op_R(U, *FGrid, *FrbGrid, *GridPtr, *GridRBPtr, pv_mass, strange_mass, pv_mass, -1.0, 1, M5);
// ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> EOFA(Strange_Op_L, Strange_Op_R, CG, OFRp, true);
// Level1.push_back(&EOFA);
// DJM: setup for EOFA ratio (Mobius)
MobiusEOFAFermionR Strange_Op_L(U, *FGrid, *FrbGrid, *GridPtr, *GridRBPtr, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Strange_Op_R(U, *FGrid, *FrbGrid, *GridPtr, *GridRBPtr, pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> EOFA(Strange_Op_L, Strange_Op_R, CG, OFRp, true);
Level1.push_back(&EOFA);
// Setup for RHMC ratio
// FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params);
// FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params);
// OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermion(StrangePauliVillarsOp,StrangeOp,OFRp);
// Level1.push_back(&StrangePseudoFermion);
////////////////////////////////////
// up down action
////////////////////////////////////
std::vector<Real> light_den;
std::vector<Real> light_num;
int n_hasenbusch = hasenbusch.size();
light_den.push_back(light_mass);
for(int h=0;h<n_hasenbusch;h++){
light_den.push_back(hasenbusch[h]);
light_num.push_back(hasenbusch[h]);
}
light_num.push_back(pv_mass);
std::vector<FermionAction *> Numerators;
std::vector<FermionAction *> Denominators;
std::vector<TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy> *> Quotients;
for(int h=0;h<n_hasenbusch+1;h++){
std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl;
Numerators.push_back (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, Params));
Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, Params));
Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],CG,CG));
}
for(int h=0;h<n_hasenbusch+1;h++){
Level1.push_back(Quotients[h]);
}
/////////////////////////////////////////////////////////////
// Gauge action
/////////////////////////////////////////////////////////////
Level2.push_back(&GaugeAction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
std::cout << GridLogMessage << " Action complete "<< std::endl;
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
std::cout << GridLogMessage << " Running the HMC "<< std::endl;
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

73
tests/smearing/Test_smearing.cc Normal file
View File

@ -0,0 +1,73 @@
/*
*
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_cg_prec.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.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>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds);
LatticeFermion src(&Grid); random(pRNG,src);
RealD nrm = norm2(src);
LatticeFermion result(&Grid); result=zero;
LatticeGaugeField Umu(&Grid);
// SU3::HotConfiguration(pRNG,Umu);
SU3::ColdConfiguration(Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
std::vector<int> site({4,4,0,0});
src=zero;
SpinColourVector scv;
scv=zero;
scv()(0)(0) = 1.0;
pokeSite(scv,src,site);
CovariantSmearing<PeriodicGimplR>::GaussianSmear(U, src, 2.0, 50, Tdir);
std::cout << src <<std::endl;
}

194
tests/solver/Test_wilsonclover_mg_lime.cc Normal file
View File

@ -0,0 +1,194 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/solver/Test_wilsonclover_mg_mp.cc
Copyright (C) 2015-2018
Author: Daniel Richtmann <daniel.richtmann@ur.de>
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
*************************************************************************************/
/* */
#include <Grid/Grid.h>
#include <Test_multigrid_common.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main(int argc, char **argv) {
Grid_init(&argc, &argv);
// clang-format off
GridCartesian *FGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
GridCartesian *FGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d);
GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f);
// clang-format on
std::vector<int> fSeeds({1, 2, 3, 4});
GridParallelRNG fPRNG(FGrid_d);
fPRNG.SeedFixedIntegers(fSeeds);
// clang-format off
LatticeFermionD src_d(FGrid_d); gaussian(fPRNG, src_d);
LatticeFermionD resultMGD_d(FGrid_d); resultMGD_d = zero;
LatticeFermionD resultMGF_d(FGrid_d); resultMGF_d = zero;
LatticeGaugeFieldD Umu_d(FGrid_d);
#if 0
{
FieldMetaData header;
std::string file("./qcdsf.769.00399.lime");
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
std::cout <<GridLogMessage<<"** Reading back ILDG conf *********"<<std::endl;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
IldgReader _IldgReader;
_IldgReader.open(file);
_IldgReader.readConfiguration(Umu_d,header);
_IldgReader.close();
}
#else
{
FieldMetaData header;
std::string file("./ckpoint_lat.IEEE64BIG.1100");
NerscIO::readConfiguration(Umu_d,header,file);
}
#endif
// SU3::HotConfiguration(fPRNG, Umu_d);
LatticeGaugeFieldF Umu_f(FGrid_f); precisionChange(Umu_f, Umu_d);
// clang-format on
RealD mass = -0.25;
RealD csw_r = 1.0;
RealD csw_t = 1.0;
MultiGridParams mgParams;
std::string inputXml{"./mg_params.xml"};
if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {
inputXml = GridCmdOptionPayload(argv, argv + argc, "--inputxml");
assert(inputXml.length() != 0);
}
{
XmlWriter writer("mg_params_template.xml");
write(writer, "Params", mgParams);
std::cout << GridLogMessage << "Written mg_params_template.xml" << std::endl;
XmlReader reader(inputXml);
read(reader, "Params", mgParams);
std::cout << GridLogMessage << "Read in " << inputXml << std::endl;
}
checkParameterValidity(mgParams);
std::cout << mgParams << std::endl;
LevelInfo levelInfo_d(FGrid_d, mgParams);
LevelInfo levelInfo_f(FGrid_f, mgParams);
// Note: We do chiral doubling, so actually only nbasis/2 full basis vectors are used
const int nbasis = 40;
WilsonCloverFermionD Dwc_d(Umu_d, *FGrid_d, *FrbGrid_d, mass, csw_r, csw_t);
WilsonCloverFermionF Dwc_f(Umu_f, *FGrid_f, *FrbGrid_f, mass, csw_r, csw_t);
MdagMLinearOperator<WilsonCloverFermionD, LatticeFermionD> MdagMOpDwc_d(Dwc_d);
MdagMLinearOperator<WilsonCloverFermionF, LatticeFermionF> MdagMOpDwc_f(Dwc_f);
std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Testing single-precision Multigrid for Wilson Clover" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl;
auto MGPreconDwc_f = createMGInstance<vSpinColourVectorF, vTComplexF, nbasis, WilsonCloverFermionF>(mgParams, levelInfo_f, Dwc_f, Dwc_f);
MGPreconDwc_f->setup();
if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
MGPreconDwc_f->runChecks(1e-6);
}
MixedPrecisionFlexibleGeneralisedMinimalResidual<LatticeFermionD, LatticeFermionF> MPFGMRESPREC(
1.0e-12, 50000, FGrid_f, *MGPreconDwc_f, 100, false);
std::cout << std::endl << "Starting with a new solver" << std::endl;
MPFGMRESPREC(MdagMOpDwc_d, src_d, resultMGF_d);
MGPreconDwc_f->reportTimings();
if(GridCmdOptionExists(argv, argv + argc, "--docomparison")) {
std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Testing double-precision Multigrid for Wilson Clover" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl;
auto MGPreconDwc_d = createMGInstance<vSpinColourVectorD, vTComplexD, nbasis, WilsonCloverFermionD>(mgParams, levelInfo_d, Dwc_d, Dwc_d);
MGPreconDwc_d->setup();
if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
MGPreconDwc_d->runChecks(1e-13);
}
FlexibleGeneralisedMinimalResidual<LatticeFermionD> FGMRESPREC(1.0e-12, 50000, *MGPreconDwc_d, 100, false);
std::cout << std::endl << "Starting with a new solver" << std::endl;
FGMRESPREC(MdagMOpDwc_d, src_d, resultMGD_d);
MGPreconDwc_d->reportTimings();
std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Comparing single-precision Multigrid with double-precision one for Wilson Clover" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl;
LatticeFermionD diffFullSolver(FGrid_d);
RealD deviationFullSolver = axpy_norm(diffFullSolver, -1.0, resultMGF_d, resultMGD_d);
// clang-format off
LatticeFermionF src_f(FGrid_f); precisionChange(src_f, src_d);
LatticeFermionF resMGF_f(FGrid_f); resMGF_f = zero;
LatticeFermionD resMGD_d(FGrid_d); resMGD_d = zero;
// clang-format on
(*MGPreconDwc_f)(src_f, resMGF_f);
(*MGPreconDwc_d)(src_d, resMGD_d);
LatticeFermionD diffOnlyMG(FGrid_d);
LatticeFermionD resMGF_d(FGrid_d);
precisionChange(resMGF_d, resMGF_f);
RealD deviationOnlyPrec = axpy_norm(diffOnlyMG, -1.0, resMGF_d, resMGD_d);
// clang-format off
std::cout << GridLogMessage << "Absolute difference between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver << std::endl;
std::cout << GridLogMessage << "Relative deviation between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver / norm2(resultMGD_d) << std::endl;
std::cout << GridLogMessage << "Absolute difference between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec << std::endl;
std::cout << GridLogMessage << "Relative deviation between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec / norm2(resMGD_d) << std::endl;
// clang-format on
}
Grid_finalize();
}