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

Merge branch 'feature/hadrons' of https://github.com/paboyle/Grid into feature/hadrons

Browse Source
This commit is contained in:
fionnoh 2018-06-28 16:13:07 +01:00
commit 8fe9a13cdd
28 changed files with 1627 additions and 264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
extras/Hadrons/Global.hpp
View File

@ -93,17 +93,15 @@ typedef typename SImpl::Field ScalarField##suffix;\
typedef typename SImpl::Field PropagatorField##suffix; typedef typename SImpl::Field PropagatorField##suffix;
#define SOLVER_TYPE_ALIASES(FImpl, suffix)\ #define SOLVER_TYPE_ALIASES(FImpl, suffix)\
typedef std::function<void(FermionField##suffix &,\ typedef Solver<FImpl> Solver##suffix;
const FermionField##suffix &)> SolverFn##suffix;
#define SINK_TYPE_ALIASES(suffix)\ #define SINK_TYPE_ALIASES(suffix)\
typedef std::function<SlicedPropagator##suffix\ typedef std::function<SlicedPropagator##suffix\
(const PropagatorField##suffix &)> SinkFn##suffix; (const PropagatorField##suffix &)> SinkFn##suffix;
#define FGS_TYPE_ALIASES(FImpl, suffix)\ #define FG_TYPE_ALIASES(FImpl, suffix)\
FERM_TYPE_ALIASES(FImpl, suffix)\ FERM_TYPE_ALIASES(FImpl, suffix)\
GAUGE_TYPE_ALIASES(FImpl, suffix)\ GAUGE_TYPE_ALIASES(FImpl, suffix)
SOLVER_TYPE_ALIASES(FImpl, suffix)
// logger // logger
class HadronsLogger: public Logger class HadronsLogger: public Logger

3
extras/Hadrons/Makefile.am
View File

@ -16,16 +16,17 @@ nobase_libHadrons_a_HEADERS = \
$(modules_hpp) \ $(modules_hpp) \
AllToAllVectors.hpp \ AllToAllVectors.hpp \
Application.hpp \ Application.hpp \
EigenPack.hpp \
Environment.hpp \ Environment.hpp \
Exceptions.hpp \ Exceptions.hpp \
Factory.hpp \ Factory.hpp \
GeneticScheduler.hpp \ GeneticScheduler.hpp \
Global.hpp \ Global.hpp \
Graph.hpp \ Graph.hpp \
EigenPack.hpp \
Module.hpp \ Module.hpp \
Modules.hpp \ Modules.hpp \
ModuleFactory.hpp \ ModuleFactory.hpp \
Solver.hpp \
VirtualMachine.hpp VirtualMachine.hpp
HadronsXmlRun_SOURCES = HadronsXmlRun.cc HadronsXmlRun_SOURCES = HadronsXmlRun.cc

92
extras/Hadrons/Modules.hpp
View File

@ -1,54 +1,56 @@
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp> #include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp> #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
#include <Grid/Hadrons/Modules/MContraction/MesonFieldGmu.hpp> #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/MesonFieldGamma.hpp>
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp> #include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp> #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
#include <Grid/Hadrons/Modules/MContraction/MesonFieldGmu.hpp>
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp> #include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp> #include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Grid/Hadrons/Modules/MSource/Point.hpp>
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
#include <Grid/Hadrons/Modules/MSink/Point.hpp>
#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp> #include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp> #include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp> #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp> #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
#include <Grid/Hadrons/Modules/MGauge/Random.hpp> #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.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/MAction/ZMobiusDWF.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>

2
extras/Hadrons/Modules/MAction/DWF.hpp
View File

@ -56,7 +56,7 @@ template <typename FImpl>
class TDWF: public Module<DWFPar> class TDWF: public Module<DWFPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TDWF(const std::string name); TDWF(const std::string name);

2
extras/Hadrons/Modules/MAction/Wilson.hpp
View File

@ -54,7 +54,7 @@ template <typename FImpl>
class TWilson: public Module<WilsonPar> class TWilson: public Module<WilsonPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TWilson(const std::string name); TWilson(const std::string name);

2
extras/Hadrons/Modules/MAction/WilsonClover.hpp
View File

@ -59,7 +59,7 @@ template <typename FImpl>
class TWilsonClover: public Module<WilsonCloverPar> class TWilsonClover: public Module<WilsonCloverPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TWilsonClover(const std::string name); TWilsonClover(const std::string name);

2
extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp
View File

@ -57,7 +57,7 @@ template <typename FImpl>
class TZMobiusDWF: public Module<ZMobiusDWFPar> class TZMobiusDWF: public Module<ZMobiusDWFPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TZMobiusDWF(const std::string name); TZMobiusDWF(const std::string name);

8
extras/Hadrons/Modules/MFermion/FreeProp.hpp
View File

@ -34,8 +34,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Module.hpp> #include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp> #include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
/****************************************************************************** /******************************************************************************
@ -57,7 +55,7 @@ template <typename FImpl>
class TFreeProp: public Module<FreePropPar> class TFreeProp: public Module<FreePropPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TFreeProp(const std::string name); TFreeProp(const std::string name);
@ -152,7 +150,7 @@ void TFreeProp<FImpl>::execute(void)
else else
{ {
PropToFerm<FImpl>(tmp, fullSrc, s, c); PropToFerm<FImpl>(tmp, fullSrc, s, c);
make_5D(tmp, source, Ls_); mat.ImportPhysicalFermionSource(tmp, source);
} }
} }
// source conversion for 5D sources // source conversion for 5D sources
@ -176,7 +174,7 @@ void TFreeProp<FImpl>::execute(void)
if (Ls_ > 1) if (Ls_ > 1)
{ {
PropagatorField &p4d = envGet(PropagatorField, getName()); PropagatorField &p4d = envGet(PropagatorField, getName());
make_4D(sol, tmp, Ls_); mat.ExportPhysicalFermionSolution(sol, tmp);
FermToProp<FImpl>(p4d, tmp, s, c); FermToProp<FImpl>(p4d, tmp, s, c);
} }
} }

2
extras/Hadrons/Modules/MFermion/GaugeProp.cc
View File

@ -32,4 +32,4 @@ using namespace Hadrons;
using namespace MFermion; using namespace MFermion;
template class Grid::Hadrons::MFermion::TGaugeProp<FIMPL>; template class Grid::Hadrons::MFermion::TGaugeProp<FIMPL>;
template class Grid::Hadrons::MFermion::TGaugeProp<ZFIMPL>;

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

@ -35,30 +35,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp> #include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp> #include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp> #include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Solver.hpp>
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* 5D -> 4D and 4D -> 5D conversions. *
******************************************************************************/
template<class vobj> // Note that 5D object is modified.
inline void make_4D(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
{
axpby_ssp_pminus(in_5d, 0., in_5d, 1., in_5d, 0, 0);
axpby_ssp_pplus(in_5d, 1., in_5d, 1., in_5d, 0, Ls-1);
ExtractSlice(out_4d, in_5d, 0, 0);
}
template<class vobj>
inline void make_5D(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
{
out_5d = zero;
InsertSlice(in_4d, out_5d, 0, 0);
InsertSlice(in_4d, out_5d, Ls-1, 0);
axpby_ssp_pplus(out_5d, 0., out_5d, 1., out_5d, 0, 0);
axpby_ssp_pminus(out_5d, 0., out_5d, 1., out_5d, Ls-1, Ls-1);
}
/****************************************************************************** /******************************************************************************
* GaugeProp * * GaugeProp *
******************************************************************************/ ******************************************************************************/
@ -76,7 +56,8 @@ template <typename FImpl>
class TGaugeProp: public Module<GaugePropPar> class TGaugeProp: public Module<GaugePropPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TGaugeProp(const std::string name); TGaugeProp(const std::string name);
@ -92,10 +73,12 @@ protected:
virtual void execute(void); virtual void execute(void);
private: private:
unsigned int Ls_; unsigned int Ls_;
SolverFn *solver_{nullptr}; Solver *solver_{nullptr};
}; };
MODULE_REGISTER_TMP(GaugeProp, TGaugeProp<FIMPL>, MFermion); MODULE_REGISTER_TMP(GaugeProp, TGaugeProp<FIMPL>, MFermion);
MODULE_REGISTER_TMP(ZGaugeProp, TGaugeProp<ZFIMPL>, MFermion);
/****************************************************************************** /******************************************************************************
* TGaugeProp implementation * * TGaugeProp implementation *
******************************************************************************/ ******************************************************************************/
@ -147,7 +130,8 @@ void TGaugeProp<FImpl>::execute(void)
std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d"); std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
auto &prop = envGet(PropagatorField, propName); auto &prop = envGet(PropagatorField, propName);
auto &fullSrc = envGet(PropagatorField, par().source); auto &fullSrc = envGet(PropagatorField, par().source);
auto &solver = envGet(SolverFn, par().solver); auto &solver = envGet(Solver, par().solver);
auto &mat = solver.getFMat();
envGetTmp(FermionField, source); envGetTmp(FermionField, source);
envGetTmp(FermionField, sol); envGetTmp(FermionField, sol);
@ -155,11 +139,12 @@ void TGaugeProp<FImpl>::execute(void)
LOG(Message) << "Inverting using solver '" << par().solver LOG(Message) << "Inverting using solver '" << par().solver
<< "' on source '" << par().source << "'" << std::endl; << "' on source '" << par().source << "'" << std::endl;
for (unsigned int s = 0; s < Ns; ++s) for (unsigned int s = 0; s < Ns; ++s)
for (unsigned int c = 0; c < FImpl::Dimension; ++c) for (unsigned int c = 0; c < FImpl::Dimension; ++c)
{ {
LOG(Message) << "Inversion for spin= " << s << ", color= " << c LOG(Message) << "Inversion for spin= " << s << ", color= " << c
<< std::endl; << std::endl;
// source conversion for 4D sources // source conversion for 4D sources
LOG(Message) << "Import source" << std::endl;
if (!env().isObject5d(par().source)) if (!env().isObject5d(par().source))
{ {
if (Ls_ == 1) if (Ls_ == 1)
@ -169,7 +154,7 @@ void TGaugeProp<FImpl>::execute(void)
else else
{ {
PropToFerm<FImpl>(tmp, fullSrc, s, c); PropToFerm<FImpl>(tmp, fullSrc, s, c);
make_5D(tmp, source, Ls_); mat.ImportPhysicalFermionSource(tmp, source);
} }
} }
// source conversion for 5D sources // source conversion for 5D sources
@ -184,14 +169,16 @@ void TGaugeProp<FImpl>::execute(void)
PropToFerm<FImpl>(source, fullSrc, s, c); PropToFerm<FImpl>(source, fullSrc, s, c);
} }
} }
LOG(Message) << "Solve" << std::endl;
sol = zero; sol = zero;
solver(sol, source); solver(sol, source);
LOG(Message) << "Export solution" << std::endl;
FermToProp<FImpl>(prop, sol, s, c); FermToProp<FImpl>(prop, sol, s, c);
// create 4D propagators from 5D one if necessary // create 4D propagators from 5D one if necessary
if (Ls_ > 1) if (Ls_ > 1)
{ {
PropagatorField &p4d = envGet(PropagatorField, getName()); PropagatorField &p4d = envGet(PropagatorField, getName());
make_4D(sol, tmp, Ls_); mat.ExportPhysicalFermionSolution(sol, tmp);
FermToProp<FImpl>(p4d, tmp, s, c); FermToProp<FImpl>(p4d, tmp, s, c);
} }
} }

