portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 17:25:37 +01:00
Code Releases Activity

QedFVol: add ChargedProp as an input to ScalarVP module, instead of calculating scalar propagator within ScalarVP.

Browse Source
This commit is contained in:
James Harrison 2017-06-08 17:43:39 +01:00
parent e38612e6fa
commit 20ac13fdf3
4 changed files with 51 additions and 249 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -23,8 +23,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(), getName()+"_0", getName()+"_D1", std::vector<std::string> out = {getName(), getName()+"_Q",
getName()+"_D1D1", getName()+"_D2"}; getName()+"_Sun", getName()+"_Tad"};
return out; return out;
} }
@ -40,9 +40,9 @@ void TChargedProp::setup(void)
} }
GFSrcName_ = "_" + getName() + "_DinvSrc"; GFSrcName_ = "_" + getName() + "_DinvSrc";
prop0Name_ = getName() + "_0"; prop0Name_ = getName() + "_0";
propD1Name_ = getName() + "_D1"; propQName_ = getName() + "_Q";
propD1D1Name_ = getName() + "_D1D1"; propSunName_ = getName() + "_Sun";
propD2Name_ = getName() + "_D2"; propTadName_ = getName() + "_Tad";
if (!env().hasRegisteredObject(freeMomPropName_)) if (!env().hasRegisteredObject(freeMomPropName_))
{ {
env().registerLattice<ScalarField>(freeMomPropName_); env().registerLattice<ScalarField>(freeMomPropName_);
@ -63,9 +63,9 @@ void TChargedProp::setup(void)
env().registerLattice<ScalarField>(prop0Name_); env().registerLattice<ScalarField>(prop0Name_);
} }
env().registerLattice<ScalarField>(getName()); env().registerLattice<ScalarField>(getName());
env().registerLattice<ScalarField>(propD1Name_); env().registerLattice<ScalarField>(propQName_);
env().registerLattice<ScalarField>(propD1D1Name_); env().registerLattice<ScalarField>(propSunName_);
env().registerLattice<ScalarField>(propD2Name_); env().registerLattice<ScalarField>(propTadName_);
} }
// execution /////////////////////////////////////////////////////////////////// // execution ///////////////////////////////////////////////////////////////////
@ -140,9 +140,9 @@ void TChargedProp::execute(void)
<< ", charge= " << par().charge << ")..." << std::endl; << ", charge= " << par().charge << ")..." << std::endl;
ScalarField &prop = *env().createLattice<ScalarField>(getName()); ScalarField &prop = *env().createLattice<ScalarField>(getName());
ScalarField &propD1 = *env().createLattice<ScalarField>(propD1Name_); ScalarField &propQ = *env().createLattice<ScalarField>(propQName_);
ScalarField &propD1D1 = *env().createLattice<ScalarField>(propD1D1Name_); ScalarField &propSun = *env().createLattice<ScalarField>(propSunName_);
ScalarField &propD2 = *env().createLattice<ScalarField>(propD2Name_); ScalarField &propTad = *env().createLattice<ScalarField>(propTadName_);
ScalarField buf(env().getGrid()); ScalarField buf(env().getGrid());
ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_; ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_;
double q = par().charge; double q = par().charge;
@ -151,22 +151,22 @@ void TChargedProp::execute(void)
buf = GFSrc; buf = GFSrc;
momD1(buf, fft); momD1(buf, fft);
buf = -G*buf; buf = -G*buf;
fft.FFT_all_dim(propD1, buf, FFT::backward); fft.FFT_all_dim(propQ, buf, FFT::backward);
// G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) // G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src)
buf = -buf; buf = -buf;
momD1(buf, fft); momD1(buf, fft);
propD1D1 = G*buf; propSun = G*buf;
fft.FFT_all_dim(propD1D1, propD1D1, FFT::backward); fft.FFT_all_dim(propSun, propSun, FFT::backward);
// -G*momD2*G*F*Src (momD2 = F*D2*Finv) // -G*momD2*G*F*Src (momD2 = F*D2*Finv)
buf = GFSrc; buf = GFSrc;
momD2(buf, fft); momD2(buf, fft);
buf = -G*buf; buf = -G*buf;
fft.FFT_all_dim(propD2, buf, FFT::backward); fft.FFT_all_dim(propTad, buf, FFT::backward);
// full charged scalar propagator // full charged scalar propagator
prop = (*prop0_) + q*propD1 + q*q*propD1D1 + q*q*propD2; prop = (*prop0_) + q*propQ + q*q*propSun + q*q*propTad;
// OUTPUT IF NECESSARY // OUTPUT IF NECESSARY
if (!par().output.empty()) if (!par().output.empty())
@ -200,29 +200,29 @@ void TChargedProp::execute(void)
} }
write(writer, "prop_0", result); write(writer, "prop_0", result);
// Write propagator D1 term // Write propagator O(q) term
sliceSum(propD1, vecBuf, Tp); sliceSum(propQ, vecBuf, Tp);
for (unsigned int t = 0; t < vecBuf.size(); ++t) for (unsigned int t = 0; t < vecBuf.size(); ++t)
{ {
result[t] = TensorRemove(vecBuf[t]); result[t] = TensorRemove(vecBuf[t]);
} }
write(writer, "prop_D1", result); write(writer, "prop_Q", result);
// Write propagator D1D1 term // Write propagator sunset term
sliceSum(propD1D1, vecBuf, Tp); sliceSum(propSun, vecBuf, Tp);
for (unsigned int t = 0; t < vecBuf.size(); ++t) for (unsigned int t = 0; t < vecBuf.size(); ++t)
{ {
result[t] = TensorRemove(vecBuf[t]); result[t] = TensorRemove(vecBuf[t]);
} }
write(writer, "prop_D1D1", result); write(writer, "prop_Sun", result);
// Write propagator D2 term // Write propagator tadpole term
sliceSum(propD2, vecBuf, Tp); sliceSum(propTad, vecBuf, Tp);
for (unsigned int t = 0; t < vecBuf.size(); ++t) for (unsigned int t = 0; t < vecBuf.size(); ++t)
{ {
result[t] = TensorRemove(vecBuf[t]); result[t] = TensorRemove(vecBuf[t]);
} }
write(writer, "prop_D2", result); write(writer, "prop_Tad", result);
} }
} }

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

