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

updated my fork and conflicts fixed

Browse Source
This commit is contained in:
Vera Guelpers 2018-02-14 13:57:56 +00:00
commit de42456171
63 changed files with 4527 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitignore vendored
View File

@ -123,3 +123,9 @@ make-bin-BUCK.sh
#####################
lib/qcd/spin/gamma-gen/*.h
lib/qcd/spin/gamma-gen/*.cc
# vs code editor files #
########################
.vscode/
.vscode/settings.json
settings.json

1
.travis.yml
View File

@ -44,3 +44,4 @@ script:
- make -j4
- ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
- make check

61
extras/Hadrons/Application.cc
View File

@ -43,12 +43,6 @@ using namespace Hadrons;
Application::Application(void)
{
initLogger();
LOG(Message) << "Modules available:" << std::endl;
auto list = ModuleFactory::getInstance().getBuilderList();
for (auto &m: list)
{
LOG(Message) << " " << m << std::endl;
}
auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim);
locVol_ = 1;
for (unsigned int d = 0; d < dim.size(); ++d)
@ -138,24 +132,27 @@ void Application::parseParameterFile(const std::string parameterFileName)
void Application::saveParameterFile(const std::string parameterFileName)
{
XmlWriter writer(parameterFileName);
ObjectId id;
const unsigned int nMod = vm().getNModule();
LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
write(writer, "parameters", getPar());
push(writer, "modules");
for (unsigned int i = 0; i < nMod; ++i)
if (env().getGrid()->IsBoss())
{
push(writer, "module");
id.name = vm().getModuleName(i);
id.type = vm().getModule(i)->getRegisteredName();
write(writer, "id", id);
vm().getModule(i)->saveParameters(writer, "options");
XmlWriter writer(parameterFileName);
ObjectId id;
const unsigned int nMod = vm().getNModule();
write(writer, "parameters", getPar());
push(writer, "modules");
for (unsigned int i = 0; i < nMod; ++i)
{
push(writer, "module");
id.name = vm().getModuleName(i);
id.type = vm().getModule(i)->getRegisteredName();
write(writer, "id", id);
vm().getModule(i)->saveParameters(writer, "options");
pop(writer);
}
pop(writer);
pop(writer);
}
pop(writer);
pop(writer);
}
// schedule computation ////////////////////////////////////////////////////////
@ -170,20 +167,24 @@ void Application::schedule(void)
void Application::saveSchedule(const std::string filename)
{
TextWriter writer(filename);
std::vector<std::string> program;
if (!scheduled_)
{
HADRON_ERROR(Definition, "Computation not scheduled");
}
LOG(Message) << "Saving current schedule to '" << filename << "'..."
<< std::endl;
for (auto address: program_)
if (env().getGrid()->IsBoss())
{
program.push_back(vm().getModuleName(address));
TextWriter writer(filename);
std::vector<std::string> program;
if (!scheduled_)
{
HADRON_ERROR(Definition, "Computation not scheduled");
}
for (auto address: program_)
{
program.push_back(vm().getModuleName(address));
}
write(writer, "schedule", program);
}
write(writer, "schedule", program);
}
void Application::loadSchedule(const std::string filename)

7
extras/Hadrons/Module.hpp
View File

@ -137,6 +137,13 @@ envTmp(type, name, Ls, env().getGrid(Ls))
#define envTmpLat(...)\
MACRO_REDIRECT(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
#define saveResult(ioStem, name, result)\
if (env().getGrid()->IsBoss())\
{\
ResultWriter _writer(RESULT_FILE_NAME(ioStem));\
write(_writer, name, result);\
}
/******************************************************************************
* Module class *
******************************************************************************/

3
extras/Hadrons/Modules.hpp
View File

@ -39,6 +39,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
@ -47,6 +48,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>

153
extras/Hadrons/Modules/MAction/WilsonClover.hpp Normal file
View File

@ -0,0 +1,153 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MAction/Wilson.hpp
Copyright (C) 2015
Copyright (C) 2016
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 */
#ifndef Hadrons_MAction_WilsonClover_hpp_
#define Hadrons_MAction_WilsonClover_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* TWilson quark action *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MAction)
class WilsonCloverPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverPar,
std::string, gauge,
double , mass,
double , csw_r,
double , csw_t,
WilsonAnisotropyCoefficients ,clover_anisotropy,
std::string, boundary
);
};
template <typename FImpl>
class TWilsonClover: public Module<WilsonCloverPar>
{
public:
FGS_TYPE_ALIASES(FImpl,);
public:
// constructor
TWilsonClover(const std::string name);
// destructor
virtual ~TWilsonClover(void) = default;
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(WilsonClover, TWilsonClover<FIMPL>, MAction);
/******************************************************************************
* TWilsonClover template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TWilsonClover<FImpl>::TWilsonClover(const std::string name)
: Module<WilsonCloverPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TWilsonClover<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().gauge};
return in;
}
template <typename FImpl>
std::vector<std::string> TWilsonClover<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWilsonClover<FImpl>::setup(void)
{
//unsigned int size;
// size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
// env().registerObject(getName(), size);
LOG(Message) << "Setting up TWilsonClover fermion matrix with m= " << par().mass
<< " using gauge field '" << par().gauge << "'" << std::endl;
LOG(Message) << "Fermion boundary conditions: " << par().boundary
<< std::endl;
LOG(Message) << "Clover term csw_r: " << par().csw_r
<< " csw_t: " << par().csw_t
<< std::endl;
auto &U = envGet(LatticeGaugeField, par().gauge);
auto &grid = *env().getGrid();
auto &gridRb = *env().getRbGrid();
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid, gridRb, par().mass,
par().csw_r,
par().csw_t,
par().clover_anisotropy,
implParams);
//FMat *fMatPt = new WilsonCloverFermion<FImpl>(U, grid, gridRb, par().mass,
// par().csw_r,
// par().csw_t,
// par().clover_anisotropy,
// implParams);
//env().setObject(getName(), fMatPt);
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWilsonClover<FImpl>::execute()
{
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_WilsonClover_hpp_

3
extras/Hadrons/Modules/MContraction/Baryon.hpp
View File

@ -122,7 +122,6 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
<< " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
<< par().q3 << "'" << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField3, par().q2);
@ -131,7 +130,7 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
// FIXME: do contractions
// write(writer, "meson", result);
// saveResult(par().output, "meson", result);
}
END_MODULE_NAMESPACE

5
extras/Hadrons/Modules/MContraction/DiscLoop.hpp
View File

@ -119,7 +119,6 @@ void TDiscLoop<FImpl>::execute(void)
<< "' using '" << par().q_loop << "' with " << par().gamma
<< " insertion." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q_loop = envGet(PropagatorField, par().q_loop);
Gamma gamma(par().gamma);
std::vector<TComplex> buf;
@ -128,15 +127,13 @@ void TDiscLoop<FImpl>::execute(void)
envGetTmp(LatticeComplex, c);
c = trace(gamma*q_loop);
sliceSum(c, buf, Tp);
result.gamma = par().gamma;
result.corr.resize(buf.size());
for (unsigned int t = 0; t < buf.size(); ++t)
{
result.corr[t] = TensorRemove(buf[t]);
}
write(writer, "disc", result);
saveResult(par().output, "disc", result);
}
END_MODULE_NAMESPACE

4
extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
View File

@ -153,7 +153,6 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
// Initialise variables. q2 and q3 are normal propagators, q1 may be
// sink smeared.
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(SlicedPropagator1, par().q1);
auto &q2 = envGet(PropagatorField2, par().q2);
auto &q3 = envGet(PropagatorField2, par().q3);
@ -175,8 +174,7 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
{
result.corr[t] = TensorRemove(buf[t]);
}
write(writer, "gamma3pt", result);
saveResult(par().output, "gamma3pt", result);
}
END_MODULE_NAMESPACE

3
extras/Hadrons/Modules/MContraction/Meson.hpp
View File

@ -172,7 +172,6 @@ void TMeson<FImpl1, FImpl2>::execute(void)
<< " quarks '" << par().q1 << "' and '" << par().q2 << "'"
<< std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
std::vector<TComplex> buf;
std::vector<Result> result;
Gamma g5(Gamma::Algebra::Gamma5);
@ -239,7 +238,7 @@ void TMeson<FImpl1, FImpl2>::execute(void)
}
}
}
write(writer, "meson", result);
saveResult(par().output, "meson", result);
}
END_MODULE_NAMESPACE

4
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
View File

@ -104,7 +104,6 @@ void TWeakHamiltonianEye::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(SlicedPropagator, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -147,5 +146,6 @@ void TWeakHamiltonianEye::execute(void)
SUM_MU(expbuf, E_body[mu]*E_loop[mu])
MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E")
write(writer, "HW_Eye", result);
// IO
saveResult(par().output, "HW_Eye", result);
}

4
extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
View File

@ -104,7 +104,6 @@ void TWeakHamiltonianNonEye::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -144,5 +143,6 @@ void TWeakHamiltonianNonEye::execute(void)
SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu])
MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W")
write(writer, "HW_NonEye", result);
// IO
saveResult(par().output, "HW_NonEye", result);
}

4
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
View File

@ -104,7 +104,6 @@ void TWeakNeutral4ptDisc::execute(void)
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
<< "'." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &q3 = envGet(PropagatorField, par().q3);
@ -138,5 +137,6 @@ void TWeakNeutral4ptDisc::execute(void)
expbuf *= curr;
MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")
write(writer, "HW_disc0", result);
// IO
saveResult(par().output, "HW_disc0", result);
}

13
extras/Hadrons/Modules/MFermion/GaugeProp.hpp
View File

@ -94,7 +94,6 @@ private:
};
MODULE_REGISTER_NS(GaugeProp, TGaugeProp<FIMPL>, MFermion);
/******************************************************************************
* TGaugeProp implementation *
******************************************************************************/
@ -154,7 +153,7 @@ void TGaugeProp<FImpl>::execute(void)
LOG(Message) << "Inverting using solver '" << par().solver
<< "' on source '" << par().source << "'" << std::endl;
for (unsigned int s = 0; s < Ns; ++s)
for (unsigned int c = 0; c < Nc; ++c)
for (unsigned int c = 0; c < FImpl::Dimension; ++c)
{
LOG(Message) << "Inversion for spin= " << s << ", color= " << c
<< std::endl;
@ -163,11 +162,11 @@ void TGaugeProp<FImpl>::execute(void)
{
if (Ls_ == 1)
{
PropToFerm(source, fullSrc, s, c);
PropToFerm<FImpl>(source, fullSrc, s, c);
}
else
{
PropToFerm(tmp, fullSrc, s, c);
PropToFerm<FImpl>(tmp, fullSrc, s, c);
make_5D(tmp, source, Ls_);
}
}
@ -180,18 +179,18 @@ void TGaugeProp<FImpl>::execute(void)
}
else
{
PropToFerm(source, fullSrc, s, c);
PropToFerm<FImpl>(source, fullSrc, s, c);
}
}
sol = zero;
solver(sol, source);
FermToProp(prop, sol, s, c);
FermToProp<FImpl>(prop, sol, s, c);
// create 4D propagators from 5D one if necessary
if (Ls_ > 1)
{
PropagatorField &p4d = envGet(PropagatorField, getName());
make_4D(sol, tmp, Ls_);
FermToProp(p4d, tmp, s, c);
FermToProp<FImpl>(p4d, tmp, s, c);
}
}
}

75
extras/Hadrons/Modules/MGauge/FundtoHirep.cc Normal file
View File