16
extras/Hadrons/Modules/MGauge/StochEm.cc
View File

@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
Copyright (C) 2015-2018 Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.ac.uk>
Author: Vera Guelpers <vmg1n14@soton.ac.uk> Author: Vera Guelpers <vmg1n14@soton.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
@ -58,12 +59,8 @@ std::vector<std::string> TStochEm::getOutput(void)
// setup /////////////////////////////////////////////////////////////////////// // setup ///////////////////////////////////////////////////////////////////////
void TStochEm::setup(void) void TStochEm::setup(void)
{ {
create_weight = false; weightDone_ = env().hasCreatedObject("_" + getName() + "_weight");
if (!env().hasCreatedObject("_" + getName() + "_weight")) envCacheLat(EmComp, "_" + getName() + "_weight");
{
envCacheLat(EmComp, "_" + getName() + "_weight");
create_weight = true;
}
envCreateLat(EmField, getName()); envCreateLat(EmField, getName());
} }
@ -72,13 +69,14 @@ void TStochEm::execute(void)
{ {
LOG(Message) << "Generating stochastic EM potential..." << std::endl; LOG(Message) << "Generating stochastic EM potential..." << std::endl;
PhotonR photon(par().gauge, par().zmScheme); std::vector<Real> improvements = strToVec<Real>(par().improvement);
PhotonR photon(par().gauge, par().zmScheme, improvements, par().G0_qedInf);
auto &a = envGet(EmField, getName()); auto &a = envGet(EmField, getName());
auto &w = envGet(EmComp, "_" + getName() + "_weight"); auto &w = envGet(EmComp, "_" + getName() + "_weight");
if (create_weight) if (!weightDone_)
{ {
LOG(Message) << "Caching stochatic EM potential weight (gauge: " LOG(Message) << "Caching stochastic EM potential weight (gauge: "
<< par().gauge << ", zero-mode scheme: " << par().gauge << ", zero-mode scheme: "
<< par().zmScheme << ")..." << std::endl; << par().zmScheme << ")..." << std::endl;
photon.StochasticWeight(w); photon.StochasticWeight(w);

9
extras/Hadrons/Modules/MGauge/StochEm.hpp
View File

@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
Copyright (C) 2015-2018 Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.ac.uk>
Author: Vera Guelpers <vmg1n14@soton.ac.uk> Author: Vera Guelpers <vmg1n14@soton.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
@ -45,7 +46,9 @@ class StochEmPar: Serializable
public: public:
GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar, GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar,
PhotonR::Gauge, gauge, PhotonR::Gauge, gauge,
PhotonR::ZmScheme, zmScheme); PhotonR::ZmScheme, zmScheme,
std::string, improvement,
Real, G0_qedInf);
}; };
class TStochEm: public Module<StochEmPar> class TStochEm: public Module<StochEmPar>
@ -61,13 +64,13 @@ public:
// dependency relation // dependency relation
virtual std::vector<std::string> getInput(void); virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void); virtual std::vector<std::string> getOutput(void);
private:
bool create_weight;
protected: protected:
// setup // setup
virtual void setup(void); virtual void setup(void);
// execution // execution
virtual void execute(void); virtual void execute(void);
private:
bool weightDone_;
}; };
MODULE_REGISTER(StochEm, TStochEm, MGauge); MODULE_REGISTER(StochEm, TStochEm, MGauge);

69
extras/Hadrons/Modules/MGauge/UnitEm.cc Normal file
View File