@ -46,7 +46,7 @@ private:
void momD2(ScalarField &s, FFT &fft); void momD2(ScalarField &s, FFT &fft);
private: private:
std::string freeMomPropName_, GFSrcName_, prop0Name_, std::string freeMomPropName_, GFSrcName_, prop0Name_,
propD1Name_, propD1D1Name_, propD2Name_; propQName_, propSunName_, propTadName_;
std::vector<std::string> phaseName_; std::vector<std::string> phaseName_;
ScalarField *freeMomProp_, *GFSrc_, *prop0_; ScalarField *freeMomProp_, *GFSrc_, *prop0_;
std::vector<ScalarField *> phase_; std::vector<ScalarField *> phase_;

240
extras/Hadrons/Modules/MScalar/ScalarVP.cc
View File

@ -1,3 +1,4 @@
#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp> #include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp> #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
@ -16,16 +17,19 @@ TScalarVP::TScalarVP(const std::string name)
// dependencies/products /////////////////////////////////////////////////////// // dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TScalarVP::getInput(void) std::vector<std::string> TScalarVP::getInput(void)
{ {
std::vector<std::string> in = {par().source, par().emField}; propQName_ = par().scalarProp + "_Q";
propSunName_ = par().scalarProp + "_Sun";
propTadName_ = par().scalarProp + "_Tad";
std::vector<std::string> in = {par().emField, propQName_, propSunName_,
propTadName_};
return in; return in;
} }
std::vector<std::string> TScalarVP::getOutput(void) std::vector<std::string> TScalarVP::getOutput(void)
{ {
std::vector<std::string> out = {getName()+"_propQ", std::vector<std::string> out;
getName()+"_propSun",
getName()+"_propTad"};
for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{ {
@ -44,12 +48,9 @@ std::vector<std::string> TScalarVP::getOutput(void)
// setup /////////////////////////////////////////////////////////////////////// // setup ///////////////////////////////////////////////////////////////////////
void TScalarVP::setup(void) void TScalarVP::setup(void)
{ {
freeMomPropName_ = FREEMOMPROP(par().mass); freeMomPropName_ = FREEMOMPROP(static_cast<TChargedProp *>(env().getModule(par().scalarProp))->par().mass);
GFSrcName_ = "_" + getName() + "_DinvSrc"; GFSrcName_ = "_" + par().scalarProp + "_DinvSrc";
prop0Name_ = getName() + "_prop0"; prop0Name_ = par().scalarProp + "_0";
propQName_ = getName() + "_propQ";
propSunName_ = getName() + "_propSun";
propTadName_ = getName() + "_propTad";
phaseName_.clear(); phaseName_.clear();
muPropQName_.clear(); muPropQName_.clear();
@ -74,174 +75,38 @@ void TScalarVP::setup(void)
freeVpTensorName_.push_back(freeVpTensorName_mu); freeVpTensorName_.push_back(freeVpTensorName_mu);
} }
if (!env().hasRegisteredObject(freeMomPropName_))
{
env().registerLattice<ScalarField>(freeMomPropName_);
}
if (!env().hasRegisteredObject(phaseName_[0]))
{
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
env().registerLattice<ScalarField>(phaseName_[mu]);
}
}
if (!env().hasRegisteredObject(GFSrcName_))
{
env().registerLattice<ScalarField>(GFSrcName_);
}
if (!env().hasRegisteredObject(prop0Name_))
{
env().registerLattice<ScalarField>(prop0Name_);
}
env().registerLattice<ScalarField>(propQName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{ {
env().registerLattice<ScalarField>(muPropQName_[mu]); env().registerLattice<ScalarField>(muPropQName_[mu]);
}
env().registerLattice<ScalarField>(propSunName_);
env().registerLattice<ScalarField>(propTadName_);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
for (unsigned int nu = 0; nu < env().getNd(); ++nu) for (unsigned int nu = 0; nu < env().getNd(); ++nu)
{ {
env().registerLattice<ScalarField>(vpTensorName_[mu][nu]); env().registerLattice<ScalarField>(vpTensorName_[mu][nu]);
env().registerLattice<ScalarField>(freeVpTensorName_[mu][nu]); env().registerLattice<ScalarField>(freeVpTensorName_[mu][nu]);
} }
} }
} }
// execution /////////////////////////////////////////////////////////////////// // execution ///////////////////////////////////////////////////////////////////
void TScalarVP::execute(void) void TScalarVP::execute(void)
{ {
// CACHING ANALYTIC EXPRESSIONS // Get objects cached by ChargedProp module
ScalarField &source = *env().getObject<ScalarField>(par().source);
Complex ci(0.0,1.0); Complex ci(0.0,1.0);
FFT fft(env().getGrid()); FFT fft(env().getGrid());
Real q = par().charge; Real q = static_cast<TChargedProp *>(env().getModule(par().scalarProp))->par().charge;
// cache momentum-space free scalar propagator freeMomProp_ = env().getObject<ScalarField>(freeMomPropName_);
if (!env().hasCreatedObject(freeMomPropName_)) for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{ {
LOG(Message) << "Caching momentum space free scalar propagator" phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
<< " (mass= " << par().mass << ")..." << std::endl;
freeMomProp_ = env().createLattice<ScalarField>(freeMomPropName_);
SIMPL::MomentumSpacePropagator(*freeMomProp_, par().mass);
} }
else GFSrc_ = env().getObject<ScalarField>(GFSrcName_);
{ prop0_ = env().getObject<ScalarField>(prop0Name_);
freeMomProp_ = env().getObject<ScalarField>(freeMomPropName_);
}
// cache phases
if (!env().hasCreatedObject(phaseName_[0]))
{
std::vector<int> &l = env().getGrid()->_fdimensions;
LOG(Message) << "Caching shift phases..." << std::endl;
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Real twoPiL = M_PI*2./l[mu];
phase_.push_back(env().createLattice<ScalarField>(phaseName_[mu]));
LatticeCoordinate(*(phase_[mu]), mu);
*(phase_[mu]) = exp(ci*twoPiL*(*(phase_[mu])));
}
}
else
{
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
}
}
// cache G*F*src
if (!env().hasCreatedObject(GFSrcName_))
{
GFSrc_ = env().createLattice<ScalarField>(GFSrcName_);
fft.FFT_all_dim(*GFSrc_, source, FFT::forward);
*GFSrc_ = (*freeMomProp_)*(*GFSrc_);
}
else
{
GFSrc_ = env().getObject<ScalarField>(GFSrcName_);
}
// cache position-space free scalar propagators
if (!env().hasCreatedObject(prop0Name_))
{
prop0_ = env().createLattice<ScalarField>(prop0Name_);
fft.FFT_all_dim(*prop0_, *GFSrc_, FFT::backward);
}
else
{
prop0_ = env().getObject<ScalarField>(prop0Name_);
}
// PROPAGATOR CALCULATION
// Propagator from unshifted source // Propagator from unshifted source
LOG(Message) << "Computing O(alpha) charged scalar propagator" ScalarField &propQ = *env().getObject<ScalarField>(propQName_);
<< " (mass= " << par().mass ScalarField &propSun = *env().getObject<ScalarField>(propSunName_);
<< ", charge= " << q << ")..." ScalarField &propTad = *env().getObject<ScalarField>(propTadName_);
<< std::endl;
ScalarField &propQ = *env().createLattice<ScalarField>(propQName_);
ScalarField &propSun = *env().createLattice<ScalarField>(propSunName_);
ScalarField &propTad = *env().createLattice<ScalarField>(propTadName_);
chargedProp(propQ, propSun, propTad, *GFSrc_, fft);
// // OUTPUT IF NECESSARY
// if (!par().output.empty())
// {
// ScalarField fullProp = (*prop0_) + q*propQ + q*q*propSun + q*q*propTad;
// std::string filename = par().output + "_prop_000." +
// std::to_string(env().getTrajectory());
// LOG(Message) << "Saving zero-momentum projection to '"
// << filename << "'..." << std::endl;
// CorrWriter writer(filename);
// std::vector<TComplex> vecBuf;
// std::vector<Complex> result;
// write(writer, "charge", q);
// // Write full propagator
// sliceSum(fullProp, vecBuf, Tp);
// result.resize(vecBuf.size());
// for (unsigned int t = 0; t < vecBuf.size(); ++t)
// {
// result[t] = TensorRemove(vecBuf[t]);
// }
// write(writer, "prop", result);
// // Write free propagator
// sliceSum(*prop0_, vecBuf, Tp);
// for (unsigned int t = 0; t < vecBuf.size(); ++t)
// {
// result[t] = TensorRemove(vecBuf[t]);
// }
// write(writer, "prop_0", result);
// // Write propagator D1 term
// sliceSum(propD1, vecBuf, Tp);
// for (unsigned int t = 0; t < vecBuf.size(); ++t)
// {
// result[t] = TensorRemove(vecBuf[t]);
// }
// write(writer, "prop_q", result);
// // Write propagator D1D1 term
// sliceSum(propD1D1, vecBuf, Tp);
// for (unsigned int t = 0; t < vecBuf.size(); ++t)
// {
// result[t] = TensorRemove(vecBuf[t]);
// }
// write(writer, "prop_sun", result);
// // Write propagator D2 term
// sliceSum(propD2, vecBuf, Tp);
// for (unsigned int t = 0; t < vecBuf.size(); ++t)
// {
// result[t] = TensorRemove(vecBuf[t]);
// }
// write(writer, "prop_tad", result);
// }
// Propagators from shifted sources // Propagators from shifted sources
LOG(Message) << "Computing O(q) charged scalar propagators..." LOG(Message) << "Computing O(q) charged scalar propagators..."
@ -357,7 +222,7 @@ void TScalarVP::execute(void)
std::vector<Complex> result; std::vector<Complex> result;
write(writer, "charge", q); write(writer, "charge", q);
write(writer, "mass", par().mass); write(writer, "mass", static_cast<TChargedProp *>(env().getModule(par().scalarProp))->par().mass);
for (unsigned int mu = 0; mu < env().getNd(); ++mu) for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{ {
@ -386,34 +251,6 @@ void TScalarVP::execute(void)
} }
} }
// Calculate O(q) and O(q^2) terms of position-space charged propagator
void TScalarVP::chargedProp(ScalarField &prop_q, ScalarField &prop_sun,
ScalarField &prop_tad, ScalarField &GFSrc,
FFT &fft)
{
Complex ci(0.0,1.0);
ScalarField &G = *freeMomProp_;
ScalarField buf(env().getGrid());
// -G*momD1*G*F*Src (momD1 = F*D1*Finv)
buf = GFSrc;
momD1(buf, fft);
buf = G*buf;
prop_q = -buf;
fft.FFT_all_dim(prop_q, prop_q, FFT::backward);
// G*momD1*G*momD1*G*F*Src
momD1(buf, fft);
prop_sun = G*buf;
fft.FFT_all_dim(prop_sun, prop_sun, FFT::backward);
// -G*momD2*G*F*Src (momD2 = F*D2*Finv)
buf = GFSrc;
momD2(buf, fft);
prop_tad = -G*buf;
fft.FFT_all_dim(prop_tad, prop_tad, FFT::backward);
}
void TScalarVP::momD1(ScalarField &s, FFT &fft) void TScalarVP::momD1(ScalarField &s, FFT &fft)
{ {
EmField &A = *env().getObject<EmField>(par().emField); EmField &A = *env().getObject<EmField>(par().emField);
@ -443,32 +280,3 @@ void TScalarVP::momD1(ScalarField &s, FFT &fft)
s = result; s = result;
} }
void TScalarVP::momD2(ScalarField &s, FFT &fft)
{
EmField &A = *env().getObject<EmField>(par().emField);
ScalarField buf(env().getGrid()), result(env().getGrid()),
Amu(env().getGrid());
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*Amu*buf;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result + .5*buf;
}
fft.FFT_all_dim(s, s, FFT::backward);
for (unsigned int mu = 0; mu < env().getNd(); ++mu)
{
Amu = peekLorentz(A, mu);
buf = Amu*Amu*s;
fft.FFT_all_dim(buf, buf, FFT::forward);
result = result + .5*adj(*phase_[mu])*buf;
}
s = result;
}

8
extras/Hadrons/Modules/MScalar/ScalarVP.hpp
View File

@ -17,9 +17,7 @@ class ScalarVPPar: Serializable
public: public:
GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar, GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar,
std::string, emField, std::string, emField,
std::string, source, std::string, scalarProp,
double, mass,
double, charge,
std::string, output); std::string, output);
}; };
@ -42,11 +40,7 @@ public:
// execution // execution
virtual void execute(void); virtual void execute(void);
private: private:
void chargedProp(ScalarField &prop_q, ScalarField &prop_sun,
ScalarField &prop_tad, ScalarField &GFSrc,
FFT &fft);
void momD1(ScalarField &s, FFT &fft); void momD1(ScalarField &s, FFT &fft);
void momD2(ScalarField &s, FFT &fft);
private: private:
std::string freeMomPropName_, GFSrcName_, std::string freeMomPropName_, GFSrcName_,
prop0Name_, propQName_, prop0Name_, propQName_,