@ -0,0 +1,75 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc
Copyright (C) 2015
Copyright (C) 2016
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/Hadrons/Modules/MGauge/FundtoHirep.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MGauge;
// constructor /////////////////////////////////////////////////////////////////
template <class Rep>
TFundtoHirep<Rep>::TFundtoHirep(const std::string name)
: Module<FundtoHirepPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <class Rep>
std::vector<std::string> TFundtoHirep<Rep>::getInput(void)
{
std::vector<std::string> in;
return in;
}
template <class Rep>
std::vector<std::string> TFundtoHirep<Rep>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename Rep>
void TFundtoHirep<Rep>::setup(void)
{
env().template registerLattice<typename Rep::LatticeField>(getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <class Rep>
void TFundtoHirep<Rep>::execute(void)
{
auto &U = *env().template getObject<LatticeGaugeField>(par().gaugeconf);
LOG(Message) << "Transforming Representation" << std::endl;
Rep TargetRepresentation(U._grid);
TargetRepresentation.update_representation(U);
typename Rep::LatticeField &URep = *env().template createLattice<typename Rep::LatticeField>(getName());
URep = TargetRepresentation.U;
}

77
extras/Hadrons/Modules/MGauge/FundtoHirep.hpp Normal file
View File

@ -0,0 +1,77 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.hpp
Copyright (C) 2015
Copyright (C) 2016
Author: David Preti <david.preti@to.infn.it>
Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MGauge_FundtoHirep_hpp_
#define Hadrons_MGauge_FundtoHirep_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Load a NERSC configuration *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MGauge)
class FundtoHirepPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(FundtoHirepPar,
std::string, gaugeconf);
};
template <class Rep>
class TFundtoHirep: public Module<FundtoHirepPar>
{
public:
// constructor
TFundtoHirep(const std::string name);
// destructor
virtual ~TFundtoHirep(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
void setup(void);
// execution
void execute(void);
};
//MODULE_REGISTER_NS(FundtoAdjoint, TFundtoHirep<AdjointRepresentation>, MGauge);
//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MGauge_FundtoHirep_hpp_

5
extras/Hadrons/Modules/MScalar/ChargedProp.cc
View File

@ -133,7 +133,6 @@ void TChargedProp::execute(void)
LOG(Message) << "Saving zero-momentum projection to '"
<< filename << "'..." << std::endl;
ResultWriter writer(RESULT_FILE_NAME(par().output));
std::vector<TComplex> vecBuf;
std::vector<Complex> result;
@ -143,8 +142,8 @@ void TChargedProp::execute(void)
{
result[t] = TensorRemove(vecBuf[t]);
}
write(writer, "charge", q);
write(writer, "prop", result);
saveResult(par().output, "charge", q);
saveResult(par().output, "prop", result);
}
}

4
extras/Hadrons/Modules/MScalar/FreeProp.cc
View File

@ -83,8 +83,6 @@ void TFreeProp::execute(void)
if (!par().output.empty())
{
TextWriter writer(par().output + "." +
std::to_string(vm().getTrajectory()));
std::vector<TComplex> buf;
std::vector<Complex> result;
@ -94,6 +92,6 @@ void TFreeProp::execute(void)
{
result[t] = TensorRemove(buf[t]);
}
write(writer, "prop", result);
saveResult(par().output, "freeprop", result);
}
}

165
extras/Hadrons/Modules/MScalarSUN/Div.hpp Normal file
View File