@ -0,0 +1,69 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
Copyright (C) 2015
Copyright (C) 2016
Author: James Harrison <j.harrison@soton.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MGauge;
/******************************************************************************
* TStochEm implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TUnitEm::TUnitEm(const std::string name)
: Module<NoPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TUnitEm::getInput(void)
{
return std::vector<std::string>();
}
std::vector<std::string> TUnitEm::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TUnitEm::setup(void)
{
envCreateLat(EmField, getName());
}
// execution ///////////////////////////////////////////////////////////////////
void TUnitEm::execute(void)
{
PhotonR photon(0, 0); // Just chose arbitrary input values here
auto &a = envGet(EmField, getName());
LOG(Message) << "Generating unit EM potential..." << std::endl;
photon.UnitField(a);
}

69
extras/Hadrons/Modules/MGauge/UnitEm.hpp Normal file
View File

@ -0,0 +1,69 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
Copyright (C) 2015
Copyright (C) 2016
Author: James Harrison <j.harrison@soton.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MGauge_UnitEm_hpp_
#define Hadrons_MGauge_UnitEm_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* StochEm *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MGauge)
class TUnitEm: public Module<NoPar>
{
public:
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
public:
// constructor
TUnitEm(const std::string name);
// destructor
virtual ~TUnitEm(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER(UnitEm, TUnitEm, MGauge);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MGauge_UnitEm_hpp_

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

@ -51,7 +51,8 @@ std::vector<std::string> TChargedProp::getInput(void)
std::vector<std::string> TChargedProp::getOutput(void) std::vector<std::string> TChargedProp::getOutput(void)
{ {
std::vector<std::string> out = {getName()}; std::vector<std::string> out = {getName(), getName()+"_0", getName()+"_Q",
getName()+"_Sun", getName()+"_Tad"};
return out; return out;
} }
@ -66,18 +67,27 @@ void TChargedProp::setup(void)
phaseName_.push_back("_shiftphase_" + std::to_string(mu)); phaseName_.push_back("_shiftphase_" + std::to_string(mu));
} }
GFSrcName_ = getName() + "_DinvSrc"; GFSrcName_ = getName() + "_DinvSrc";
prop0Name_ = getName() + "_0";
propQName_ = getName() + "_Q";
propSunName_ = getName() + "_Sun";
propTadName_ = getName() + "_Tad";
fftName_ = getName() + "_fft"; fftName_ = getName() + "_fft";
freeMomPropDone_ = env().hasCreatedObject(freeMomPropName_); freeMomPropDone_ = env().hasCreatedObject(freeMomPropName_);
GFSrcDone_ = env().hasCreatedObject(GFSrcName_); GFSrcDone_ = env().hasCreatedObject(GFSrcName_);
phasesDone_ = env().hasCreatedObject(phaseName_[0]); phasesDone_ = env().hasCreatedObject(phaseName_[0]);
prop0Done_ = env().hasCreatedObject(prop0Name_);
envCacheLat(ScalarField, freeMomPropName_); envCacheLat(ScalarField, freeMomPropName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{ {
envCacheLat(ScalarField, phaseName_[mu]); envCacheLat(ScalarField, phaseName_[mu]);
} }
envCacheLat(ScalarField, GFSrcName_); envCacheLat(ScalarField, GFSrcName_);
envCacheLat(ScalarField, prop0Name_);
envCreateLat(ScalarField, getName()); envCreateLat(ScalarField, getName());
envCreateLat(ScalarField, propQName_);
envCreateLat(ScalarField, propSunName_);
envCreateLat(ScalarField, propTadName_);
envTmpLat(ScalarField, "buf"); envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "result"); envTmpLat(ScalarField, "result");
envTmpLat(ScalarField, "Amu"); envTmpLat(ScalarField, "Amu");
@ -95,79 +105,125 @@ void TChargedProp::execute(void)
<< " (mass= " << par().mass << " (mass= " << par().mass
<< ", charge= " << par().charge << ")..." << std::endl; << ", charge= " << par().charge << ")..." << std::endl;
auto &prop = envGet(ScalarField, getName()); auto &prop = envGet(ScalarField, getName());
auto &GFSrc = envGet(ScalarField, GFSrcName_); auto &prop0 = envGet(ScalarField, prop0Name_);
auto &G = envGet(ScalarField, freeMomPropName_); auto &propQ = envGet(ScalarField, propQName_);
auto &fft = envGet(FFT, fftName_); auto &propSun = envGet(ScalarField, propSunName_);
double q = par().charge; auto &propTad = envGet(ScalarField, propTadName_);
envGetTmp(ScalarField, result); auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &G = envGet(ScalarField, freeMomPropName_);
auto &fft = envGet(FFT, fftName_);
double q = par().charge;
envGetTmp(ScalarField, buf); envGetTmp(ScalarField, buf);
// G*F*Src // -G*momD1*G*F*Src (momD1 = F*D1*Finv)
prop = GFSrc; propQ = GFSrc;
momD1(propQ, fft);
propQ = -G*propQ;
propSun = -propQ;
fft.FFT_dim(propQ, propQ, env().getNd()-1, FFT::backward);
// - q*G*momD1*G*F*Src (momD1 = F*D1*Finv) // G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
buf = GFSrc; momD1(propSun, fft);
momD1(buf, fft); propSun = G*propSun;
buf = G*buf; fft.FFT_dim(propSun, propSun, env().getNd()-1, FFT::backward);
prop = prop - q*buf;
// + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) // -G*momD2*G*F*Src (momD2 = F*D2*Finv)
momD1(buf, fft); propTad = GFSrc;
prop = prop + q*q*G*buf; momD2(propTad, fft);
propTad = -G*propTad;
fft.FFT_dim(propTad, propTad, env().getNd()-1, FFT::backward);
// - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv) // full charged scalar propagator
buf = GFSrc; fft.FFT_dim(buf, GFSrc, env().getNd()-1, FFT::backward);
momD2(buf, fft); prop = buf + q*propQ + q*q*propSun + q*q*propTad;
prop = prop - q*q*G*buf;
// final FT
fft.FFT_all_dim(prop, prop, FFT::backward);
// OUTPUT IF NECESSARY // OUTPUT IF NECESSARY
if (!par().output.empty()) if (!par().output.empty())
{ {
std::string filename = par().output + "." + Result result;
std::to_string(vm().getTrajectory()); TComplex site;
std::vector<int> siteCoor;
LOG(Message) << "Saving zero-momentum projection to '" LOG(Message) << "Saving momentum-projected propagator to '"
<< filename << "'..." << std::endl; << RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
std::vector<TComplex> vecBuf; result.projection.resize(par().outputMom.size());
std::vector<Complex> result; result.lattice_size = env().getGrid()->_fdimensions;
result.mass = par().mass;
sliceSum(prop, vecBuf, Tp); result.charge = q;
result.resize(vecBuf.size()); siteCoor.resize(env().getNd());
for (unsigned int t = 0; t < vecBuf.size(); ++t) for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{ {
result[t] = TensorRemove(vecBuf[t]); result.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
LOG(Message) << "Calculating (" << par().outputMom[i_p]
<< ") momentum projection" << std::endl;
result.projection[i_p].corr_0.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Q.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Sun.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
result.projection[i_p].corr_Tad.resize(env().getGrid()->_fdimensions[env().getNd()-1]);
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
siteCoor[j] = result.projection[i_p].momentum[j];
}
for (unsigned int t = 0; t < result.projection[i_p].corr.size(); ++t)
{
siteCoor[env().getNd()-1] = t;
peekSite(site, prop, siteCoor);
result.projection[i_p].corr[t]=TensorRemove(site);
peekSite(site, buf, siteCoor);
result.projection[i_p].corr_0[t]=TensorRemove(site);
peekSite(site, propQ, siteCoor);
result.projection[i_p].corr_Q[t]=TensorRemove(site);
peekSite(site, propSun, siteCoor);
result.projection[i_p].corr_Sun[t]=TensorRemove(site);
peekSite(site, propTad, siteCoor);
result.projection[i_p].corr_Tad[t]=TensorRemove(site);
}
} }
saveResult(par().output, "charge", q);
saveResult(par().output, "prop", result); saveResult(par().output, "prop", result);
} }
std::vector<int> mask(env().getNd(),1);
mask[env().getNd()-1] = 0;
fft.FFT_dim_mask(prop, prop, mask, FFT::backward);
fft.FFT_dim_mask(propQ, propQ, mask, FFT::backward);
fft.FFT_dim_mask(propSun, propSun, mask, FFT::backward);
fft.FFT_dim_mask(propTad, propTad, mask, FFT::backward);
} }
void TChargedProp::makeCaches(void) void TChargedProp::makeCaches(void)
{ {
auto &freeMomProp = envGet(ScalarField, freeMomPropName_); auto &freeMomProp = envGet(ScalarField, freeMomPropName_);
auto &GFSrc = envGet(ScalarField, GFSrcName_); auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &prop0 = envGet(ScalarField, prop0Name_);
auto &fft = envGet(FFT, fftName_); auto &fft = envGet(FFT, fftName_);
if (!freeMomPropDone_) if (!freeMomPropDone_)
{ {
LOG(Message) << "Caching momentum space free scalar propagator" LOG(Message) << "Caching momentum-space free scalar propagator"
<< " (mass= " << par().mass << ")..." << std::endl; << " (mass= " << par().mass << ")..." << std::endl;
SIMPL::MomentumSpacePropagator(freeMomProp, par().mass); SIMPL::MomentumSpacePropagator(freeMomProp, par().mass);
} }
if (!GFSrcDone_) if (!GFSrcDone_)
{ {
FFT fft(env().getGrid());
auto &source = envGet(ScalarField, par().source); auto &source = envGet(ScalarField, par().source);
LOG(Message) << "Caching G*F*src..." << std::endl; LOG(Message) << "Caching G*F*src..." << std::endl;
fft.FFT_all_dim(GFSrc, source, FFT::forward); fft.FFT_all_dim(GFSrc, source, FFT::forward);
GFSrc = freeMomProp*GFSrc; GFSrc = freeMomProp*GFSrc;
} }
if (!prop0Done_)
{
LOG(Message) << "Caching position-space free scalar propagator..."
<< std::endl;
fft.FFT_all_dim(prop0, GFSrc, FFT::backward);
}
if (!phasesDone_) if (!phasesDone_)
{ {
std::vector<int> &l = env().getGrid()->_fdimensions; std::vector<int> &l = env().getGrid()->_fdimensions;
@ -184,6 +240,14 @@ void TChargedProp::makeCaches(void)
phase_.push_back(&phmu); phase_.push_back(&phmu);
} }
} }
else
{
phase_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
}
}
} }
void TChargedProp::momD1(ScalarField &s, FFT &fft) void TChargedProp::momD1(ScalarField &s, FFT &fft)

30
extras/Hadrons/Modules/MScalar/ChargedProp.hpp
View File

@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MScalar/ChargedProp.hpp
Copyright (C) 2015-2018 Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <j.harrison@soton.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -47,7 +48,8 @@ public:
std::string, source, std::string, source,
double, mass, double, mass,
double, charge, double, charge,
std::string, output); std::string, output,
std::vector<std::string>, outputMom);
}; };
class TChargedProp: public Module<ChargedPropPar> class TChargedProp: public Module<ChargedPropPar>
@ -56,6 +58,26 @@ public:
SCALAR_TYPE_ALIASES(SIMPL,); SCALAR_TYPE_ALIASES(SIMPL,);
typedef PhotonR::GaugeField EmField; typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp; typedef PhotonR::GaugeLinkField EmComp;
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<Complex>, corr,
std::vector<Complex>, corr_0,
std::vector<Complex>, corr_Q,
std::vector<Complex>, corr_Sun,
std::vector<Complex>, corr_Tad);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
double, charge,
std::vector<Projection>, projection);
};
public: public:
// constructor // constructor
TChargedProp(const std::string name); TChargedProp(const std::string name);
@ -74,8 +96,10 @@ private:
void momD1(ScalarField &s, FFT &fft); void momD1(ScalarField &s, FFT &fft);
void momD2(ScalarField &s, FFT &fft); void momD2(ScalarField &s, FFT &fft);
private: private:
bool freeMomPropDone_, GFSrcDone_, phasesDone_; bool freeMomPropDone_, GFSrcDone_, prop0Done_,
std::string freeMomPropName_, GFSrcName_, fftName_; phasesDone_;
std::string freeMomPropName_, GFSrcName_, prop0Name_,
propQName_, propSunName_, propTadName_, fftName_;
std::vector<std::string> phaseName_; std::vector<std::string> phaseName_;
std::vector<ScalarField *> phase_; std::vector<ScalarField *> phase_;
}; };

536
extras/Hadrons/Modules/MScalar/ScalarVP.cc Normal file
View File

@ -0,0 +1,536 @@
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalar;
/*
* Scalar QED vacuum polarisation up to O(alpha)
*
* Conserved vector 2-point function diagram notation:
* _______
* / \
* U_nu * * U_mu
* \_______/
*
* ( adj(S(a\hat{nu}|x)) U_mu(x) S(0|x+a\hat{mu}) U_nu(0) )
* = 2 Re( - )
* ( adj(S(a\hat{nu}|x+a\hat{mu})) adj(U_mu(x)) S(0|x) U_nu(0) )
*
*
* _______
* / \
* free = 1 * * 1
* \_______/
*
*
*
* _______
* / \
* S = iA_nu * * iA_mu
* \_______/
*
*
* Delta_1
* ___*___
* / \
* X = 1 * * 1
* \___*___/
* Delta_1
*
* Delta_1 Delta_1
* ___*___ ___*___
* / \ / \
* 1 * * iA_mu + iA_nu * * 1
* \_______/ \_______/
* 4C = _______ _______
* / \ / \
* + 1 * * iA_mu + iA_nu * * 1
* \___*___/ \___*___/
* Delta_1 Delta_1
*
* Delta_1 Delta_1
* _*___*_ _______
* / \ / \
* 2E = 1 * * 1 + 1 * * 1
* \_______/ \_*___*_/
* Delta_1 Delta_1
*
* Delta_2
* ___*___ _______
* / \ / \
* 2T = 1 * * 1 + 1 * * 1
* \_______/ \___*___/
* Delta_2
*
*
* _______
* / \
* srcT = -A_nu^2/2 * * 1
* \_______/
*
*
*
* _______
* / \
* snkT = 1 * * -A_mu^2/2
* \_______/
*
* Full VP to O(alpha) = free + q^2*(S+X+4C+2E+2T+srcT+snkT)
*/
/******************************************************************************
* TScalarVP implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TScalarVP::TScalarVP(const std::string name)
: Module<ScalarVPPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TScalarVP::getInput(void)
{
prop0Name_ = par().scalarProp + "_0";
propQName_ = par().scalarProp + "_Q";
propSunName_ = par().scalarProp + "_Sun";
propTadName_ = par().scalarProp + "_Tad";
std::vector<std::string> in = {par().emField, prop0Name_, propQName_,
propSunName_, propTadName_};
return in;
}
std::vector<std::string> TScalarVP::getOutput(void)
{
std::vector<std::string> out;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
// out.push_back(getName() + "_propQ_" + std::to_string(mu));
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
out.push_back(getName() + "_" + std::to_string(mu)
+ "_" + std::to_string(nu));
}
}
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TScalarVP::setup(void)
{
freeMomPropName_ = FREEMOMPROP(static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass);
GFSrcName_ = par().scalarProp + "_DinvSrc";
fftName_ = par().scalarProp + "_fft";
phaseName_.clear();
muPropQName_.clear();
vpTensorName_.clear();
momPhaseName_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
muPropQName_.push_back(getName() + "_propQ_" + std::to_string(mu));
std::vector<std::string> vpTensorName_mu;
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
vpTensorName_mu.push_back(getName() + "_" + std::to_string(mu)
+ "_" + std::to_string(nu));
}
vpTensorName_.push_back(vpTensorName_mu);
}
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
}
}
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
envCreateLat(ScalarField, muPropQName_[mu]);
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
envCreateLat(ScalarField, vpTensorName_[mu][nu]);
}
}
if (!par().output.empty())
{
momPhasesDone_ = env().hasCreatedObject(momPhaseName_[0]);
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
envCacheLat(ScalarField, momPhaseName_[i_p]);
}
}
envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "result");
envTmpLat(ScalarField, "Amu");
envTmpLat(ScalarField, "Usnk");
envTmpLat(ScalarField, "tmpProp");
}
// execution ///////////////////////////////////////////////////////////////////
void TScalarVP::execute(void)
{
// CACHING ANALYTIC EXPRESSIONS
makeCaches();
Complex ci(0.0,1.0);
Real q = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().charge;
auto &prop0 = envGet(ScalarField, prop0Name_);
auto &propQ = envGet(ScalarField, propQName_);
auto &propSun = envGet(ScalarField, propSunName_);
auto &propTad = envGet(ScalarField, propTadName_);
auto &GFSrc = envGet(ScalarField, GFSrcName_);
auto &G = envGet(ScalarField, freeMomPropName_);
auto &fft = envGet(FFT, fftName_);
phase_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
auto &phmu = envGet(ScalarField, phaseName_[mu]);
phase_.push_back(&phmu);
}
// PROPAGATORS FROM SHIFTED SOURCES
LOG(Message) << "Computing O(q) charged scalar propagators..."
<< std::endl;
std::vector<ScalarField *> muPropQ;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
auto &propmu = envGet(ScalarField, muPropQName_[mu]);
// -G*momD1*G*F*tau_mu*Src (momD1 = F*D1*Finv)
propmu = adj(*phase_[mu])*GFSrc;
momD1(propmu, fft);
propmu = -G*propmu;
fft.FFT_all_dim(propmu, propmu, FFT::backward);
muPropQ.push_back(&propmu);
}
// CONTRACTIONS
auto &A = envGet(EmField, par().emField);
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, result);
envGetTmp(ScalarField, Amu);
envGetTmp(ScalarField, Usnk);
envGetTmp(ScalarField, tmpProp);
TComplex Anu0, Usrc;
std::vector<int> coor0 = {0, 0, 0, 0};
std::vector<std::vector<ScalarField *> > vpTensor;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
std::vector<ScalarField *> vpTensor_mu;
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
auto &vpmunu = envGet(ScalarField, vpTensorName_[mu][nu]);
vpTensor_mu.push_back(&vpmunu);
}
vpTensor.push_back(vpTensor_mu);
}
// Prepare output data structure if necessary
Result outputData;
if (!par().output.empty())
{
outputData.projection.resize(par().outputMom.size());
outputData.lattice_size = env().getGrid()->_fdimensions;
outputData.mass = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass;
outputData.charge = q;
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
outputData.projection[i_p].pi.resize(env().getNd());
outputData.projection[i_p].pi_free.resize(env().getNd());
outputData.projection[i_p].pi_2E.resize(env().getNd());
outputData.projection[i_p].pi_2T.resize(env().getNd());
outputData.projection[i_p].pi_S.resize(env().getNd());
outputData.projection[i_p].pi_4C.resize(env().getNd());
outputData.projection[i_p].pi_X.resize(env().getNd());
outputData.projection[i_p].pi_srcT.resize(env().getNd());
outputData.projection[i_p].pi_snkT.resize(env().getNd());
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
outputData.projection[i_p].pi[nu].resize(env().getNd());
outputData.projection[i_p].pi_free[nu].resize(env().getNd());
outputData.projection[i_p].pi_2E[nu].resize(env().getNd());
outputData.projection[i_p].pi_2T[nu].resize(env().getNd());
outputData.projection[i_p].pi_S[nu].resize(env().getNd());
outputData.projection[i_p].pi_4C[nu].resize(env().getNd());
outputData.projection[i_p].pi_X[nu].resize(env().getNd());
outputData.projection[i_p].pi_srcT[nu].resize(env().getNd());
outputData.projection[i_p].pi_snkT[nu].resize(env().getNd());
}
}
}
// Do contractions
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
peekSite(Anu0, peekLorentz(A, nu), coor0);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
LOG(Message) << "Computing Pi[" << mu << "][" << nu << "]..."
<< std::endl;
Amu = peekLorentz(A, mu);
// free
tmpProp = Cshift(prop0, nu, -1); // S_0(0|x-a\hat{\nu})
// = S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, prop0, tmpProp, Usrc, mu);
*vpTensor[mu][nu] = result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_free[mu][nu], result,
i_p);
}
}
tmpProp = result; // Just using tmpProp as a temporary ScalarField
// here (buf is modified by calls to writeVP())
// srcT
result = tmpProp * (-0.5)*Anu0*Anu0;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_srcT[mu][nu], result,
i_p);
}
}
// snkT
result = tmpProp * (-0.5)*Amu*Amu;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_snkT[mu][nu], result,
i_p);
}
}
// S
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = ci*Anu0;
Usnk = ci*Amu;
vpContraction(result, prop0, tmpProp, Usrc, Usnk, mu);
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_S[mu][nu], result,
i_p);
}
}
// 4C
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
Usnk = ci*Amu;
vpContraction(result, propQ, tmpProp, Usrc, Usnk, mu);
Usrc = ci*Anu0;
vpContraction(buf, propQ, tmpProp, Usrc, mu);
result += buf;
vpContraction(buf, prop0, *muPropQ[nu], Usrc, mu);
result += buf;
Usrc = Complex(1.0,0.0);
Usnk = ci*Amu;
vpContraction(buf, prop0, *muPropQ[nu], Usrc, Usnk, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_4C[mu][nu], result,
i_p);
}
}
// X
Usrc = Complex(1.0,0.0);
vpContraction(result, propQ, *muPropQ[nu], Usrc, mu);
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_X[mu][nu], result,
i_p);
}
}
// 2E
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, propSun, tmpProp, Usrc, mu);
tmpProp = Cshift(propSun, nu, -1); // S_\Sigma(0|x-a\hat{\nu})
//(Note: <S(0|x-a\hat{\nu})> = <S(a\hat{\nu}|x)>)
vpContraction(buf, prop0, tmpProp, Usrc, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_2E[mu][nu], result,
i_p);
}
}
// 2T
tmpProp = Cshift(prop0, nu, -1); // S_0(a\hat{\nu}|x)
Usrc = Complex(1.0,0.0);
vpContraction(result, propTad, tmpProp, Usrc, mu);
tmpProp = Cshift(propTad, nu, -1); // S_T(0|x-a\hat{\nu})
vpContraction(buf, prop0, tmpProp, Usrc, mu);
result += buf;
*vpTensor[mu][nu] += q*q*result;
// Do momentum projections if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi_2T[mu][nu], result,
i_p);
}
}
// Do momentum projections of full VP if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pi[mu][nu],
*vpTensor[mu][nu], i_p);
}
}
}
}
// OUTPUT IF NECESSARY
if (!par().output.empty())
{
LOG(Message) << "Saving momentum-projected HVP to '"
<< RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
saveResult(par().output, "HVP", outputData);
}
}
void TScalarVP::makeCaches(void)
{
envGetTmp(ScalarField, buf);
if ( (!par().output.empty()) && (!momPhasesDone_) )
{
LOG(Message) << "Caching phases for momentum projections..."
<< std::endl;
std::vector<int> &l = env().getGrid()->_fdimensions;
Complex ci(0.0,1.0);
// Calculate phase factors
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
std::vector<int> mom = strToVec<int>(par().outputMom[i_p]);
auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
momph_ip = zero;
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
Real twoPiL = M_PI*2./l[j];
LatticeCoordinate(buf, j);
buf = mom[j]*twoPiL*buf;
momph_ip = momph_ip + buf;
}
momph_ip = exp(-ci*momph_ip);
momPhase_.push_back(&momph_ip);
}
}
}
void TScalarVP::vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, ScalarField &u_snk, int mu)
{
// Note: this function assumes a point source is used.
vp = adj(prop_nu_x) * u_snk * Cshift(prop_0_x, mu, 1) * u_src;
vp -= Cshift(adj(prop_nu_x), mu, 1) * adj(u_snk) * prop_0_x * u_src;
vp = 2.0*real(vp);
}
void TScalarVP::vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, int mu)
{
// Note: this function assumes a point source is used.
vp = adj(prop_nu_x) * Cshift(prop_0_x, mu, 1) * u_src;
vp -= Cshift(adj(prop_nu_x), mu, 1) * prop_0_x * u_src;
vp = 2.0*real(vp);
}
void TScalarVP::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
{
std::vector<TComplex> vecBuf;
envGetTmp(ScalarField, buf);
buf = vp*(*momPhase_[i_p]);
sliceSum(buf, vecBuf, Tp);
projection.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
projection[t] = TensorRemove(vecBuf[t]);
}
}
void TScalarVP::momD1(ScalarField &s, FFT &fft)
{
auto &A = envGet(EmField, par().emField);
Complex ci(0.0,1.0);
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, result);
envGetTmp(ScalarField, Amu);
result = zero;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Amu = peekLorentz(A, mu);
buf = (*phase_[mu])*s;
fft.FFT_all_dim(buf, buf, FFT::backward);
buf = Amu*buf;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result - ci*buf;
}
fft.FFT_all_dim(s, s, FFT::backward);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Amu = peekLorentz(A, mu);
buf = Amu*s;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result + ci*adj(*phase_[mu])*buf;
}
s = result;
}

101
extras/Hadrons/Modules/MScalar/ScalarVP.hpp Normal file
View File

@ -0,0 +1,101 @@
#ifndef Hadrons_MScalar_ScalarVP_hpp_
#define Hadrons_MScalar_ScalarVP_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* Scalar vacuum polarisation *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalar)
class ScalarVPPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar,
std::string, emField,
std::string, scalarProp,
std::string, output,
std::vector<std::string>, outputMom);
};
class TScalarVP: public Module<ScalarVPPar>
{
public:
SCALAR_TYPE_ALIASES(SIMPL,);
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<std::vector<std::vector<Complex>>>, pi,
std::vector<std::vector<std::vector<Complex>>>, pi_free,
std::vector<std::vector<std::vector<Complex>>>, pi_2E,
std::vector<std::vector<std::vector<Complex>>>, pi_2T,
std::vector<std::vector<std::vector<Complex>>>, pi_S,
std::vector<std::vector<std::vector<Complex>>>, pi_4C,
std::vector<std::vector<std::vector<Complex>>>, pi_X,
std::vector<std::vector<std::vector<Complex>>>, pi_srcT,
std::vector<std::vector<std::vector<Complex>>>, pi_snkT);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
double, charge,
std::vector<Projection>, projection);
};
public:
// constructor
TScalarVP(const std::string name);
// destructor
virtual ~TScalarVP(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
void makeCaches(void);
// conserved vector two-point contraction
void vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, ScalarField &u_snk, int mu);
// conserved vector two-point contraction with unit gauge link at sink
void vpContraction(ScalarField &vp,
ScalarField &prop_0_x, ScalarField &prop_nu_x,
TComplex u_src, int mu);
// write momentum-projected vacuum polarisation to file(s)
void project(std::vector<Complex> &projection, const ScalarField &vp,
int i_p);
// momentum-space Delta_1 insertion
void momD1(ScalarField &s, FFT &fft);
private:
bool momPhasesDone_;
std::string freeMomPropName_, GFSrcName_,
prop0Name_, propQName_,
propSunName_, propTadName_,
fftName_;
std::vector<std::string> phaseName_, muPropQName_,
momPhaseName_;
std::vector<std::vector<std::string> > vpTensorName_;
std::vector<ScalarField *> phase_, momPhase_;
};
MODULE_REGISTER(ScalarVP, TScalarVP, MScalar);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalar_ScalarVP_hpp_

232
extras/Hadrons/Modules/MScalar/VPCounterTerms.cc Normal file
View File

@ -0,0 +1,232 @@
#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MScalar;
/******************************************************************************
* TVPCounterTerms implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TVPCounterTerms::TVPCounterTerms(const std::string name)
: Module<VPCounterTermsPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TVPCounterTerms::getInput(void)
{
std::vector<std::string> in = {par().source};
return in;
}
std::vector<std::string> TVPCounterTerms::getOutput(void)
{
std::vector<std::string> out;
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
void TVPCounterTerms::setup(void)
{
freeMomPropName_ = FREEMOMPROP(par().mass);
phaseName_.clear();
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
}
GFSrcName_ = getName() + "_DinvSrc";
phatsqName_ = getName() + "_pHatSquared";
prop0Name_ = getName() + "_freeProp";
twoscalarName_ = getName() + "_2scalarProp";
psquaredName_ = getName() + "_psquaredProp";
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
}
}
envCreateLat(ScalarField, freeMomPropName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
envCreateLat(ScalarField, phaseName_[mu]);
}
envCreateLat(ScalarField, phatsqName_);
envCreateLat(ScalarField, GFSrcName_);
envCreateLat(ScalarField, prop0Name_);
envCreateLat(ScalarField, twoscalarName_);
envCreateLat(ScalarField, psquaredName_);
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
envCacheLat(ScalarField, momPhaseName_[i_p]);
}
}
envTmpLat(ScalarField, "buf");
envTmpLat(ScalarField, "tmp_vp");
envTmpLat(ScalarField, "vpPhase");
}
// execution ///////////////////////////////////////////////////////////////////
void TVPCounterTerms::execute(void)
{
auto &source = envGet(ScalarField, par().source);
Complex ci(0.0,1.0);
FFT fft(env().getGrid());
envGetTmp(ScalarField, buf);
envGetTmp(ScalarField, tmp_vp);
// Momentum-space free scalar propagator
auto &G = envGet(ScalarField, freeMomPropName_);
SIMPL::MomentumSpacePropagator(G, par().mass);
// Phases and hat{p}^2
auto &phatsq = envGet(ScalarField, phatsqName_);
std::vector<int> &l = env().getGrid()->_fdimensions;
LOG(Message) << "Calculating shift phases..." << std::endl;
phatsq = zero;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Real twoPiL = M_PI*2./l[mu];
auto &phmu = envGet(ScalarField, phaseName_[mu]);
LatticeCoordinate(buf, mu);
phmu = exp(ci*twoPiL*buf);
phase_.push_back(&phmu);
buf = 2.*sin(.5*twoPiL*buf);
phatsq = phatsq + buf*buf;
}
// G*F*src
auto &GFSrc = envGet(ScalarField, GFSrcName_);
fft.FFT_all_dim(GFSrc, source, FFT::forward);
GFSrc = G*GFSrc;
// Position-space free scalar propagator
auto &prop0 = envGet(ScalarField, prop0Name_);
prop0 = GFSrc;
fft.FFT_all_dim(prop0, prop0, FFT::backward);
// Propagators for counter-terms
auto &twoscalarProp = envGet(ScalarField, twoscalarName_);
auto &psquaredProp = envGet(ScalarField, psquaredName_);
twoscalarProp = G*GFSrc;
fft.FFT_all_dim(twoscalarProp, twoscalarProp, FFT::backward);
psquaredProp = G*phatsq*GFSrc;
fft.FFT_all_dim(psquaredProp, psquaredProp, FFT::backward);
// Prepare output data structure if necessary
Result outputData;
if (!par().output.empty())
{
outputData.projection.resize(par().outputMom.size());
outputData.lattice_size = env().getGrid()->_fdimensions;
outputData.mass = par().mass;
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
outputData.projection[i_p].twoScalar.resize(env().getNd());
outputData.projection[i_p].threeScalar.resize(env().getNd());
outputData.projection[i_p].pSquaredInsertion.resize(env().getNd());
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
outputData.projection[i_p].twoScalar[nu].resize(env().getNd());
outputData.projection[i_p].threeScalar[nu].resize(env().getNd());
outputData.projection[i_p].pSquaredInsertion[nu].resize(env().getNd());
}
// Calculate phase factors
auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
momph_ip = zero;
for (unsigned int j = 0; j < env().getNd()-1; ++j)
{
Real twoPiL = M_PI*2./l[j];
LatticeCoordinate(buf, j);
buf = outputData.projection[i_p].momentum[j]*twoPiL*buf;
momph_ip = momph_ip + buf;
}
momph_ip = exp(-ci*momph_ip);
momPhase_.push_back(&momph_ip);
}
}
// Contractions
for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{
buf = adj(Cshift(prop0, nu, -1));
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
// Two-scalar loop
tmp_vp = buf * Cshift(prop0, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * prop0;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].twoScalar[mu][nu],
tmp_vp, i_p);
}
}
// Three-scalar loop (no vertex)
tmp_vp = buf * Cshift(twoscalarProp, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * twoscalarProp;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].threeScalar[mu][nu],
tmp_vp, i_p);
}
}
// Three-scalar loop (hat{p}^2 insertion)
tmp_vp = buf * Cshift(psquaredProp, mu, 1);
tmp_vp -= Cshift(buf, mu, 1) * psquaredProp;
tmp_vp = 2.0*real(tmp_vp);
// Output if necessary
if (!par().output.empty())
{
for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
{
project(outputData.projection[i_p].pSquaredInsertion[mu][nu],
tmp_vp, i_p);
}
}
}
}
// OUTPUT IF NECESSARY
if (!par().output.empty())
{
LOG(Message) << "Saving momentum-projected correlators to '"
<< RESULT_FILE_NAME(par().output) << "'..."
<< std::endl;
saveResult(par().output, "scalar_loops", outputData);
}
}
void TVPCounterTerms::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
{
std::vector<TComplex> vecBuf;
envGetTmp(ScalarField, vpPhase);
vpPhase = vp*(*momPhase_[i_p]);
sliceSum(vpPhase, vecBuf, Tp);
projection.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
projection[t] = TensorRemove(vecBuf[t]);
}
}

75
extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp Normal file
View File

@ -0,0 +1,75 @@
#ifndef Hadrons_MScalar_VPCounterTerms_hpp_
#define Hadrons_MScalar_VPCounterTerms_hpp_
#include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/******************************************************************************
* VPCounterTerms *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MScalar)
class VPCounterTermsPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(VPCounterTermsPar,
std::string, source,
double, mass,
std::string, output,
std::vector<std::string>, outputMom);
};
class TVPCounterTerms: public Module<VPCounterTermsPar>
{
public:
SCALAR_TYPE_ALIASES(SIMPL,);
class Result: Serializable
{
public:
class Projection: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
std::vector<int>, momentum,
std::vector<std::vector<std::vector<Complex>>>, twoScalar,
std::vector<std::vector<std::vector<Complex>>>, threeScalar,
std::vector<std::vector<std::vector<Complex>>>, pSquaredInsertion);
};
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<int>, lattice_size,
double, mass,
std::vector<Projection>, projection);
};
public:
// constructor
TVPCounterTerms(const std::string name);
// destructor
virtual ~TVPCounterTerms(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
void project(std::vector<Complex> &projection, const ScalarField &vp, int i_p);
private:
std::string freeMomPropName_, GFSrcName_, phatsqName_, prop0Name_,
twoscalarName_, twoscalarVertexName_,
psquaredName_, psquaredVertexName_;
std::vector<std::string> phaseName_, momPhaseName_;
std::vector<ScalarField *> phase_, momPhase_;
};
MODULE_REGISTER(VPCounterTerms, TVPCounterTerms, MScalar);
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MScalar_VPCounterTerms_hpp_

10
extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
View File

@ -32,6 +32,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/Global.hpp> #include <Grid/Hadrons/Global.hpp>
#include <Grid/Hadrons/Module.hpp> #include <Grid/Hadrons/Module.hpp>
#include <Grid/Hadrons/ModuleFactory.hpp> #include <Grid/Hadrons/ModuleFactory.hpp>
#include <Grid/Hadrons/Solver.hpp>
#include <Grid/Hadrons/EigenPack.hpp> #include <Grid/Hadrons/EigenPack.hpp>
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
@ -55,7 +56,8 @@ template <typename FImpl, int nBasis>
class TRBPrecCG: public Module<RBPrecCGPar> class TRBPrecCG: public Module<RBPrecCGPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
SOLVER_TYPE_ALIASES(FImpl,);
typedef FermionEigenPack<FImpl> EPack; typedef FermionEigenPack<FImpl> EPack;
typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack; typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack;
typedef std::shared_ptr<Guesser<FermionField>> GuesserPt; typedef std::shared_ptr<Guesser<FermionField>> GuesserPt;
@ -176,14 +178,12 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
schurSolver(mat, source, sol, *guesser); schurSolver(mat, source, sol, *guesser);
}; };
}; };
auto solver = makeSolver(false); auto solver = makeSolver(false);
envCreate(SolverFn, getName(), Ls, solver); envCreate(Solver, getName(), Ls, solver, mat);
auto solver_subtract = makeSolver(true); auto solver_subtract = makeSolver(true);
envCreate(SolverFn, getName() + "_subtract", Ls, solver_subtract); envCreate(Solver, getName() + "_subtract", Ls, solver_subtract, mat);
} }
// execution /////////////////////////////////////////////////////////////////// // execution ///////////////////////////////////////////////////////////////////
template <typename FImpl, int nBasis> template <typename FImpl, int nBasis>
void TRBPrecCG<FImpl, nBasis>::execute(void) void TRBPrecCG<FImpl, nBasis>::execute(void)

2
extras/Hadrons/Modules/MSource/SeqGamma.hpp
View File

@ -71,7 +71,7 @@ template <typename FImpl>
class TSeqGamma: public Module<SeqGammaPar> class TSeqGamma: public Module<SeqGammaPar>
{ {
public: public:
FGS_TYPE_ALIASES(FImpl,); FG_TYPE_ALIASES(FImpl,);
public: public:
// constructor // constructor
TSeqGamma(const std::string name); TSeqGamma(const std::string name);

62
extras/Hadrons/Solver.hpp Normal file
View File

@ -0,0 +1,62 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/EigenPack.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_Solver_hpp_
#define Hadrons_Solver_hpp_
#include <Grid/Hadrons/Global.hpp>
BEGIN_HADRONS_NAMESPACE
template <typename FImpl>
class Solver
{
public:
typedef typename FImpl::FermionField FermionField;
typedef FermionOperator<FImpl> FMat;
typedef std::function<void(FermionField &,
const FermionField &)> SolverFn;
public:
Solver(SolverFn fn, FMat &mat): mat_(mat), fn_(fn) {}
void operator()(FermionField &sol, const FermionField &src)
{
fn_(sol, src);
}
FMat & getFMat(void)
{
return mat_;
}
private:
FMat &mat_;
SolverFn fn_;
};
END_HADRONS_NAMESPACE
#endif // Hadrons_Solver_hpp_

174
extras/Hadrons/modules.inc
View File

@ -1,107 +1,113 @@
modules_cc =\ modules_cc =\
Modules/MSource/SeqConserved.cc \ Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MSource/SeqGamma.cc \ Modules/MContraction/Baryon.cc \
Modules/MSource/Point.cc \ Modules/MContraction/MesonFieldGamma.cc \
Modules/MSource/Z2.cc \ Modules/MContraction/Meson.cc \
Modules/MSource/Wall.cc \ Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MIO/LoadCoarseEigenPack.cc \
Modules/MIO/LoadBinary.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/Div.cc \
Modules/MScalarSUN/ShiftProbe.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/TimeMomProbe.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MAction/DWF.cc \
Modules/MAction/ZMobiusDWF.cc \
Modules/MAction/Wilson.cc \
Modules/MAction/WilsonClover.cc \
Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MContraction/WeakHamiltonianNonEye.cc \
Modules/MContraction/WardIdentity.cc \ Modules/MContraction/WardIdentity.cc \
Modules/MContraction/WeakHamiltonianEye.cc \
Modules/MContraction/DiscLoop.cc \ Modules/MContraction/DiscLoop.cc \
Modules/MContraction/Baryon.cc \
Modules/MContraction/MesonFieldGmu.cc \
Modules/MContraction/Gamma3pt.cc \ Modules/MContraction/Gamma3pt.cc \
Modules/MContraction/WeakNeutral4ptDisc.cc \
Modules/MContraction/Meson.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/FreeProp.cc \
Modules/MUtilities/TestSeqGamma.cc \
Modules/MUtilities/TestSeqConserved.cc \
Modules/MFermion/FreeProp.cc \ Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \ Modules/MFermion/GaugeProp.cc \
Modules/MSource/Point.cc \
Modules/MSource/Wall.cc \
Modules/MSource/SeqConserved.cc \
Modules/MSource/SeqGamma.cc \
Modules/MSource/Z2.cc \
Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \
Modules/MSolver/A2AVectors.cc \
Modules/MSolver/RBPrecCG.cc \ Modules/MSolver/RBPrecCG.cc \
Modules/MSolver/LocalCoherenceLanczos.cc \ Modules/MSolver/LocalCoherenceLanczos.cc \
Modules/MSolver/A2AVectors.cc \
Modules/MLoop/NoiseLoop.cc \
Modules/MGauge/Unit.cc \ Modules/MGauge/Unit.cc \
Modules/MGauge/Random.cc \ Modules/MGauge/UnitEm.cc \
Modules/MGauge/StochEm.cc \ Modules/MGauge/StochEm.cc \
Modules/MGauge/FundtoHirep.cc Modules/MGauge/Random.cc \
Modules/MGauge/FundtoHirep.cc \
Modules/MUtilities/TestSeqGamma.cc \
Modules/MUtilities/TestSeqConserved.cc \
Modules/MLoop/NoiseLoop.cc \
Modules/MScalar/FreeProp.cc \
Modules/MScalar/ChargedProp.cc \
Modules/MScalar/ScalarVP.cc \
Modules/MScalar/VPCounterTerms.cc \
Modules/MAction/Wilson.cc \
Modules/MAction/ZMobiusDWF.cc \
Modules/MAction/WilsonClover.cc \
Modules/MAction/DWF.cc \
Modules/MScalarSUN/TrPhi.cc \
Modules/MScalarSUN/Grad.cc \
Modules/MScalarSUN/TimeMomProbe.cc \
Modules/MScalarSUN/TrMag.cc \
Modules/MScalarSUN/TrKinetic.cc \
Modules/MScalarSUN/EMT.cc \
Modules/MScalarSUN/ShiftProbe.cc \
Modules/MScalarSUN/TransProj.cc \
Modules/MScalarSUN/StochFreeField.cc \
Modules/MScalarSUN/TwoPoint.cc \
Modules/MScalarSUN/TwoPointNPR.cc \
Modules/MScalarSUN/Div.cc \
Modules/MIO/LoadEigenPack.cc \
Modules/MIO/LoadBinary.cc \
Modules/MIO/LoadNersc.cc \
Modules/MIO/LoadCoarseEigenPack.cc
modules_hpp =\ modules_hpp =\
Modules/MSource/SeqConserved.hpp \ Modules/MContraction/Baryon.hpp \
Modules/MSource/SeqGamma.hpp \ Modules/MContraction/MesonFieldGamma.hpp \
Modules/MSource/Z2.hpp \
Modules/MSource/Point.hpp \
Modules/MSource/Wall.hpp \
Modules/MSink/Smear.hpp \
Modules/MSink/Point.hpp \
Modules/MIO/LoadBinary.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/TrPhi.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/ShiftProbe.hpp \
Modules/MScalarSUN/TimeMomProbe.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/EMT.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/DWF.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Meson.hpp \ Modules/MContraction/Meson.hpp \
Modules/MContraction/WeakHamiltonian.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \
Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/WeakNeutral4ptDisc.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/WardIdentity.hpp \ Modules/MContraction/WardIdentity.hpp \
Modules/MContraction/WeakHamiltonianEye.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \
Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/WeakHamiltonianNonEye.hpp \
Modules/MContraction/MesonFieldGmu.hpp \
Modules/MContraction/Baryon.hpp \
Modules/MContraction/WeakNeutral4ptDisc.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MUtilities/TestSeqGamma.hpp \
Modules/MFermion/FreeProp.hpp \ Modules/MFermion/FreeProp.hpp \
Modules/MFermion/GaugeProp.hpp \ Modules/MFermion/GaugeProp.hpp \
Modules/MSource/SeqGamma.hpp \
Modules/MSource/Point.hpp \
Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp \
Modules/MSource/SeqConserved.hpp \
Modules/MSink/Smear.hpp \
Modules/MSink/Point.hpp \
Modules/MSolver/A2AVectors.hpp \ Modules/MSolver/A2AVectors.hpp \
Modules/MSolver/RBPrecCG.hpp \
Modules/MSolver/LocalCoherenceLanczos.hpp \ Modules/MSolver/LocalCoherenceLanczos.hpp \
Modules/MLoop/NoiseLoop.hpp \ Modules/MSolver/RBPrecCG.hpp \
Modules/MGauge/StochEm.hpp \ Modules/MGauge/Random.hpp \
Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/FundtoHirep.hpp \
Modules/MGauge/StochEm.hpp \
Modules/MGauge/Unit.hpp \ Modules/MGauge/Unit.hpp \
Modules/MGauge/Random.hpp Modules/MGauge/UnitEm.hpp \
Modules/MUtilities/TestSeqGamma.hpp \
Modules/MUtilities/TestSeqConserved.hpp \
Modules/MLoop/NoiseLoop.hpp \
Modules/MScalar/ChargedProp.hpp \
Modules/MScalar/FreeProp.hpp \
Modules/MScalar/Scalar.hpp \
Modules/MScalar/ScalarVP.hpp \
Modules/MScalar/VPCounterTerms.hpp \
Modules/MAction/DWF.hpp \
Modules/MAction/Wilson.hpp \
Modules/MAction/WilsonClover.hpp \
Modules/MAction/ZMobiusDWF.hpp \
Modules/MScalarSUN/StochFreeField.hpp \
Modules/MScalarSUN/TwoPointNPR.hpp \
Modules/MScalarSUN/ShiftProbe.hpp \
Modules/MScalarSUN/Div.hpp \
Modules/MScalarSUN/TimeMomProbe.hpp \
Modules/MScalarSUN/TrMag.hpp \
Modules/MScalarSUN/EMT.hpp \
Modules/MScalarSUN/TwoPoint.hpp \
Modules/MScalarSUN/TrPhi.hpp \
Modules/MScalarSUN/Utils.hpp \
Modules/MScalarSUN/TransProj.hpp \
Modules/MScalarSUN/Grad.hpp \
Modules/MScalarSUN/TrKinetic.hpp \
Modules/MIO/LoadEigenPack.hpp \
Modules/MIO/LoadNersc.hpp \
Modules/MIO/LoadCoarseEigenPack.hpp \
Modules/MIO/LoadBinary.hpp

3
lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
View File

@ -57,7 +57,8 @@ void basisRotate(std::vector<Field> &basis,Eigen::MatrixXd& Qt,int j0, int j1, i
parallel_region parallel_region
{ {
std::vector < vobj > B(Nm); // Thread private
std::vector < vobj , commAllocator<vobj> > B(Nm); // Thread private
parallel_for_internal(int ss=0;ss < grid->oSites();ss++){ parallel_for_internal(int ss=0;ss < grid->oSites();ss++){
for(int j=j0; j<j1; ++j) B[j]=0.; for(int j=j0; j<j1; ++j) B[j]=0.;

11
lib/lattice/Lattice_rng.h
View File

@ -158,10 +158,19 @@ namespace Grid {
// tens of seconds per trajectory so this is clean in all reasonable cases, // tens of seconds per trajectory so this is clean in all reasonable cases,
// and margin of safety is orders of magnitude. // and margin of safety is orders of magnitude.
// We could hack Sitmo to skip in the higher order words of state if necessary // We could hack Sitmo to skip in the higher order words of state if necessary
//
// Replace with 2^30 ; avoid problem on large volumes
//
///////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////
// uint64_t skip = site+1; // Old init Skipped then drew. Checked compat with faster init // uint64_t skip = site+1; // Old init Skipped then drew. Checked compat with faster init
const int shift = 30;
uint64_t skip = site; uint64_t skip = site;
skip = skip<<40;
skip = skip<<shift;
assert((skip >> shift)==site); // check for overflow
eng.discard(skip); eng.discard(skip);
// std::cout << " Engine " <<site << " state " <<eng<<std::endl; // std::cout << " Engine " <<site << " state " <<eng<<std::endl;
} }

167
lib/qcd/action/gauge/Photon.h
View File

@ -58,9 +58,12 @@ namespace QCD{
public: public:
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3); GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3);
GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2); GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2, qedInf, 3);
public: public:
Photon(Gauge gauge, ZmScheme zmScheme); Photon(Gauge gauge, ZmScheme zmScheme);
Photon(Gauge gauge, ZmScheme zmScheme, std::vector<Real> improvements);
Photon(Gauge gauge, ZmScheme zmScheme, Real G0);
Photon(Gauge gauge, ZmScheme zmScheme, std::vector<Real> improvements, Real G0);
virtual ~Photon(void) = default; virtual ~Photon(void) = default;
void FreePropagator(const GaugeField &in, GaugeField &out); void FreePropagator(const GaugeField &in, GaugeField &out);
void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out);
@ -68,19 +71,42 @@ namespace QCD{
void StochasticField(GaugeField &out, GridParallelRNG &rng); void StochasticField(GaugeField &out, GridParallelRNG &rng);
void StochasticField(GaugeField &out, GridParallelRNG &rng, void StochasticField(GaugeField &out, GridParallelRNG &rng,
const GaugeLinkField &weight); const GaugeLinkField &weight);
void UnitField(GaugeField &out);
private: private:
void infVolPropagator(GaugeLinkField &out);
void invKHatSquared(GaugeLinkField &out); void invKHatSquared(GaugeLinkField &out);
void zmSub(GaugeLinkField &out); void zmSub(GaugeLinkField &out);
private: private:
Gauge gauge_; Gauge gauge_;
ZmScheme zmScheme_; ZmScheme zmScheme_;
std::vector<Real> improvement_;
Real G0_;
}; };
typedef Photon<QedGimplR> PhotonR; typedef Photon<QedGimplR> PhotonR;
template<class Gimpl> template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme) Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme)
: gauge_(gauge), zmScheme_(zmScheme) : gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector<Real>()),
G0_(0.15493339023106021408483720810737508876916113364521)
{}
template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme,
std::vector<Real> improvements)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements),
G0_(0.15493339023106021408483720810737508876916113364521)
{}
template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme, Real G0)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector<Real>()), G0_(G0)
{}
template<class Gimpl>
Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme,
std::vector<Real> improvements, Real G0)
: gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), G0_(G0)
{} {}
template<class Gimpl> template<class Gimpl>
@ -96,6 +122,34 @@ namespace QCD{
theFFT.FFT_all_dim(out,prop_k,FFT::backward); theFFT.FFT_all_dim(out,prop_k,FFT::backward);
} }
template<class Gimpl>
void Photon<Gimpl>::infVolPropagator(GaugeLinkField &out)
{
auto *grid = dynamic_cast<GridCartesian *>(out._grid);
LatticeReal xmu(grid);
GaugeLinkField one(grid);
const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
std::vector<int> x0(nd,0);
TComplex Tone = Complex(1.0,0.0);
TComplex Tzero = Complex(G0_,0.0);
FFT fft(grid);
one = Complex(1.0,0.0);
out = zero;
for(int mu = 0; mu < nd; mu++)
{
LatticeCoordinate(xmu,mu);
Real lo2 = l[mu]/2.0;
xmu = where(xmu < lo2, xmu, xmu-double(l[mu]));
out = out + toComplex(4*M_PI*M_PI*xmu*xmu);
}
pokeSite(Tone, out, x0);
out = one/out;
pokeSite(Tzero, out, x0);
fft.FFT_all_dim(out, out, FFT::forward);
}
template<class Gimpl> template<class Gimpl>
void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out) void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out)
{ {
@ -127,6 +181,7 @@ namespace QCD{
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
const unsigned int nd = grid->_ndimension; const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
switch (zmScheme_) switch (zmScheme_)
{ {
@ -148,11 +203,17 @@ namespace QCD{
for(int d = 0; d < grid->_ndimension - 1; d++) for(int d = 0; d < grid->_ndimension - 1; d++)
{ {
LatticeCoordinate(coor,d); LatticeCoordinate(coor,d);
coor = where(coor < Integer(l[d]/2), coor, coor-Integer(l[d]));
spNrm = spNrm + coor*coor; spNrm = spNrm + coor*coor;
} }
out = where(spNrm == Integer(0), 0.*out, out); out = where(spNrm == Integer(0), 0.*out, out);
break; // IR improvement
for(int i = 0; i < improvement_.size(); i++)
{
Real f = sqrt(improvement_[i]+1);
out = where(spNrm == Integer(i+1), f*out, out);
}
} }
default: default:
break; break;
@ -163,13 +224,28 @@ namespace QCD{
void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in, void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in,
GaugeField &out) GaugeField &out)
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
LatticeComplex k2Inv(grid); LatticeComplex momProp(grid);
invKHatSquared(k2Inv); switch (zmScheme_)
zmSub(k2Inv); {
case ZmScheme::qedTL:
case ZmScheme::qedL:
{
invKHatSquared(momProp);
zmSub(momProp);
break;
}
case ZmScheme::qedInf:
{
infVolPropagator(momProp);
break;
}
default:
break;
}
out = in*k2Inv; out = in*momProp;
} }
template<class Gimpl> template<class Gimpl>
@ -179,14 +255,30 @@ namespace QCD{
const unsigned int nd = grid->_ndimension; const unsigned int nd = grid->_ndimension;
std::vector<int> latt_size = grid->_fdimensions; std::vector<int> latt_size = grid->_fdimensions;
Integer vol = 1; switch (zmScheme_)
for(int d = 0; d < nd; d++)
{ {
vol = vol * latt_size[d]; case ZmScheme::qedTL:
case ZmScheme::qedL:
{
Integer vol = 1;
for(int d = 0; d < nd; d++)
{
vol = vol * latt_size[d];
}
invKHatSquared(weight);
weight = sqrt(vol)*sqrt(weight);
zmSub(weight);
break;
}
case ZmScheme::qedInf:
{
infVolPropagator(weight);
weight = sqrt(real(weight));
break;
}
default:
break;
} }
invKHatSquared(weight);
weight = sqrt(vol)*sqrt(weight);
zmSub(weight);
} }
template<class Gimpl> template<class Gimpl>
@ -209,13 +301,52 @@ namespace QCD{
GaugeField aTilde(grid); GaugeField aTilde(grid);
FFT fft(grid); FFT fft(grid);
switch (zmScheme_)
{
case ZmScheme::qedTL:
case ZmScheme::qedL:
{
for(int mu = 0; mu < nd; mu++)
{
gaussian(rng, r);
r = weight*r;
pokeLorentz(aTilde, r, mu);
}
break;
}
case ZmScheme::qedInf:
{
Complex shift(1., 1.); // This needs to be a GaugeLink element?
for(int mu = 0; mu < nd; mu++)
{
bernoulli(rng, r);
r = weight*(2.*r - shift);
pokeLorentz(aTilde, r, mu);
}
break;
}
default:
break;
}
fft.FFT_all_dim(out, aTilde, FFT::backward);
out = real(out);
}
template<class Gimpl>
void Photon<Gimpl>::UnitField(GaugeField &out)
{
auto *grid = dynamic_cast<GridCartesian *>(out._grid);
const unsigned int nd = grid->_ndimension;
GaugeLinkField r(grid);
r = Complex(1.0,0.0);
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
gaussian(rng, r); pokeLorentz(out, r, mu);
r = weight*r;
pokeLorentz(aTilde, r, mu);
} }
fft.FFT_all_dim(out, aTilde, FFT::backward);
out = real(out); out = real(out);
} }

5
lib/qcd/hmc/checkpointers/ScidacCheckpointer.h
View File

@ -86,10 +86,7 @@ class ScidacHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
_ScidacWriter.writeScidacFieldRecord(U, MData); _ScidacWriter.writeScidacFieldRecord(U, MData);
_ScidacWriter.close(); _ScidacWriter.close();
std::cout << GridLogMessage << "Written Scidac Configuration on " << config std::cout << GridLogMessage << "Written Scidac Configuration on " << config << std::endl;
<< " checksum " << std::hex << nersc_csum<<"/"
<< scidac_csuma<<"/" << scidac_csumb
<< std::dec << std::endl;
} }
}; };