@ -0,0 +1,165 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MScalarSUN/Div.hpp
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 */
#ifndef Hadrons_MScalarSUN_Div_hpp_
#define Hadrons_MScalarSUN_Div_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Div *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalarSUN)
class DivPar: Serializable
{
public:
GRID_SERIALIZABLE_ENUM(DiffType, undef, forward, 1, backward, 2, central, 3);
GRID_SERIALIZABLE_CLASS_MEMBERS(DivPar,
std::vector<std::string>, op,
DiffType, type,
std::string, output);
};
template <typename SImpl>
class TDiv: public Module<DivPar>
{
public:
typedef typename SImpl::Field Field;
typedef typename SImpl::ComplexField ComplexField;
class Result: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
DivPar::DiffType, type,
Complex, value);
};
public:
// constructor
TDiv(const std::string name);
// destructor
virtual ~TDiv(void) = default;
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_NS(DivSU2, TDiv<ScalarNxNAdjImplR<2>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU3, TDiv<ScalarNxNAdjImplR<3>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU4, TDiv<ScalarNxNAdjImplR<4>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
MODULE_REGISTER_NS(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
/******************************************************************************
* TDiv implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename SImpl>
TDiv<SImpl>::TDiv(const std::string name)
: Module<DivPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename SImpl>
std::vector<std::string> TDiv<SImpl>::getInput(void)
{
return par().op;
}
template <typename SImpl>
std::vector<std::string> TDiv<SImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename SImpl>
void TDiv<SImpl>::setup(void)
{
if (par().op.size() != env().getNd())
{
HADRON_ERROR(Size, "the number of components differs from number of dimensions");
}
envCreateLat(ComplexField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
template <typename SImpl>
void TDiv<SImpl>::execute(void)
{
const auto nd = env().getNd();
LOG(Message) << "Computing the " << par().type << " divergence of [";
for (unsigned int mu = 0; mu < nd; ++mu)
{
std::cout << par().op[mu] << ((mu == nd - 1) ? "]" : ", ");
}
std::cout << std::endl;
auto &div = envGet(ComplexField, getName());
div = zero;
for (unsigned int mu = 0; mu < nd; ++mu)
{
auto &op = envGet(ComplexField, par().op[mu]);
switch(par().type)
{
case DivPar::DiffType::backward:
div += op - Cshift(op, mu, -1);
break;
case DivPar::DiffType::forward:
div += Cshift(op, mu, 1) - op;
break;
case DivPar::DiffType::central:
div += 0.5*(Cshift(op, mu, 1) - Cshift(op, mu, -1));
break;
}
}
if (!par().output.empty())
{
Result r;
r.type = par().type;
r.value = TensorRemove(sum(div));
saveResult(par().output, "div", r);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalarSUN_Div_hpp_

3
extras/Hadrons/Modules/MScalarSUN/TrMag.hpp
View File

@ -120,7 +120,6 @@ void TTrMag<SImpl>::execute(void)
<< "..." << std::endl;
std::vector<Result> result;
ResultWriter writer(RESULT_FILE_NAME(par().output));
auto &phi = envGet(Field, par().field);
auto m2 = sum(phi), mn = m2;
@ -136,7 +135,7 @@ void TTrMag<SImpl>::execute(void)
r.value = TensorRemove(trace(mn)).real();
result.push_back(r);
}
write(writer, "trmag", result);
saveResult(par().output, "trmag", result);
}
END_MODULE_NAMESPACE

4
extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp
View File

@ -162,9 +162,7 @@ void TTrPhi<SImpl>::execute(void)
}
if (result.size() > 0)
{
ResultWriter writer(RESULT_FILE_NAME(par().output));
write(writer, "trphi", result);
saveResult(par().output, "trphi", result);
}
}

3
extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp
View File

@ -129,7 +129,6 @@ void TTwoPoint<SImpl>::execute(void)
LOG(Message) << " '" << o << "'" << std::endl;
}
ResultWriter writer(RESULT_FILE_NAME(par().output));
const unsigned int nd = env().getDim().size();
std::vector<Result> result;
@ -150,7 +149,7 @@ void TTwoPoint<SImpl>::execute(void)
r.data = makeTwoPoint(slicedOp[i], slicedOp[j]);
result.push_back(r);
}
write(writer, "twopt", result);
saveResult(par().output, "twopt", result);
}
// make 2-pt function //////////////////////////////////////////////////////////

6
extras/Hadrons/modules.inc
View File

@ -1,10 +1,13 @@
modules_cc =\
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/FreeProp.cc \
Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MContraction/WeakHamiltonianNonEye.cc \
Modules/MGauge/Unit.cc \
Modules/MGauge/StochEm.cc \
Modules/MGauge/Random.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MScalar/FreeProp.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MIO/LoadNersc.cc
@ -31,6 +34,7 @@ modules_hpp =\
Modules/MGauge/Unit.hpp \
Modules/MGauge/Random.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/FundtoHirep.hpp \
Modules/MUtilities/TestSeqGamma.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MLoop/NoiseLoop.hpp \
@ -39,6 +43,8 @@ modules_hpp =\
Modules/MScalar/ChargedProp.hpp \
Modules/MAction/DWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TrPhi.hpp \

10
lib/algorithms/LinearOperator.h
View File

@ -183,11 +183,13 @@ namespace Grid {
virtual RealD Mpc (const Field &in, Field &out) =0;
virtual RealD MpcDag (const Field &in, Field &out) =0;
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
Field tmp(in._grid);
Field tmp(in._grid);
tmp.checkerboard = in.checkerboard;
ni=Mpc(in,tmp);
no=MpcDag(tmp,out);
}
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
out.checkerboard = in.checkerboard;
MpcDagMpc(in,out,n1,n2);
}
virtual void HermOp(const Field &in, Field &out){
@ -215,13 +217,15 @@ namespace Grid {
public:
SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in._grid);
// std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl;
Field tmp(in._grid);
tmp.checkerboard = !in.checkerboard;
//std::cout <<"grid pointers: in._grid="<< in._grid << " out._grid=" << out._grid << " _Mat.Grid=" << _Mat.Grid() << " _Mat.RedBlackGrid=" << _Mat.RedBlackGrid() << std::endl;
_Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
//std::cout << "cb in " << in.checkerboard << " cb out " << out.checkerboard << std::endl;
_Mat.Mooee(in,out);
return axpy_norm(out,-1.0,tmp,out);
}

2
lib/cartesian/Cartesian_base.h
View File

@ -79,6 +79,8 @@ public:
std::vector<int> _lstart; // local start of array in gcoors _processor_coor[d]*_ldimensions[d]
std::vector<int> _lend ; // local end of array in gcoors _processor_coor[d]*_ldimensions[d]+_ldimensions_[d]-1
bool _isCheckerBoarded;
public:
////////////////////////////////////////////////////////////////

1
lib/cartesian/Cartesian_full.h
View File

@ -97,6 +97,7 @@ public:
///////////////////////
// Grid information
///////////////////////
_isCheckerBoarded = false;
_ndimension = dimensions.size();
_fdimensions.resize(_ndimension);

5
lib/cartesian/Cartesian_red_black.h
View File

@ -171,9 +171,8 @@ public:
const std::vector<int> &checker_dim_mask,
int checker_dim)
{
///////////////////////
// Grid information
///////////////////////
_isCheckerBoarded = true;
_checker_dim = checker_dim;
assert(checker_dim_mask[checker_dim] == 1);
_ndimension = dimensions.size();

49
lib/qcd/QCD.h
View File

@ -39,6 +39,7 @@ namespace QCD {
static const int Zdir = 2;
static const int Tdir = 3;
static const int Xp = 0;
static const int Yp = 1;
static const int Zp = 2;
@ -420,15 +421,16 @@ namespace QCD {
//////////////////////////////////////////////
// Fermion <-> propagator assignements
//////////////////////////////////////////////
template <class Prop, class Ferm>
void FermToProp(Prop &p, const Ferm &f, const int s, const int c)
//template <class Prop, class Ferm>
template <class Fimpl>
void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c)
{
for(int j = 0; j < Ns; ++j)
for(int j = 0; j < Ns; ++j)
{
auto pjs = peekSpin(p, j, s);
auto fj = peekSpin(f, j);
for(int i = 0; i < Nc; ++i)
for(int i = 0; i < Fimpl::Dimension; ++i)
{
pokeColour(pjs, peekColour(fj, i), i, c);
}
@ -436,15 +438,16 @@ namespace QCD {
}
}
template <class Prop, class Ferm>
void PropToFerm(Ferm &f, const Prop &p, const int s, const int c)
//template <class Prop, class Ferm>
template <class Fimpl>
void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c)
{
for(int j = 0; j < Ns; ++j)
{
auto pjs = peekSpin(p, j, s);
auto fj = peekSpin(f, j);
for(int i = 0; i < Nc; ++i)
for(int i = 0; i < Fimpl::Dimension; ++i)
{
pokeColour(fj, peekColour(pjs, i, c), i);
}
@ -503,38 +506,6 @@ namespace QCD {
} //namespace QCD
} // Grid
/*
<<<<<<< HEAD
#include <Grid/qcd/utils/SpaceTimeGrid.h>
#include <Grid/qcd/spin/Dirac.h>
#include <Grid/qcd/spin/TwoSpinor.h>
#include <Grid/qcd/utils/LinalgUtils.h>
#include <Grid/qcd/utils/CovariantCshift.h>
// Include representations
#include <Grid/qcd/utils/SUn.h>
#include <Grid/qcd/utils/SUnAdjoint.h>
#include <Grid/qcd/utils/SUnTwoIndex.h>
#include <Grid/qcd/representations/hmc_types.h>
// Scalar field
#include <Grid/qcd/utils/ScalarObjs.h>
#include <Grid/qcd/action/Actions.h>
#include <Grid/qcd/smearing/Smearing.h>
#include <Grid/qcd/hmc/integrators/Integrator.h>
#include <Grid/qcd/hmc/integrators/Integrator_algorithm.h>
#include <Grid/qcd/observables/hmc_observable.h>
#include <Grid/qcd/hmc/HMC.h>
//#include <Grid/qcd/modules/mods.h>
=======
>>>>>>> develop
*/
#endif

29
lib/qcd/action/fermion/Fermion.h
View File

@ -50,11 +50,13 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
////////////////////////////////////////////
#include <Grid/qcd/action/fermion/WilsonFermion.h> // 4d wilson like
#include <Grid/qcd/action/fermion/WilsonTMFermion.h> // 4d wilson like
#include <Grid/qcd/action/fermion/WilsonTMFermion.h> // 4d wilson like
#include <Grid/qcd/action/fermion/WilsonCloverFermion.h> // 4d wilson clover fermions
#include <Grid/qcd/action/fermion/WilsonFermion5D.h> // 5d base used by all 5d overlap types
//#include <Grid/qcd/action/fermion/CloverFermion.h>
#include <Grid/qcd/action/fermion/ImprovedStaggeredFermion.h>
#include <Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h>
#include <Grid/qcd/action/fermion/CayleyFermion5D.h> // Cayley types
#include <Grid/qcd/action/fermion/DomainWallFermion.h>
#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
@ -104,10 +106,33 @@ typedef WilsonFermion<WilsonTwoIndexSymmetricImplR> WilsonTwoIndexSymmetricFermi
typedef WilsonFermion<WilsonTwoIndexSymmetricImplF> WilsonTwoIndexSymmetricFermionF;
typedef WilsonFermion<WilsonTwoIndexSymmetricImplD> WilsonTwoIndexSymmetricFermionD;
typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonTwoIndexAntiSymmetricFermionR;
typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonTwoIndexAntiSymmetricFermionF;
typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonTwoIndexAntiSymmetricFermionD;
// Twisted mass fermion
typedef WilsonTMFermion<WilsonImplR> WilsonTMFermionR;
typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
// Clover fermions
typedef WilsonCloverFermion<WilsonImplR> WilsonCloverFermionR;
typedef WilsonCloverFermion<WilsonImplF> WilsonCloverFermionF;
typedef WilsonCloverFermion<WilsonImplD> WilsonCloverFermionD;
typedef WilsonCloverFermion<WilsonAdjImplR> WilsonCloverAdjFermionR;
typedef WilsonCloverFermion<WilsonAdjImplF> WilsonCloverAdjFermionF;
typedef WilsonCloverFermion<WilsonAdjImplD> WilsonCloverAdjFermionD;
typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
// Domain Wall fermions
typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;
typedef DomainWallFermion<WilsonImplF> DomainWallFermionF;
typedef DomainWallFermion<WilsonImplD> DomainWallFermionD;

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

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

145
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -164,6 +164,7 @@ namespace QCD {
public:
static const int Dimension = Representation::Dimension;
static const bool isFundamental = Representation::isFundamental;
static const bool LsVectorised=false;
static const int Nhcs = Options::Nhcs;
@ -261,8 +262,22 @@ namespace QCD {
GaugeLinkField link(mat._grid);
link = TraceIndex<SpinIndex>(outerProduct(Btilde,A));
PokeIndex<LorentzIndex>(mat,link,mu);
}
}
inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){
mat = outerProduct(B,A);
}
inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
mat = TraceIndex<SpinIndex>(P);
}
inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
for (int mu = 0; mu < Nd; mu++)
mat[mu] = PeekIndex<LorentzIndex>(Uds, mu);
}
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
int Ls=Btilde._grid->_fdimensions[0];
@ -284,27 +299,28 @@ namespace QCD {
////////////////////////////////////////////////////////////////////////////////////
// Single flavour four spinors with colour index, 5d redblack
////////////////////////////////////////////////////////////////////////////////////
template<class S,int Nrepresentation=Nc, class Options=CoeffReal>
class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > {
template<class S,class Representation = FundamentalRepresentation, class Options=CoeffReal>
class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Representation::Dimension> > {
public:
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > Gimpl;
INHERIT_GIMPL_TYPES(Gimpl);
static const int Dimension = Nrepresentation;
static const int Dimension = Representation::Dimension;
static const bool isFundamental = Representation::isFundamental;
static const bool LsVectorised=true;
static const int Nhcs = Options::Nhcs;
typedef typename Options::_Coeff_t Coeff_t;
typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
template <typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
template <typename vtype> using iImplPropagator = iScalar<iMatrix<iMatrix<vtype, Nrepresentation>, Ns> >;
template <typename vtype> using iImplHalfSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
template <typename vtype> using iImplHalfCommSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhcs> >;
template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
template <typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
template <typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
template <typename vtype> using iImplPropagator = iScalar<iMatrix<iMatrix<vtype, Dimension>, Ns> >;
template <typename vtype> using iImplHalfSpinor = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
template <typename vtype> using iImplHalfCommSpinor = iScalar<iVector<iVector<vtype, Dimension>, Nhcs> >;
template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nd>;
template <typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
typedef iImplSpinor<Simd> SiteSpinor;
typedef iImplPropagator<Simd> SitePropagator;
@ -340,8 +356,8 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
StencilImpl &St) {
SiteGaugeLink UU;
for (int i = 0; i < Nrepresentation; i++) {
for (int j = 0; j < Nrepresentation; j++) {
for (int i = 0; i < Dimension; i++) {
for (int j = 0; j < Dimension; j++) {
vsplat(UU()()(i, j), U(mu)()(i, j));
}
}
@ -353,8 +369,8 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
const SitePropagator &chi,
int mu) {
SiteGaugeLink UU;
for (int i = 0; i < Nrepresentation; i++) {
for (int j = 0; j < Nrepresentation; j++) {
for (int i = 0; i < Dimension; i++) {
for (int j = 0; j < Dimension; j++) {
vsplat(UU()()(i, j), U(mu)()(i, j));
}
}
@ -393,6 +409,19 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
assert(0);
}
inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
assert(0);
}
inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
assert(0);
}
inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
assert(0);
}
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
assert(0);
@ -445,25 +474,26 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
////////////////////////////////////////////////////////////////////////////////////////
// Flavour doubled spinors; is Gparity the only? what about C*?
////////////////////////////////////////////////////////////////////////////////////////
template <class S, int Nrepresentation, class Options=CoeffReal>
class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
public:
static const int Dimension = Nrepresentation;
static const int Dimension = Representation::Dimension;
static const bool isFundamental = Representation::isFundamental;
static const int Nhcs = Options::Nhcs;
static const bool LsVectorised=false;
typedef ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> > Gimpl;
typedef ConjugateGaugeImpl< GaugeImplTypes<S,Dimension> > Gimpl;
INHERIT_GIMPL_TYPES(Gimpl);
typedef typename Options::_Coeff_t Coeff_t;
typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
template <typename vtype> using iImplSpinor = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp>;
template <typename vtype> using iImplPropagator = iVector<iMatrix<iMatrix<vtype, Nrepresentation>, Ns>, Ngp>;
template <typename vtype> using iImplHalfSpinor = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp>;
template <typename vtype> using iImplHalfCommSpinor = iVector<iVector<iVector<vtype, Nrepresentation>, Nhcs>, Ngp>;
template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
template <typename vtype> using iImplSpinor = iVector<iVector<iVector<vtype, Dimension>, Ns>, Ngp>;
template <typename vtype> using iImplPropagator = iVector<iMatrix<iMatrix<vtype, Dimension>, Ns>, Ngp>;
template <typename vtype> using iImplHalfSpinor = iVector<iVector<iVector<vtype, Dimension>, Nhs>, Ngp>;
template <typename vtype> using iImplHalfCommSpinor = iVector<iVector<iVector<vtype, Dimension>, Nhcs>, Ngp>;
template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>, Ngp>;
typedef iImplSpinor<Simd> SiteSpinor;
typedef iImplPropagator<Simd> SitePropagator;
@ -636,6 +666,25 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
return;
}
inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
//mat = outerProduct(Btilde, A);
assert(0);
}
inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
assert(0);
/*
auto tmp = TraceIndex<SpinIndex>(P);
parallel_for(auto ss = tmp.begin(); ss < tmp.end(); ss++) {
mat[ss]() = tmp[ss](0, 0) + conjugate(tmp[ss](1, 1));
}
*/
}
inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
assert(0);
}
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
int Ls = Btilde._grid->_fdimensions[0];
@ -665,6 +714,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
typedef RealD _Coeff_t ;
static const int Dimension = Representation::Dimension;
static const bool isFundamental = Representation::isFundamental;
static const bool LsVectorised=false;
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
@ -776,8 +826,8 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
GaugeLinkField link(mat._grid);
link = TraceIndex<SpinIndex>(outerProduct(Btilde,A));
PokeIndex<LorentzIndex>(mat,link,mu);
}
}
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
assert (0);
// Must never hit
@ -793,6 +843,7 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
public:
static const int Dimension = Representation::Dimension;
static const bool isFundamental = Representation::isFundamental;
static const bool LsVectorised=true;
typedef RealD Coeff_t ;
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
@ -983,29 +1034,33 @@ typedef WilsonImpl<vComplex, TwoIndexSymmetricRepresentation, CoeffReal > Wilso
typedef WilsonImpl<vComplexF, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplF; // Float
typedef WilsonImpl<vComplexD, TwoIndexSymmetricRepresentation, CoeffReal > WilsonTwoIndexSymmetricImplD; // Double
typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffReal> DomainWallVec5dImplF; // Float
typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffReal> DomainWallVec5dImplD; // Double
typedef WilsonImpl<vComplex, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplR; // Real.. whichever prec
typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF; // Float
typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD; // Double
typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplR; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplF; // Float
typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffReal> DomainWallVec5dImplD; // Double
typedef DomainWallVec5dImpl<vComplex ,Nc, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,Nc, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
typedef DomainWallVec5dImpl<vComplexD,Nc, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplRL; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplFH; // Float
typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation, CoeffRealHalfComms> DomainWallVec5dImplDF; // Double
typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplex> ZDomainWallVec5dImplF; // Float
typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplex> ZDomainWallVec5dImplD; // Double
typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplR; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplF; // Float
typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplex> ZDomainWallVec5dImplD; // Double
typedef DomainWallVec5dImpl<vComplex ,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
typedef DomainWallVec5dImpl<vComplexD,Nc,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
typedef DomainWallVec5dImpl<vComplex ,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplRL; // Real.. whichever prec
typedef DomainWallVec5dImpl<vComplexF,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplFH; // Float
typedef DomainWallVec5dImpl<vComplexD,FundamentalRepresentation,CoeffComplexHalfComms> ZDomainWallVec5dImplDF; // Double
typedef GparityWilsonImpl<vComplex , Nc,CoeffReal> GparityWilsonImplR; // Real.. whichever prec
typedef GparityWilsonImpl<vComplexF, Nc,CoeffReal> GparityWilsonImplF; // Float
typedef GparityWilsonImpl<vComplexD, Nc,CoeffReal> GparityWilsonImplD; // Double
typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffReal> GparityWilsonImplR; // Real.. whichever prec
typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffReal> GparityWilsonImplF; // Float
typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffReal> GparityWilsonImplD; // Double
typedef GparityWilsonImpl<vComplex , Nc,CoeffRealHalfComms> GparityWilsonImplRL; // Real.. whichever prec
typedef GparityWilsonImpl<vComplexF, Nc,CoeffRealHalfComms> GparityWilsonImplFH; // Float
typedef GparityWilsonImpl<vComplexD, Nc,CoeffRealHalfComms> GparityWilsonImplDF; // Double
typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplRL; // Real.. whichever prec
typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplFH; // Float
typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplDF; // Double
typedef StaggeredImpl<vComplex, FundamentalRepresentation > StaggeredImplR; // Real.. whichever prec
typedef StaggeredImpl<vComplexF, FundamentalRepresentation > StaggeredImplF; // Float

243
lib/qcd/action/fermion/WilsonCloverFermion.cc Normal file
View File

@ -0,0 +1,243 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
Copyright (C) 2017
Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/Eigen/Dense>
#include <Grid/qcd/spin/Dirac.h>
namespace Grid
{
namespace QCD
{
// *NOT* EO
template <class Impl>
RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
{
FermionField temp(out._grid);
// Wilson term
out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerNo);
// Clover term
Mooee(in, temp);
out += temp;
return norm2(out);
}
template <class Impl>
RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
{
FermionField temp(out._grid);
// Wilson term
out.checkerboard = in.checkerboard;
this->Dhop(in, out, DaggerYes);
// Clover term
MooeeDag(in, temp);
out += temp;
return norm2(out);
}
template <class Impl>
void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
{
WilsonFermion<Impl>::ImportGauge(_Umu);
GridBase *grid = _Umu._grid;
typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
// Compute the field strength terms mu>nu
WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
// Compute the Clover Operator acting on Colour and Spin
// multiply here by the clover coefficients for the anisotropy
CloverTerm = fillCloverYZ(Bx) * csw_r;
CloverTerm += fillCloverXZ(By) * csw_r;
CloverTerm += fillCloverXY(Bz) * csw_r;
CloverTerm += fillCloverXT(Ex) * csw_t;
CloverTerm += fillCloverYT(Ey) * csw_t;
CloverTerm += fillCloverZT(Ez) * csw_t;
CloverTerm += diag_mass;
int lvol = _Umu._grid->lSites();
int DimRep = Impl::Dimension;
Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
std::vector<int> lcoor;
typename SiteCloverType::scalar_object Qx = zero, Qxinv = zero;
for (int site = 0; site < lvol; site++)
{
grid->LocalIndexToLocalCoor(site, lcoor);
EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
peekLocalSite(Qx, CloverTerm, lcoor);
Qxinv = zero;
//if (csw!=0){
for (int j = 0; j < Ns; j++)
for (int k = 0; k < Ns; k++)
for (int a = 0; a < DimRep; a++)
for (int b = 0; b < DimRep; b++)
EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b);
// if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
EigenInvCloverOp = EigenCloverOp.inverse();
//std::cout << EigenInvCloverOp << std::endl;
for (int j = 0; j < Ns; j++)
for (int k = 0; k < Ns; k++)
for (int a = 0; a < DimRep; a++)
for (int b = 0; b < DimRep; b++)
Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
// if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
// }
pokeLocalSite(Qxinv, CloverTermInv, lcoor);
}
// Separate the even and odd parts
pickCheckerboard(Even, CloverTermEven, CloverTerm);
pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm));
pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm));
pickCheckerboard(Even, CloverTermInvEven, CloverTermInv);
pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv);
pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
}
template <class Impl>
void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerNo, InverseNo);
}
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerYes, InverseNo);
}
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerNo, InverseYes);
}
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
{
this->MooeeInternal(in, out, DaggerYes, InverseYes);
}
template <class Impl>
void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
{
out.checkerboard = in.checkerboard;
CloverFieldType *Clover;
assert(in.checkerboard == Odd || in.checkerboard == Even);
if (dag)
{
if (in._grid->_isCheckerBoarded)
{
if (in.checkerboard == Odd)
{
Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd;
}
else
{
Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
}
out = *Clover * in;
}
else
{
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = adj(*Clover) * in;
}
}
else
{
if (in._grid->_isCheckerBoarded)
{
if (in.checkerboard == Odd)
{
// std::cout << "Calling clover term Odd" << std::endl;
Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
}
else
{
// std::cout << "Calling clover term Even" << std::endl;
Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
}
out = *Clover * in;
// std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl;
}
else
{
Clover = (inv) ? &CloverTermInv : &CloverTerm;
out = *Clover * in;
}
}
} // MooeeInternal
// Derivative parts
template <class Impl>
void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
{
assert(0);
}
// Derivative parts
template <class Impl>
void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
{
assert(0); // not implemented yet
}
FermOpTemplateInstantiate(WilsonCloverFermion);
AdjointFermOpTemplateInstantiate(WilsonCloverFermion);
TwoIndexFermOpTemplateInstantiate(WilsonCloverFermion);
//GparityFermOpTemplateInstantiate(WilsonCloverFermion);
}
}

366
lib/qcd/action/fermion/WilsonCloverFermion.h Normal file
View File

@ -0,0 +1,366 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
Author: David Preti <>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef GRID_QCD_WILSON_CLOVER_FERMION_H
#define GRID_QCD_WILSON_CLOVER_FERMION_H
#include <Grid/Grid.h>
namespace Grid
{
namespace QCD
{
///////////////////////////////////////////////////////////////////
// Wilson Clover
//
// Operator ( with anisotropy coefficients):
//
// Q = 1 + (Nd-1)/xi_0 + m
// + W_t + (nu/xi_0) * W_s
// - 1/2*[ csw_t * sum_s (sigma_ts F_ts) + (csw_s/xi_0) * sum_ss (sigma_ss F_ss) ]
//
// s spatial, t temporal directions.
// where W_t and W_s are the temporal and spatial components of the
// Wilson Dirac operator
//
// csw_r = csw_t to recover the isotropic version
//////////////////////////////////////////////////////////////////
template <class Impl>
class WilsonCloverFermion : public WilsonFermion<Impl>
{
public:
// Types definitions
INHERIT_IMPL_TYPES(Impl);
template <typename vtype>
using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
typedef iImplClover<Simd> SiteCloverType;
typedef Lattice<SiteCloverType> CloverFieldType;
public:
typedef WilsonFermion<Impl> WilsonBase;
virtual void Instantiatable(void){};
// Constructors
WilsonCloverFermion(GaugeField &_Umu, GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid,
const RealD _mass,
const RealD _csw_r = 0.0,
const RealD _csw_t = 0.0,
const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(),
const ImplParams &impl_p = ImplParams()) : WilsonFermion<Impl>(_Umu,
Fgrid,
Hgrid,
_mass, impl_p, clover_anisotropy),
CloverTerm(&Fgrid),
CloverTermInv(&Fgrid),
CloverTermEven(&Hgrid),
CloverTermOdd(&Hgrid),
CloverTermInvEven(&Hgrid),
CloverTermInvOdd(&Hgrid),
CloverTermDagEven(&Hgrid),
CloverTermDagOdd(&Hgrid),
CloverTermInvDagEven(&Hgrid),
CloverTermInvDagOdd(&Hgrid)
{
assert(Nd == 4); // require 4 dimensions
if (clover_anisotropy.isAnisotropic)
{
csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0;
diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
}
else
{
csw_r = _csw_r * 0.5;
diag_mass = 4.0 + _mass;
}
csw_t = _csw_t * 0.5;
if (csw_r == 0)
std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl;
if (csw_t == 0)
std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl;
ImportGauge(_Umu);
}
virtual RealD M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out);
virtual void Mooee(const FermionField &in, FermionField &out);
virtual void MooeeDag(const FermionField &in, FermionField &out);
virtual void MooeeInv(const FermionField &in, FermionField &out);
virtual void MooeeInvDag(const FermionField &in, FermionField &out);
virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv);
//virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
virtual void MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
virtual void MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
void ImportGauge(const GaugeField &_Umu);
// Derivative parts unpreconditioned pseudofermions
void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
{
conformable(X._grid, Y._grid);
conformable(X._grid, force._grid);
GaugeLinkField force_mu(force._grid), lambda(force._grid);
GaugeField clover_force(force._grid);
PropagatorField Lambda(force._grid);
// Guido: Here we are hitting some performance issues:
// need to extract the components of the DoubledGaugeField
// for each call
// Possible solution
// Create a vector object to store them? (cons: wasting space)
std::vector<GaugeLinkField> U(Nd, this->Umu._grid);
Impl::extractLinkField(U, this->Umu);
force = zero;
// Derivative of the Wilson hopping term
this->DhopDeriv(force, X, Y, dag);
///////////////////////////////////////////////////////////
// Clover term derivative
///////////////////////////////////////////////////////////
Impl::outerProductImpl(Lambda, X, Y);
//std::cout << "Lambda:" << Lambda << std::endl;
Gamma::Algebra sigma[] = {
Gamma::Algebra::SigmaXY,
Gamma::Algebra::SigmaXZ,
Gamma::Algebra::SigmaXT,
Gamma::Algebra::MinusSigmaXY,
Gamma::Algebra::SigmaYZ,
Gamma::Algebra::SigmaYT,
Gamma::Algebra::MinusSigmaXZ,
Gamma::Algebra::MinusSigmaYZ,
Gamma::Algebra::SigmaZT,
Gamma::Algebra::MinusSigmaXT,
Gamma::Algebra::MinusSigmaYT,
Gamma::Algebra::MinusSigmaZT};
/*
sigma_{\mu \nu}=
| 0 sigma[0] sigma[1] sigma[2] |
| sigma[3] 0 sigma[4] sigma[5] |
| sigma[6] sigma[7] 0 sigma[8] |
| sigma[9] sigma[10] sigma[11] 0 |
*/
int count = 0;
clover_force = zero;
for (int mu = 0; mu < 4; mu++)
{
force_mu = zero;
for (int nu = 0; nu < 4; nu++)
{
if (mu == nu)
continue;
RealD factor;
if (nu == 4 || mu == 4)
{
factor = 2.0 * csw_t;
}
else
{
factor = 2.0 * csw_r;
}
PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok
force_mu -= factor*Cmunu(U, lambda, mu, nu); // checked
count++;
}
pokeLorentz(clover_force, U[mu] * force_mu, mu);
}
//clover_force *= csw;
force += clover_force;
}
// Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
{
conformable(lambda._grid, U[0]._grid);
GaugeLinkField out(lambda._grid), tmp(lambda._grid);
// insertion in upper staple
// please check redundancy of shift operations
// C1+
tmp = lambda * U[nu];
out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
// C2+
tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu);
out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
// C3+
tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu);
out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
// C4+
out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
// insertion in lower staple
// C1-
out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
// C2-
tmp = adj(lambda) * U[nu];
out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
// C3-
tmp = lambda * U[nu];
out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
// C4-
out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
return out;
}
private:
// here fixing the 4 dimensions, make it more general?
RealD csw_r; // Clover coefficient - spatial
RealD csw_t; // Clover coefficient - temporal
RealD diag_mass; // Mass term
CloverFieldType CloverTerm, CloverTermInv; // Clover term
CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO
CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO
CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO
CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
// eventually these can be compressed into 6x6 blocks instead of the 12x12
// using the DeGrand-Rossi basis for the gamma matrices
CloverFieldType fillCloverYZ(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
}
return T;
}
CloverFieldType fillCloverXZ(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = -F._odata[i]()();
T._odata[i]()(1, 0) = F._odata[i]()();
T._odata[i]()(2, 3) = -F._odata[i]()();
T._odata[i]()(3, 2) = F._odata[i]()();
}
return T;
}
CloverFieldType fillCloverXY(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()());
T._odata[i]()(1, 1) = timesI(F._odata[i]()());
T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 3) = timesI(F._odata[i]()());
}
return T;
}
CloverFieldType fillCloverXT(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = timesI(F._odata[i]()());
T._odata[i]()(1, 0) = timesI(F._odata[i]()());
T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()());
}
return T;
}
CloverFieldType fillCloverYT(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 1) = -(F._odata[i]()());
T._odata[i]()(1, 0) = (F._odata[i]()());
T._odata[i]()(2, 3) = (F._odata[i]()());
T._odata[i]()(3, 2) = -(F._odata[i]()());
}
return T;
}
CloverFieldType fillCloverZT(const GaugeLinkField &F)
{
CloverFieldType T(F._grid);
T = zero;
PARALLEL_FOR_LOOP
for (int i = 0; i < CloverTerm._grid->oSites(); i++)
{
T._odata[i]()(0, 0) = timesI(F._odata[i]()());
T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()());
T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());
T._odata[i]()(3, 3) = timesI(F._odata[i]()());
}
return T;
}
};
}
}
#endif // GRID_QCD_WILSON_CLOVER_FERMION_H

43
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -47,7 +47,8 @@ int WilsonFermionStatic::HandOptDslash;
template <class Impl>
WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid, RealD _mass,
const ImplParams &p)
const ImplParams &p,
const WilsonAnisotropyCoefficients &anis)
: Kernels(p),
_grid(&Fgrid),
_cbgrid(&Hgrid),
@ -60,16 +61,41 @@ WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
Umu(&Fgrid),
UmuEven(&Hgrid),
UmuOdd(&Hgrid),
_tmp(&Hgrid)
_tmp(&Hgrid),
anisotropyCoeff(anis)
{
// Allocate the required comms buffer
ImportGauge(_Umu);
if (anisotropyCoeff.isAnisotropic){
diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0);
} else {
diag_mass = 4.0 + mass;
}
}
template <class Impl>
void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
GaugeField HUmu(_Umu._grid);
HUmu = _Umu * (-0.5);
//Here multiply the anisotropy coefficients
if (anisotropyCoeff.isAnisotropic)
{
for (int mu = 0; mu < Nd; mu++)
{
GaugeLinkField U_dir = (-0.5)*PeekIndex<LorentzIndex>(_Umu, mu);
if (mu != anisotropyCoeff.t_direction)
U_dir *= (anisotropyCoeff.nu / anisotropyCoeff.xi_0);
PokeIndex<LorentzIndex>(HUmu, U_dir, mu);
}
}
else
{
HUmu = _Umu * (-0.5);
}
Impl::DoubleStore(GaugeGrid(), Umu, HUmu);
pickCheckerboard(Even, UmuEven, Umu);
pickCheckerboard(Odd, UmuOdd, Umu);
@ -83,14 +109,14 @@ template <class Impl>
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard;
Dhop(in, out, DaggerNo);
return axpy_norm(out, 4 + mass, in, out);
return axpy_norm(out, diag_mass, in, out);
}
template <class Impl>
RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard;
Dhop(in, out, DaggerYes);
return axpy_norm(out, 4 + mass, in, out);
return axpy_norm(out, diag_mass, in, out);
}
template <class Impl>
@ -114,7 +140,7 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
template <class Impl>
void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard;
typename FermionField::scalar_type scal(4.0 + mass);
typename FermionField::scalar_type scal(diag_mass);
out = scal * in;
}
@ -127,7 +153,7 @@ void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
template<class Impl>
void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard;
out = (1.0/(4.0+mass))*in;
out = (1.0/(diag_mass))*in;
}
template<class Impl>
@ -204,7 +230,7 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
FermionField Btilde(B._grid);
FermionField Atilde(B._grid);
Atilde = A;
Atilde = A;//redundant
st.HaloExchange(B, compressor);
@ -462,7 +488,6 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
{
conformable(_grid, q_in._grid);
conformable(_grid, q_out._grid);
Complex i(0.0,1.0);
PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid);
unsigned int tshift = (mu == Tp) ? 1 : 0;
unsigned int LLt = GridDefaultLatt()[Tp];

27
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -44,6 +44,21 @@ class WilsonFermionStatic {
static const int npoint = 8;
};
struct WilsonAnisotropyCoefficients: Serializable
{
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonAnisotropyCoefficients,
bool, isAnisotropic,
int, t_direction,
double, xi_0,
double, nu);
WilsonAnisotropyCoefficients():
isAnisotropic(false),
t_direction(Nd-1),
xi_0(1.0),
nu(1.0){}
};
template <class Impl>
class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
public:
@ -65,8 +80,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
// override multiply; cut number routines if pass dagger argument
// and also make interface more uniformly consistent
//////////////////////////////////////////////////////////////////
RealD M(const FermionField &in, FermionField &out);
RealD Mdag(const FermionField &in, FermionField &out);
virtual RealD M(const FermionField &in, FermionField &out);
virtual RealD Mdag(const FermionField &in, FermionField &out);
/////////////////////////////////////////////////////////
// half checkerboard operations
@ -117,8 +132,9 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
// Constructor
WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid, RealD _mass,
const ImplParams &p = ImplParams());
GridRedBlackCartesian &Hgrid, RealD _mass,
const ImplParams &p = ImplParams(),
const WilsonAnisotropyCoefficients &anis = WilsonAnisotropyCoefficients() );
// DoubleStore impl dependent
void ImportGauge(const GaugeField &_Umu);
@ -130,6 +146,7 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
// protected:
public:
RealD mass;
RealD diag_mass;
GridBase *_grid;
GridBase *_cbgrid;
@ -146,6 +163,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
LebesgueOrder Lebesgue;
LebesgueOrder LebesgueEvenOdd;
WilsonAnisotropyCoefficients anisotropyCoeff;
///////////////////////////////////////////////////////////////
// Conserved current utilities

1
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -884,7 +884,6 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
conformable(q_in._grid, FermionGrid());
conformable(q_in._grid, q_out._grid);
PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid());
Complex i(0.0, 1.0);
unsigned int tshift = (mu == Tp) ? 1 : 0;
unsigned int LLs = q_in._grid->_rdimensions[0];
unsigned int LLt = GridDefaultLatt()[Tp];

8
lib/qcd/action/fermion/WilsonKernels.h
View File

@ -55,7 +55,7 @@ template<class Impl> class WilsonKernels : public FermionOperator<Impl> , public
public:
template <bool EnableBool = true>
typename std::enable_if<Impl::Dimension == 3 && Nc == 3 &&EnableBool, void>::type
typename std::enable_if<Impl::isFundamental==true && Nc == 3 &&EnableBool, void>::type
DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1)
{
@ -99,7 +99,7 @@ public:
}
template <bool EnableBool = true>
typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool, void>::type
DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1 ) {
// no kernel choice
@ -116,7 +116,7 @@ public:
}
template <bool EnableBool = true>
typename std::enable_if<Impl::Dimension == 3 && Nc == 3 && EnableBool,void>::type
typename std::enable_if<Impl::isFundamental==true && Nc == 3 && EnableBool,void>::type
DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,
int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1)
{
@ -161,7 +161,7 @@ public:
}
template <bool EnableBool = true>
typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool,void>::type
typename std::enable_if<(Impl::isFundamental==false || (Impl::isFundamental==true && Nc != 3)) && EnableBool,void>::type
DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,SiteHalfSpinor * buf,
int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out,int interior=1,int exterior=1) {

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

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

13
lib/qcd/modules/ObservableModules.h
View File

@ -92,6 +92,19 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
PlaquetteMod(): ObsBase(NoParameters()){}
};
template < class Impl >
class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
typedef ObservableModule<PolyakovLogger<Impl>, NoParameters> ObsBase;
using ObsBase::ObsBase; // for constructors
// acquire resource
virtual void initialize(){
this->ObservablePtr.reset(new PolyakovLogger<Impl>());
}
public:
PolyakovMod(): ObsBase(NoParameters()){}
};
template < class Impl >
class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{

2
lib/qcd/observables/hmc_observable.h
View File

@ -45,5 +45,7 @@ class HmcObservable {
#include "plaquette.h"
#include "topological_charge.h"
#include "polyakov_loop.h"
#endif // HMC_OBSERVABLE_H

68
lib/qcd/observables/polyakov_loop.h Normal file
View File

@ -0,0 +1,68 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/modules/polyakov_line.h
Copyright (C) 2017
Author: David Preti <david.preti@csic.es>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#ifndef HMC_POLYAKOV_H
#define HMC_POLYAKOV_H
namespace Grid {
namespace QCD {
// this is only defined for a gauge theory
template <class Impl>
class PolyakovLogger : public HmcObservable<typename Impl::Field> {
public:
// here forces the Impl to be of gauge fields
// if not the compiler will complain
INHERIT_GIMPL_TYPES(Impl);
// necessary for HmcObservable compatibility
typedef typename Impl::Field Field;
void TrajectoryComplete(int traj,
Field &U,
GridSerialRNG &sRNG,
GridParallelRNG &pRNG) {
ComplexD polyakov = WilsonLoops<Impl>::avgPolyakovLoop(U);
int def_prec = std::cout.precision();
std::cout << GridLogMessage
<< std::setprecision(std::numeric_limits<Real>::digits10 + 1)
<< "Polyakov Loop: [ " << traj << " ] "<< polyakov << std::endl;
std::cout.precision(def_prec);
}
};
} // namespace QCD
} // namespace Grid
#endif // HMC_POLYAKOV_H

1
lib/qcd/representations/adjoint.h
View File

@ -23,6 +23,7 @@ class AdjointRep {
typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeMatrix;
typedef typename SU_Adjoint<ncolour>::LatticeAdjField LatticeField;
static const int Dimension = ncolour * ncolour - 1;
static const bool isFundamental = false;
LatticeField U;

1
lib/qcd/representations/fundamental.h
View File

@ -19,6 +19,7 @@ template <int ncolour>
class FundamentalRep {
public:
static const int Dimension = ncolour;
static const bool isFundamental = true;
// typdef to be used by the Representations class in HMC to get the
// types for the higher representation fields

1
lib/qcd/representations/two_index.h
View File

@ -29,6 +29,7 @@ class TwoIndexRep {
typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexMatrix LatticeMatrix;
typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexField LatticeField;
static const int Dimension = ncolour * (ncolour + S) / 2;
static const bool isFundamental = false;
LatticeField U;

37
lib/qcd/utils/WilsonLoops.h
View File

@ -123,6 +123,28 @@ public:
return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
}
//////////////////////////////////////////////////
// average over all x,y,z the temporal loop
//////////////////////////////////////////////////
static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4
GaugeMat Ut(Umu._grid), P(Umu._grid);
ComplexD out;
int T = Umu._grid->GlobalDimensions()[3];
int X = Umu._grid->GlobalDimensions()[0];
int Y = Umu._grid->GlobalDimensions()[1];
int Z = Umu._grid->GlobalDimensions()[2];
Ut = peekLorentz(Umu,3); //Select temporal direction
P = Ut;
for (int t=1;t<T;t++){
P = Gimpl::CovShiftForward(Ut,3,P);
}
RealD norm = 1.0/(Nc*X*Y*Z*T);
out = sum(trace(P))*norm;
return out;
}
//////////////////////////////////////////////////
// average over traced single links
//////////////////////////////////////////////////
@ -291,9 +313,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
}
}
//////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu, lower part
//////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
static void StapleLower(GaugeMat &staple, const GaugeLorentz &Umu, int mu,
int nu) {
if (nu != mu) {
@ -315,7 +337,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
//
staple = Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
mu);
}
}
@ -325,7 +349,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
static void FieldStrength(GaugeMat &FS, const GaugeLorentz &Umu, int mu, int nu){
// Fmn +--<--+ Ut +--<--+
// | | | |
// (x)+-->--+ +-->--+(x)
// (x)+-->--+ +-->--+(x) - h.c.
// | | | |
// +--<--+ +--<--+
@ -335,7 +359,9 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
GaugeMat v = Vup - Vdn;
GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu); // some redundant copies
GaugeMat vu = v*u;
FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
//FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
FS = (u*v + Cshift(vu, mu, -1));
FS = 0.125*(FS - adj(FS));
}
static Real TopologicalCharge(GaugeLorentz &U){
@ -360,6 +386,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
return TensorRemove(Tq).real();
}
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////

357
tests/core/Test_wilson_clover.cc Normal file
View File

@ -0,0 +1,357 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_wilson.cc
Copyright (C) 2015
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>
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);
GridRedBlackCartesian RBGrid(&Grid);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
std::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl;
std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl;
std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl;
std::vector<int> seeds({1, 2, 3, 4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(seeds);
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
typedef typename WilsonCloverFermionR::FermionField FermionField;
typename WilsonCloverFermionR::ImplParams params;
WilsonAnisotropyCoefficients anis;
FermionField src(&Grid);
random(pRNG, src);
FermionField result(&Grid);
result = zero;
FermionField result2(&Grid);
result2 = zero;
FermionField ref(&Grid);
ref = zero;
FermionField tmp(&Grid);
tmp = zero;
FermionField err(&Grid);
err = zero;
FermionField err2(&Grid);
err2 = zero;
FermionField phi(&Grid);
random(pRNG, phi);
FermionField chi(&Grid);
random(pRNG, chi);
LatticeGaugeField Umu(&Grid);
SU3::HotConfiguration(pRNG, Umu);
std::vector<LatticeColourMatrix> U(4, &Grid);
double volume = 1;
for (int mu = 0; mu < Nd; mu++)
{
volume = volume * latt_size[mu];
}
RealD mass = 0.1;
RealD csw_r = 1.0;
RealD csw_t = 1.0;
WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
//Dwc.ImportGauge(Umu); // not necessary, included in the constructor
std::cout << GridLogMessage << "==========================================================" << std::endl;
std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
std::cout << GridLogMessage << "==========================================================" << std::endl;
FermionField src_e(&RBGrid);
FermionField src_o(&RBGrid);
FermionField r_e(&RBGrid);
FermionField r_o(&RBGrid);
FermionField r_eo(&Grid);
pickCheckerboard(Even, src_e, src);
pickCheckerboard(Odd, src_o, src);
Dwc.Meooe(src_e, r_o);
std::cout << GridLogMessage << "Applied Meo" << std::endl;
Dwc.Meooe(src_o, r_e);
std::cout << GridLogMessage << "Applied Moe" << std::endl;
Dwc.Dhop(src, ref, DaggerNo);
setCheckerboard(r_eo, r_o);
setCheckerboard(r_eo, r_e);
err = ref - r_eo;
std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl;
std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
FermionField chi_e(&RBGrid);
FermionField chi_o(&RBGrid);
FermionField dchi_e(&RBGrid);
FermionField dchi_o(&RBGrid);
FermionField phi_e(&RBGrid);
FermionField phi_o(&RBGrid);
FermionField dphi_e(&RBGrid);
FermionField dphi_o(&RBGrid);
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
Dwc.Meooe(chi_e, dchi_o);
Dwc.Meooe(chi_o, dchi_e);
Dwc.MeooeDag(phi_e, dphi_o);
Dwc.MeooeDag(phi_o, dphi_e);
ComplexD pDce = innerProduct(phi_e, dchi_e);
ComplexD pDco = innerProduct(phi_o, dchi_o);
ComplexD cDpe = innerProduct(chi_e, dphi_e);
ComplexD cDpo = innerProduct(chi_o, dphi_o);
std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl;
std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl;
std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl;
std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
std::cout << GridLogMessage << "= Test MeeInv Mee = 1 (if csw!=0) " << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
Dwc.Mooee(chi_e, src_e);
Dwc.MooeeInv(src_e, phi_e);
Dwc.Mooee(chi_o, src_o);
Dwc.MooeeInv(src_o, phi_o);
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
err = phi - chi;
std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1 (if csw!=0) " << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
Dwc.MooeeDag(chi_e, src_e);
Dwc.MooeeInvDag(src_e, phi_e);
Dwc.MooeeDag(chi_o, src_o);
Dwc.MooeeInvDag(src_o, phi_o);
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
err = phi - chi;
std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1 (if csw!=0) " << std::endl;
std::cout << GridLogMessage << "==============================================================" << std::endl;
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
Dwc.MooeeDag(chi_e, src_e);
Dwc.MooeeInv(src_e, phi_e);
Dwc.MooeeDag(chi_o, src_o);
Dwc.MooeeInv(src_o, phi_o);
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
err = phi - chi;
std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
std::cout << GridLogMessage << "================================================================" << std::endl;
std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning " << std::endl;
std::cout << GridLogMessage << "================================================================" << std::endl;
chi = zero;
phi = zero;
tmp = zero;
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
Dwc.Mooee(src_e, chi_e);
Dwc.Mooee(src_o, chi_o);
setCheckerboard(chi, chi_e);
setCheckerboard(chi, chi_o);
setCheckerboard(src, src_e);
setCheckerboard(src, src_o);
////////////////////// Gauge Transformation
std::vector<int> seeds2({5, 6, 7, 8});
GridParallelRNG pRNG2(&Grid);
pRNG2.SeedFixedIntegers(seeds2);
LatticeColourMatrix Omega(&Grid);
LatticeColourMatrix ShiftedOmega(&Grid);
LatticeGaugeField U_prime(&Grid);
U_prime = zero;
LatticeColourMatrix U_prime_mu(&Grid);
U_prime_mu = zero;
SU<Nc>::LieRandomize(pRNG2, Omega, 1.0);
for (int mu = 0; mu < Nd; mu++)
{
U[mu] = peekLorentz(Umu, mu);
ShiftedOmega = Cshift(Omega, mu, 1);
U_prime_mu = Omega * U[mu] * adj(ShiftedOmega);
pokeLorentz(U_prime, U_prime_mu, mu);
}
/////////////////
WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
Dwc_prime.ImportGauge(U_prime);
tmp = Omega * src;
pickCheckerboard(Even, src_e, tmp);
pickCheckerboard(Odd, src_o, tmp);
Dwc_prime.Mooee(src_e, phi_e);
Dwc_prime.Mooee(src_o, phi_o);
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
err = chi - adj(Omega) * phi;
std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
std::cout << GridLogMessage << "=================================================================" << std::endl;
std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning " << std::endl;
std::cout << GridLogMessage << "================================================================" << std::endl;
chi = zero;
phi = zero;
WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
Dw.ImportGauge(Umu);
Dw.M(src, result);
Dwc.M(src, chi);
Dwc_prime.M(Omega * src, phi);
WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
Dw_prime.ImportGauge(U_prime);
Dw_prime.M(Omega * src, result2);
err = chi - adj(Omega) * phi;
err2 = result - adj(Omega) * result2;
std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl;
std::cout << GridLogMessage << "norm diff WilsonClover " << norm2(err2) << std::endl;
std::cout << GridLogMessage << "==========================================================" << std::endl;
std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson " << std::endl;
std::cout << GridLogMessage << "==========================================================" << std::endl;
chi = zero;
phi = zero;
err = zero;
WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0
Dwc_csw0.ImportGauge(Umu);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
Dw.Mooee(src_e, chi_e);
Dw.Mooee(src_o, chi_o);
Dwc_csw0.Mooee(src_e, phi_e);
Dwc_csw0.Mooee(src_o, phi_o);
setCheckerboard(chi, chi_e);
setCheckerboard(chi, chi_o);
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
setCheckerboard(src, src_e);
setCheckerboard(src, src_o);
err = chi - phi;
std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
std::cout << GridLogMessage << "==========================================================" << std::endl;
std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec " << std::endl;
std::cout << GridLogMessage << "==========================================================" << std::endl;
chi = zero;
phi = zero;
err = zero;
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
// M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o)
Dwc.M(src, ref); // Reference result from the unpreconditioned operator
// EO matrix
Dwc.Mooee(src_e, chi_e);
Dwc.Mooee(src_o, chi_o);
Dwc.Meooe(src_o, phi_e);
Dwc.Meooe(src_e, phi_o);
phi_o += chi_o;
phi_e += chi_e;
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
err = ref - phi;
std::cout << GridLogMessage << "ref (unpreconditioned operator) diff :" << norm2(ref) << std::endl;
std::cout << GridLogMessage << "phi (EO decomposition) diff :" << norm2(phi) << std::endl;
std::cout << GridLogMessage << "norm diff :" << norm2(err) << std::endl;
Grid_finalize();
}

7
tests/forces/Test_wilson_force.cc
View File

@ -50,7 +50,12 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
std::vector<int> vrand(4);
std::srand(std::time(0));
std::generate(vrand.begin(), vrand.end(), std::rand);
std::cout << GridLogMessage << vrand << std::endl;
pRNG.SeedFixedIntegers(vrand);
//pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
LatticeFermion phi (&Grid); gaussian(pRNG,phi);
LatticeFermion Mphi (&Grid);

194
tests/forces/Test_wilsonclover_force.cc Normal file
View File

@ -0,0 +1,194 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_force.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);
GridRedBlackCartesian RBGrid(&Grid);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
std::vector<int> seeds({1, 2, 30, 50});
GridParallelRNG pRNG(&Grid);
std::vector<int> vrand(4);
std::srand(std::time(0));
std::generate(vrand.begin(), vrand.end(), std::rand);
std::cout << GridLogMessage << vrand << std::endl;
pRNG.SeedFixedIntegers(vrand);
//pRNG.SeedFixedIntegers(seeds);
LatticeFermion phi(&Grid);
gaussian(pRNG, phi);
LatticeFermion Mphi(&Grid);
LatticeFermion MphiPrime(&Grid);
LatticeGaugeField U(&Grid);
std::vector<int> site = {0, 0, 0, 0};
SU3::HotConfiguration(pRNG, U);
//SU3::ColdConfiguration(pRNG, U);// Clover term zero
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD mass = 0.1;
Real csw = 1.0;
WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw, csw);
Dw.ImportGauge(U);
Dw.M(phi, Mphi);
ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(&Grid);
LatticeGaugeField tmp(&Grid);
////////////////////////////////////////////
Dw.MDeriv(tmp, Mphi, phi, DaggerNo);
UdSdU = tmp;
Dw.MDeriv(tmp, phi, Mphi, DaggerYes);
UdSdU += tmp;
/////////////////////////////////////////////
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.00005;
RealD Hmom = 0.0;
RealD Hmomprime = 0.0;
RealD Hmompp = 0.0;
LatticeColourMatrix mommu(&Grid);
LatticeColourMatrix forcemu(&Grid);
LatticeGaugeField mom(&Grid);
LatticeGaugeField Uprime(&Grid);
for (int mu = 0; mu < Nd; mu++)
{
// Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
Hmom -= real(sum(trace(mommu * mommu)));
PokeIndex<LorentzIndex>(mom, mommu, mu);
parallel_for(int ss = 0; ss < mom._grid->oSites(); ss++)
{
Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]);
}
}
std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
// New action
Dw.ImportGauge(Uprime);
Dw.M(phi, MphiPrime);
ComplexD Sprime = innerProduct(MphiPrime, MphiPrime);
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(&Grid);
dS = zero;
LatticeComplex dSmom(&Grid);
dSmom = zero;
LatticeComplex dSmom2(&Grid);
dSmom2 = zero;
for (int mu = 0; mu < Nd; mu++)
{
mommu = PeekIndex<LorentzIndex>(UdSdU, mu); // P_mu =
mommu = Ta(mommu) * 2.0; // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag)
PokeIndex<LorentzIndex>(UdSdU, mommu, mu); // UdSdU_mu = Mom
}
std::cout << GridLogMessage << "Antihermiticity tests" << std::endl;
for (int mu = 0; mu < Nd; mu++)
{
mommu = PeekIndex<LorentzIndex>(mom, mu);
std::cout << GridLogMessage << " Mommu " << norm2(mommu) << std::endl;
mommu = mommu + adj(mommu);
std::cout << GridLogMessage << " Mommu + Mommudag " << norm2(mommu) << std::endl;
mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
std::cout << GridLogMessage << " dsdumu " << norm2(mommu) << std::endl;
mommu = mommu + adj(mommu);
std::cout << GridLogMessage << " dsdumu + dag " << norm2(mommu) << std::endl;
std::cout << "" << std::endl;
}
/////////////////////////////////////////////////////
for (int mu = 0; mu < Nd; mu++)
{
forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);
mommu = PeekIndex<LorentzIndex>(mom, mu);
// Update PF action density
dS = dS + trace(mommu * forcemu) * dt;
dSmom = dSmom - trace(mommu * forcemu) * dt;
dSmom2 = dSmom2 - trace(forcemu * forcemu) * (0.25 * dt * dt);
// Update mom action density
mommu = mommu + forcemu * (dt * 0.5);
Hmomprime -= real(sum(trace(mommu * mommu)));
}
ComplexD dSpred = sum(dS);
ComplexD dSm = sum(dSmom);
ComplexD dSm2 = sum(dSmom2);
std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
std::cout << GridLogMessage << "Final mom hamiltonian is " << Hmomprime << std::endl;
std::cout << GridLogMessage << "Delta mom hamiltonian is " << Hmomprime - Hmom << std::endl;
std::cout << GridLogMessage << " S " << S << std::endl;
std::cout << GridLogMessage << " Sprime " << Sprime << std::endl;
std::cout << GridLogMessage << "dS (S' - S) :" << Sprime - S << std::endl;
std::cout << GridLogMessage << "predict dS (force) :" << dSpred << std::endl;
std::cout << GridLogMessage << "dSm " << dSm << std::endl;
std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl;
std::cout << GridLogMessage << "Total dS " << Hmomprime - Hmom + Sprime - S << std::endl;
assert(fabs(real(Sprime - S - dSpred)) < 1.0);
std::cout << GridLogMessage << "Done" << std::endl;
Grid_finalize();
}

168
tests/hadrons/Test_hadrons_2AS_spectrum.cc Normal file
View File

@ -0,0 +1,168 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons_spectrum.cc
Copyright (C) 2015
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.
*******************************************************************************/
#include <Grid/Hadrons/Application.hpp>
using namespace Grid;
using namespace Hadrons;
BEGIN_HADRONS_NAMESPACE
BEGIN_MODULE_NAMESPACE(MFermion)
MODULE_REGISTER_NS(GaugeProp2AS, TGaugeProp<WilsonTwoIndexAntiSymmetricImplR>, MFermion);
END_MODULE_NAMESPACE
BEGIN_MODULE_NAMESPACE(MSource)
MODULE_REGISTER_NS(Point2AS, TPoint<WilsonTwoIndexAntiSymmetricImplR>, MSource);
END_MODULE_NAMESPACE
BEGIN_MODULE_NAMESPACE(MContraction)
MODULE_REGISTER_NS(Meson2AS, ARG(TMeson<WilsonTwoIndexAntiSymmetricImplR, WilsonTwoIndexAntiSymmetricImplR>), MContraction);
// MODULE_REGISTER_NS(BaryonMultirep, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
END_MODULE_NAMESPACE
BEGIN_MODULE_NAMESPACE(MSink)
MODULE_REGISTER_NS(ScalarPoint2AS, TPoint<WilsonTwoIndexAntiSymmetricImplR>, MSink);
END_MODULE_NAMESPACE
BEGIN_MODULE_NAMESPACE(MSolver)
MODULE_REGISTER_NS(RBPrecCG2AS, TRBPrecCG<WilsonTwoIndexAntiSymmetricImplR>, MSolver);
END_MODULE_NAMESPACE
BEGIN_MODULE_NAMESPACE(MAction)
MODULE_REGISTER_NS(WilsonClover2AS, TWilsonClover<WilsonTwoIndexAntiSymmetricImplR>, MAction);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
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"};
std::vector<double> mass = {-0.01, -0.04};
double csw = 1.0;
// global parameters
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
application.setPar(globalPar);
// gauge field
application.createModule<MGauge::Unit>("gauge");
MSource::Point2AS::Par ptPar;
ptPar.position = "0 0 0 0";
application.createModule<MSource::Point2AS>("pt", ptPar);
// sink
MSink::ScalarPoint2AS::Par sinkPar;
sinkPar.mom = "0 0 0";
application.createModule<MSink::ScalarPoint2AS>("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::WilsonClover2AS::Par actionPar;
actionPar.gauge = "gauge";
actionPar.mass = mass[i];
actionPar.csw_r = csw;
actionPar.csw_t = csw;
actionPar.clover_anisotropy.isAnisotropic= false;
actionPar.clover_anisotropy.t_direction = Nd-1 ;
actionPar.clover_anisotropy.xi_0 = 1.0 ;
actionPar.clover_anisotropy.nu = 1.0 ;
actionPar.boundary = boundary;
application.createModule<MAction::WilsonClover2AS>("WilsonClover2AS_" + flavour[i], actionPar);
// solvers
MSolver::RBPrecCG2AS::Par solverPar;
solverPar.action = "WilsonClover2AS_" + flavour[i];
solverPar.residual = 1.0e-8;
application.createModule<MSolver::RBPrecCG2AS>("CG_" + flavour[i],
solverPar);
// propagators
MFermion::GaugeProp2AS::Par quarkPar;
quarkPar.solver = "CG_" + flavour[i];
quarkPar.source = "pt";
application.createModule<MFermion::GaugeProp2AS>("Qpt_" + flavour[i], quarkPar);
quarkPar.source = "z2";
application.createModule<MFermion::GaugeProp2AS>("QZ2_" + flavour[i], quarkPar);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
MContraction::Meson2AS::Par mesPar;
mesPar.output = "mesons2AS/pt_" + flavour[i] + flavour[j];
mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "all";
mesPar.sink = "sink";
application.createModule<MContraction::Meson2AS>("meson_pt_"
+ flavour[i] + flavour[j],
mesPar);
// mesPar.output = "mesons2AS/Z2_" + flavour[i] + flavour[j];
// mesPar.q1 = "QZ2_" + flavour[i];
// mesPar.q2 = "QZ2_" + flavour[j];
// mesPar.gammas = "all";
// mesPar.sink = "sink";
// application.createModule<MContraction::Meson2AS>("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];
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;
}

4
tests/hadrons/Test_hadrons_quark.cc
View File

@ -130,8 +130,8 @@ int main(int argc, char **argv)
for (int c = 0; c < Nc; ++c)
{
ref = prop;
PropToFerm(ferm, prop, s, c);
FermToProp(prop, ferm, s, c);
PropToFerm<WilsonImplR>(ferm, prop, s, c);
FermToProp<WilsonImplR>(prop, ferm, s, c);
std::cout << "Spin = " << s << ", Colour = " << c << std::endl;
ref -= prop;

157
tests/hadrons/Test_hadrons_wilsonFund.cc Normal file
View File

@ -0,0 +1,157 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons_spectrum.cc
Copyright (C) 2015
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.
*******************************************************************************/
#include <Grid/Hadrons/Application.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.seed = "1 2 3 4";
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;
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;
}

139
tests/hmc/Test_hmc_EOWilsonCloverFermionGauge.cc Normal file
View File

@ -0,0 +1,139 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
Copyright (C) 2016
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#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 GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplR FermionImplPolicy;
typedef WilsonCloverFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Possibile to create the module by hand
// hardcoding parameters or using a Reader
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 5.6 ;
WilsonGaugeActionR Waction(beta);
// temporarily need a gauge field
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
LatticeGaugeField U(GridPtr);
Real mass = 0.01;
Real csw = 1.0;
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw);
ConjugateGradient<FermionField> CG(1.0e-8, 2000);
TwoFlavourEvenOddPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
// Set smearing (true/false), default: false
Nf2.is_smeared = false;
// Collect actions
ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
/*
double rho = 0.1; // smearing parameter
int Nsmear = 2; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
UGrid, Nsmear, Stout);
*/
// HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main

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

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

212
tests/hmc/Test_hmc_WC2SFG_Production.cc Normal file
View File

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

210
tests/hmc/Test_hmc_WCFG_Production.cc Normal file
View File

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

224
tests/hmc/Test_hmc_WCMixedRepFG_Production.cc Normal file
View File

@ -0,0 +1,224 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
Copyright (C) 2015
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp>
Author: paboyle <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"
namespace Grid{
struct FermionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
double, mass,
double, csw,
double, StoppingCondition,
int, MaxCGIterations,
bool, ApplySmearing);
};
struct WilsonCloverHMCParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonCloverHMCParameters,
double, gauge_beta,
FermionParameters, WilsonCloverFund,
FermionParameters, WilsonCloverAS)
template <class ReaderClass >
WilsonCloverHMCParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
struct SmearingParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
double, rho,
Integer, Nsmear)
template <class ReaderClass >
SmearingParameters(Reader<ReaderClass>& Reader){
read(Reader, "StoutSmearing", *this);
}
};
}
int main(int argc, char **argv) {
using namespace Grid;
using namespace Grid::QCD;
// Here change the allowed (higher) representations
typedef Representations< FundamentalRepresentation, TwoIndexAntiSymmetricRepresentation> TheRepresentations;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef GenericHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
typedef WilsonImplR FundImplPolicy;
typedef WilsonCloverFermionR FundFermionAction;
typedef typename FundFermionAction::FermionField FundFermionField;
typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy;
typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction;
typedef typename ASymmFermionAction::FermionField ASymmFermionField;
typedef Grid::JSONReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.ReadCommandLine(argc, argv);
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
WilsonCloverHMCParameters MyParams(Reader);
// Apply smearing to the fermionic action
bool ApplySmearingFund = MyParams.WilsonCloverFund.ApplySmearing;
bool ApplySmearingAS = MyParams.WilsonCloverAS.ApplySmearing;
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams(Reader);
/*
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
*/
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar(Reader);
/*
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
*/
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
TopologyObsParameters TopParams(Reader);
TheHMC.Resources.AddObservable<QObs>(TopParams);
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
//RealD beta = 5.6;
WilsonGaugeActionR Waction(MyParams.gauge_beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
FundamentalRepresentation::LatticeField UF(GridPtr);
TwoIndexAntiSymmetricRepresentation::LatticeField UAS(GridPtr);
Real Fundmass = MyParams.WilsonCloverFund.mass;
Real Fundcsw = MyParams.WilsonCloverFund.csw;
Real ASmass = MyParams.WilsonCloverAS.mass;
Real AScsw = MyParams.WilsonCloverAS.csw;
std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl;
std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl;
FundFermionAction FundFermOp(UF, *GridPtr, *GridRBPtr, Fundmass, Fundcsw, Fundcsw);
ConjugateGradient<FundFermionField> CG_Fund(MyParams.WilsonCloverFund.StoppingCondition, MyParams.WilsonCloverFund.MaxCGIterations);
TwoFlavourPseudoFermionAction<FundImplPolicy> Nf2_Fund(FundFermOp, CG_Fund, CG_Fund);
ASymmFermionAction ASFermOp(UAS, *GridPtr, *GridRBPtr, ASmass, AScsw, AScsw);
ConjugateGradient<ASymmFermionField> CG_AS(MyParams.WilsonCloverAS.StoppingCondition, MyParams.WilsonCloverAS.MaxCGIterations);
TwoFlavourPseudoFermionAction<ASymmImplPolicy> Nf2_AS(ASFermOp, CG_AS, CG_AS);
Nf2_Fund.is_smeared = ApplySmearingFund;
Nf2_AS.is_smeared = ApplySmearingAS;
// Collect actions
ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1);
Level1.push_back(&Nf2_Fund);
Level1.push_back(&Nf2_AS);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
TheHMC.Parameters.initialize(Reader);
//TheHMC.Parameters.MD.MDsteps = 20;
//TheHMC.Parameters.MD.trajL = 1.0;
/*
if (ApplySmearingFund || ApplySmearingAS){
SmearingParameters SmPar(Reader);
//double rho = 0.1; // smearing parameter
//int Nsmear = 3; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
TheHMC.Run(SmearingPolicy); // for smearing
} else {
TheHMC.Run(); // no smearing
}
*/
TheHMC.Run();
//TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
//TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main

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

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

114
tests/hmc/Test_hmc_WG_Production.cc Normal file
View File

@ -0,0 +1,114 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
Copyright (C) 2015
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
namespace Grid{
struct ActionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(ActionParameters,
double, beta)
template <class ReaderClass >
ActionParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
}
int main(int argc, char **argv) {
using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv);
GridLogLayout();
// Typedefs to simplify notation
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
HMCWrapper TheHMC;
typedef Grid::JSONReader Serialiser;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
// Reader, file should come from command line
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
// Checkpointer definition
CheckpointerParameters CPparams(Reader);
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar(Reader);
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
TheHMC.Resources.AddObservable<PlaqObs>();
TopologyObsParameters TopParams(Reader);
TheHMC.Resources.AddObservable<QObs>(TopParams);
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
ActionParameters WilsonPar(Reader);
//RealD beta = 6.4 ;
WilsonGaugeActionR Waction(WilsonPar.beta);
ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&Waction);
//Level1.push_back(WGMod.getPtr());
TheHMC.TheAction.push_back(Level1);
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
TheHMC.Parameters.initialize(Reader);
//TheHMC.Parameters.MD.MDsteps = 17;
//TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

129
tests/hmc/Test_hmc_WilsonCloverFermionGauge.cc Normal file
View File

@ -0,0 +1,129 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
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 GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplR FermionImplPolicy;
typedef WilsonCloverFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 5.6;
WilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
Real mass = 0.01;
Real csw = 1.0;
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, csw);
ConjugateGradient<FermionField> CG(1.0e-8, 5000);
TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
// Set smearing (true/false), default: false
Nf2.is_smeared = false;
// Collect actions
ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
/*
double rho = 0.1; // smearing parameter
int Nsmear = 2; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
UGrid, Nsmear, Stout);
*/
// HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run(); // no smearing
// TheHMC.Run(SmearingPolicy); // for smearing
Grid_finalize();
} // main

178
tests/lanczos/Test_WCMultiRep_lanczos.cc Normal file
View File

@ -0,0 +1,178 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_lanczos.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;
//typedef WilsonCloverFermionR FermionOp;
//typedef typename WilsonFermionR::FermionField FermionField;
typedef WilsonImplR FundImplPolicy;
typedef WilsonCloverFermionR FundFermionAction;
typedef typename FundFermionAction::FermionField FundFermionField;
typedef WilsonTwoIndexAntiSymmetricImplR ASymmImplPolicy;
typedef WilsonCloverTwoIndexAntiSymmetricFermionR ASymmFermionAction;
typedef typename ASymmFermionAction::FermionField ASymmFermionField;
RealD AllZero(RealD x) { return 0.; }
int main(int argc, char** argv) {
Grid_init(&argc, &argv);
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
GridDefaultMpi());
GridRedBlackCartesian* UrbGrid =
SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian* FGrid = UGrid;
GridRedBlackCartesian* FrbGrid = UrbGrid;
printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid,
FrbGrid);
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5rb(FrbGrid);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG pRNG(UGrid);
GridSerialRNG sRNG;
FundamentalRepresentation::LatticeField Umu(UGrid);
TwoIndexAntiSymmetricRepresentation HiRep(UGrid);
TwoIndexAntiSymmetricRepresentation::LatticeField UmuAS(UGrid);
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.format = "IEEE64BIG";
//NerscHmcCheckpointer<PeriodicGimplR> Checkpoint(std::string("ckpoint_lat"),
// std::string("ckpoint_rng"), 1);
NerscHmcCheckpointer<PeriodicGimplR> Checkpoint(CPparams);
int CNFGSTART=1;
int CNFGEND=2;
int CNFGSTEP=1;
Real Fundmass = -0.1;
Real Fundcsw = 1.0;
Real ASmass = -0.1;
Real AScsw = 1.0;
std::cout << "Fund: mass and csw" << Fundmass << " and " << Fundcsw << std::endl;
std::cout << "AS : mass and csw" << ASmass << " and " << AScsw << std::endl;
const int Nstop = 30;
const int Nk = 40;
const int Np = 40;
const int Nm = Nk + Np;
const int MaxIt = 10000;
RealD resid = 1.0e-8;
for (int cnfg=CNFGSTART;cnfg<=CNFGEND;cnfg+=CNFGSTEP){
Checkpoint.CheckpointRestore(cnfg,Umu, sRNG, pRNG);
//SU4::HotConfiguration(RNG4, Umu); // temporary, then read.
HiRep.update_representation(Umu);
UmuAS = HiRep.U;
FundFermionAction FundFermOp(Umu,*FGrid,*FrbGrid, Fundmass, Fundcsw, Fundcsw);
MdagMLinearOperator<FundFermionAction,FundFermionField> HermOpFund(FundFermOp); /// <-----
ASymmFermionAction ASFermOp(UmuAS,*FGrid,*FrbGrid, ASmass, AScsw, AScsw);
MdagMLinearOperator<ASymmFermionAction,ASymmFermionField> HermOpAS(ASFermOp); /// <-----
std::vector<double> Coeffs{0, -1.};
Polynomial<FundFermionField> FundPolyX(Coeffs);
//Chebyshev<FundFermionField> FundCheb(0.0, 10., 12);
FunctionHermOp<FundFermionField> FundPolyXOp(FundPolyX,HermOpFund);
PlainHermOp<FundFermionField> FundOp (HermOpFund);
ImplicitlyRestartedLanczos<FundFermionField> IRL_Fund(FundOp, FundPolyXOp, Nstop, Nk, Nm,
resid, MaxIt);
Polynomial<ASymmFermionField> ASPolyX(Coeffs);
//Chebyshev<ASymmFermionField> ASCheb(0.0, 10., 12);
FunctionHermOp<ASymmFermionField> ASPolyXOp(ASPolyX,HermOpAS);
PlainHermOp<ASymmFermionField> ASOp (HermOpAS);
ImplicitlyRestartedLanczos<ASymmFermionField> IRL_AS(ASOp, ASPolyXOp, Nstop, Nk, Nm,
resid, MaxIt);
std::vector<RealD> Fundeval(Nm);
std::vector<RealD> ASeval(Nm);
FundFermionField Fundsrc(FGrid);
ASymmFermionField ASsrc(FGrid);
gaussian(RNG5, Fundsrc);
gaussian(RNG5, ASsrc);
std::vector<FundFermionField> Fundevec(Nm, FGrid);
std::vector<ASymmFermionField> ASevec(Nm, FGrid);
for (int i = 0; i < 1; i++) {
std::cout << i << " / " << Nm << "Fund: grid pointer " << Fundevec[i]._grid
<< std::endl;
};
for (int i = 0; i < 1; i++) {
std::cout << i << " / " << Nm << "AS: grid pointer " << ASevec[i]._grid
<< std::endl;
};
int FundNconv, ASNconv;
IRL_Fund.calc(Fundeval, Fundevec, Fundsrc, FundNconv);
IRL_AS.calc(ASeval, ASevec, ASsrc, ASNconv);
for (int i=0;i<FundNconv;i++){
std::cout << "Fund: eval[" << i << "] = " << Fundeval[i] << std::endl;
}
for (int i=0;i<ASNconv;i++){
std::cout << "2Index: eval[" << i << "] = " << ASeval[i] << std::endl;
}
}
Grid_finalize();
}

3
tests/qdpxx/Makefile.am
View File

@ -1,4 +1,5 @@
AM_CXXFLAGS += `chroma-config --cxxflags`
AM_LDFLAGS += `chroma-config --ldflags` `chroma-config --libs`
AM_LDFLAGS += `chroma-config --ldflags`
LIBS += `chroma-config --libs`
include Make.inc

3
tests/qdpxx/Makefile.am.qdpxx
View File

@ -1,6 +1,7 @@
# additional include paths necessary to compile the C++ library
AM_CXXFLAGS = -I$(top_srcdir)/include `chroma-config --cxxflags`
AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags` `chroma-config --libs`
AM_LDFLAGS = -L$(top_builddir)/lib `chroma-config --ldflags`
AM_LIBS = `chroma-config --libs`
include Make.inc

519
tests/qdpxx/Test_qdpxx_wilson.cc Normal file
View File

@ -0,0 +1,519 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc
Copyright (C) 2017
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
// Mass
double mq = 0.1;
// Define Wilson Types
typedef Grid::QCD::WilsonImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField;
enum ChromaAction
{
Wilson, // Wilson
WilsonClover // CloverFermions
};
namespace Chroma
{
class ChromaWrapper
{
public:
typedef multi1d<LatticeColorMatrix> U;
typedef LatticeFermion T4;
static void ImportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(LCM, gr, x);
for (int mu = 0; mu < 4; mu++)
{
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
cc = LCM(mu)()(i, j);
c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc)));
QDP::pokeColor(cm, c, i, j);
}
}
QDP::pokeSite(ch[mu], cm, cx);
}
}
}
}
}
}
static void ExportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
for (int mu = 0; mu < 4; mu++)
{
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
cm = QDP::peekSite(ch[mu], cx);
c = QDP::peekColor(cm, i, j);
cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c)));
LCM(mu)
()(i, j) = cc;
}
}
}
Grid::pokeSite(LCM, gr, x);
}
}
}
}
}
// Specific for Wilson Fermions
static void ImportFermion(Grid::QCD::LatticeFermion &gr,
QDP::LatticeFermion &ch)
{
Grid::QCD::SpinColourVector F;
Grid::Complex c;
QDP::Fermion cF;
QDP::SpinVector cS;
QDP::Complex cc;
std::vector<int> x(4); // explicit 4d fermions in Grid
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(F, gr, x);
for (int j = 0; j < 3; j++)
{
for (int sp = 0; sp < 4; sp++)
{
c = F()(sp)(j);
cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c)));
QDP::pokeSpin(cS, cc, sp);
}
QDP::pokeColor(cF, cS, j);
}
QDP::pokeSite(ch, cF, cx);
}
}
}
}
}
// Specific for 4d Wilson fermions
static void ExportFermion(Grid::QCD::LatticeFermion &gr,
QDP::LatticeFermion &ch)
{
Grid::QCD::SpinColourVector F;
Grid::Complex c;
QDP::Fermion cF;
QDP::SpinVector cS;
QDP::Complex cc;
std::vector<int> x(4); // 4d fermions
QDP::multi1d<int> cx(4);
std::vector<int> gd = gr._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
cF = QDP::peekSite(ch, cx);
for (int sp = 0; sp < 4; sp++)
{
for (int j = 0; j < 3; j++)
{
cS = QDP::peekColor(cF, j);
cc = QDP::peekSpin(cS, sp);
c = Grid::Complex(QDP::toDouble(QDP::real(cc)),
QDP::toDouble(QDP::imag(cc)));
F()
(sp)(j) = c;
}
}
Grid::pokeSite(F, gr, x);
}
}
}
}
}
static Handle<Chroma::UnprecLinearOperator<T4, U, U>> GetLinOp(U &u, ChromaAction params)
{
QDP::Real _mq(mq);
QDP::multi1d<int> bcs(QDP::Nd);
// Boundary conditions
bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;
if (params == Wilson)
{
Chroma::WilsonFermActParams p;
p.Mass = _mq;
AnisoParam_t _apar;
_apar.anisoP = true;
_apar.t_dir = 3; // in 4d
_apar.xi_0 = 2.0;
_apar.nu = 1.0;
p.anisoParam = _apar;
Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
Chroma::UnprecWilsonFermAct S_f(cfs, p);
Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
return S_f.linOp(ffs);
}
if (params == WilsonClover)
{
Chroma::CloverFermActParams p;
p.Mass = _mq;
p.clovCoeffR = QDP::Real(1.0);
p.clovCoeffT = QDP::Real(2.0);
p.u0 = QDP::Real(1.0);
AnisoParam_t _apar;
_apar.anisoP = true;
_apar.t_dir = 3; // in 4d
_apar.xi_0 = 2.0;
_apar.nu = 1.0;
p.anisoParam = _apar;
Chroma::Handle<Chroma::FermBC<T4, U, U>> fbc(new Chroma::SimpleFermBC<T4, U, U>(bcs));
Chroma::Handle<Chroma::CreateFermState<T4, U, U>> cfs(new Chroma::CreateSimpleFermState<T4, U, U>(fbc));
Chroma::UnprecCloverFermAct S_f(cfs, p);
Chroma::Handle<Chroma::FermState<T4, U, U>> ffs(S_f.createState(u));
return S_f.linOp(ffs);
}
}
};
} // namespace Chroma
void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, FermionField &res, int dag)
{
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
Chroma::ChromaWrapper::ImportGauge(lat, u);
QDP::LatticeFermion check;
QDP::LatticeFermion result;
QDP::LatticeFermion psi;
Chroma::ChromaWrapper::ImportFermion(src, psi);
for (int mu = 0; mu < 4; mu++)
{
std::cout << "Imported Gauge norm [" << mu << "] " << QDP::norm2(u[mu]) << std::endl;
}
std::cout << "Imported Fermion norm " << QDP::norm2(psi) << std::endl;
typedef QDP::LatticeFermion T;
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
auto linop = Chroma::ChromaWrapper::GetLinOp(u, action);
printf("Calling Chroma Linop\n");
fflush(stdout);
if (dag)
(*linop)(check, psi, Chroma::MINUS);
else
(*linop)(check, psi, Chroma::PLUS);
printf("Called Chroma Linop\n");
fflush(stdout);
// std::cout << "Calling Chroma Linop " << std::endl;
// linop->evenEvenLinOp(tmp, psi, isign);
// check[rb[0]] = tmp;
// linop->oddOddLinOp(tmp, psi, isign);
// check[rb[1]] = tmp;
// linop->evenOddLinOp(tmp, psi, isign);
// check[rb[0]] += tmp;
// linop->oddEvenLinOp(tmp, psi, isign);
// check[rb[1]] += tmp;
Chroma::ChromaWrapper::ExportFermion(res, check);
}
void make_gauge(GaugeField &Umu, FermionField &src)
{
using namespace Grid;
using namespace Grid::QCD;
std::vector<int> seeds4({1, 2, 3, 4});
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4, Umu);
// Fermion field
Grid::gaussian(RNG4, src);
/*
Grid::QCD::SpinColourVector F;
Grid::Complex c;
std::vector<int> x(4); // 4d fermions
std::vector<int> gd = src._grid->GlobalDimensions();
for (x[0] = 0; x[0] < gd[0]; x[0]++)
{
for (x[1] = 0; x[1] < gd[1]; x[1]++)
{
for (x[2] = 0; x[2] < gd[2]; x[2]++)
{
for (x[3] = 0; x[3] < gd[3]; x[3]++)
{
for (int sp = 0; sp < 4; sp++)
{
for (int j = 0; j < 3; j++) // colours
{
F()(sp)(j) = Grid::Complex(0.0,0.0);
if (((sp == 0)|| (sp==3)) && (j==2))
{
c = Grid::Complex(1.0, 0.0);
F()(sp)(j) = c;
}
}
}
Grid::pokeSite(F, src, x);
}
}
}
}
*/
}
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)
{
using namespace Grid;
using namespace Grid::QCD;
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::RealD _mass = mq;
if (action == Wilson)
{
WilsonAnisotropyCoefficients anis;
anis.isAnisotropic = true;
anis.t_direction = 3;
anis.xi_0 = 2.0;
anis.nu = 1.0;
WilsonImplParams iParam;
Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl;
if (dag)
Wf.Mdag(src, res);
else
Wf.M(src, res);
return;
}
if (action == WilsonClover)
{
Grid::RealD _csw_r = 1.0;
Grid::RealD _csw_t = 2.0;
WilsonAnisotropyCoefficients anis;
anis.isAnisotropic = true;
anis.t_direction = 3;
anis.xi_0 = 2.0;
anis.nu = 1.0;
WilsonImplParams CloverImplParam;
Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam);
Wf.ImportGauge(Umu);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;
if (dag)
Wf.Mdag(src, res);
else
Wf.M(src, res);
return;
}
assert(0);
}
int main(int argc, char **argv)
{
/********************************************************
* Setup QDP
*********************************************************/
Chroma::initialize(&argc, &argv);
Chroma::WilsonTypeFermActs4DEnv::registerAll();
/********************************************************
* Setup Grid
*********************************************************/
Grid::Grid_init(&argc, &argv);
Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd);
for (int mu = 0; mu < 4; mu++)
nrow[mu] = gd[mu];
QDP::Layout::setLattSize(nrow);
QDP::Layout::create();
GaugeField Ug(UGrid);
FermionField src(UGrid);
FermionField res_chroma(UGrid);
FermionField res_grid(UGrid);
FermionField only_wilson(UGrid);
FermionField difference(UGrid);
std::vector<ChromaAction> ActionList({Wilson, WilsonClover});
std::vector<std::string> ActionName({"Wilson", "WilsonClover"});
{
for (int i = 0; i < ActionList.size(); i++)
{
std::cout << "*****************************" << std::endl;
std::cout << "Action " << ActionName[i] << std::endl;
std::cout << "*****************************" << std::endl;
make_gauge(Ug, src); // fills the gauge field and the fermion field with random numbers
for (int dag = 0; dag < 2; dag++)
{
{
std::cout << "Dag = " << dag << std::endl;
calc_chroma(ActionList[i], Ug, src, res_chroma, dag);
// Remove the normalisation of Chroma Gauge links ????????
std::cout << "Norm of Chroma " << ActionName[i] << " multiply " << Grid::norm2(res_chroma) << std::endl;
calc_grid(ActionList[i], Ug, src, res_grid, dag);
std::cout << "Norm of gauge " << Grid::norm2(Ug) << std::endl;
std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl;
difference = res_chroma - res_grid;
std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl;
}
}
std::cout << "Finished test " << std::endl;
Chroma::finalize();
}
}
}