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

Output all terms of scalar propagator separately

Browse Source
This commit is contained in:
James Harrison 2017-03-24 17:13:55 +00:00
parent 0c006fbfaa
commit 85516e9c7c
2 changed files with 78 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -23,7 +23,8 @@ std::vector<std::string> TChargedProp::getInput(void)
std::vector<std::string> TChargedProp::getOutput(void)
{
std::vector<std::string> out = {getName()};
std::vector<std::string> out = {getName(), getName()+"_0", getName()+"_D1",
getName()+"_D1D1", getName()+"_D2"};
return out;
}
@ -38,6 +39,10 @@ void TChargedProp::setup(void)
phaseName_.push_back("_shiftphase_" + std::to_string(mu));
}
GFSrcName_ = "_" + getName() + "_DinvSrc";
prop0Name_ = getName() + "_0";
propD1Name_ = getName() + "_D1";
propD1D1Name_ = getName() + "_D1D1";
propD2Name_ = getName() + "_D2";
if (!env().hasRegisteredObject(freeMomPropName_))
{
env().registerLattice<ScalarField>(freeMomPropName_);
@ -53,7 +58,14 @@ void TChargedProp::setup(void)
{
env().registerLattice<ScalarField>(GFSrcName_);
}
if (!env().hasRegisteredObject(prop0Name_))
{
env().registerLattice<ScalarField>(prop0Name_);
}
env().registerLattice<ScalarField>(getName());
env().registerLattice<ScalarField>(propD1Name_);
env().registerLattice<ScalarField>(propD1D1Name_);
env().registerLattice<ScalarField>(propD2Name_);
}
// execution ///////////////////////////////////////////////////////////////////
@ -64,7 +76,7 @@ void TChargedProp::execute(void)
Complex ci(0.0,1.0);
FFT fft(env().getGrid());
// cache free scalar propagator
// cache momentum-space free scalar propagator
if (!env().hasCreatedObject(freeMomPropName_))
{
LOG(Message) << "Caching momentum space free scalar propagator"
@ -88,6 +100,17 @@ void TChargedProp::execute(void)
{
GFSrc_ = env().getObject<ScalarField>(GFSrcName_);
}
// cache free scalar propagator
if (!env().hasCreatedObject(prop0Name_))
{
prop0_ = env().createLattice<ScalarField>(prop0Name_);
*prop0_ = *GFSrc_;
fft.FFT_all_dim(*prop0_, *prop0_, FFT::backward);
}
else
{
prop0_ = env().getObject<ScalarField>(prop0Name_);
}
// cache phases
if (!env().hasCreatedObject(phaseName_[0]))
{
@ -117,30 +140,33 @@ void TChargedProp::execute(void)
<< ", charge= " << par().charge << ")..." << std::endl;
ScalarField &prop = *env().createLattice<ScalarField>(getName());
ScalarField &propD1 = *env().createLattice<ScalarField>(propD1Name_);
ScalarField &propD1D1 = *env().createLattice<ScalarField>(propD1D1Name_);
ScalarField &propD2 = *env().createLattice<ScalarField>(propD2Name_);
ScalarField buf(env().getGrid());
ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_;
double q = par().charge;
// G*F*Src
prop = GFSrc;
// - q*G*momD1*G*F*Src (momD1 = F*D1*Finv)
// -G*momD1*G*F*Src (momD1 = F*D1*Finv)
buf = GFSrc;
momD1(buf, fft);
buf = G*buf;
prop = prop - q*buf;
buf = -G*buf;
fft.FFT_all_dim(propD1, buf, FFT::backward);
// + q^2*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;
momD1(buf, fft);
prop = prop + q*q*G*buf;
propD1D1 = G*buf;
fft.FFT_all_dim(propD1D1, propD1D1, FFT::backward);
// - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv)
// -G*momD2*G*F*Src (momD2 = F*D2*Finv)
buf = GFSrc;
momD2(buf, fft);
prop = prop - q*q*G*buf;
buf = -G*buf;
fft.FFT_all_dim(propD2, buf, FFT::backward);
// final FT
fft.FFT_all_dim(prop, prop, FFT::backward);
// full charged scalar propagator
prop = (*prop0_) + q*propD1 + q*q*propD1D1 + q*q*propD2;
// OUTPUT IF NECESSARY
if (!par().output.empty())
@ -155,14 +181,48 @@ void TChargedProp::execute(void)
std::vector<TComplex> vecBuf;
std::vector<Complex> result;
write(writer, "charge", q);
// Write full propagator
sliceSum(prop, vecBuf, Tp);
result.resize(vecBuf.size());
for (unsigned int t = 0; t < vecBuf.size(); ++t)
{
result[t] = TensorRemove(vecBuf[t]);
}
write(writer, "charge", q);
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_D1", 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_D1D1", 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_D2", result);
}
}

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

@ -45,9 +45,10 @@ private:
void momD1(ScalarField &s, FFT &fft);
void momD2(ScalarField &s, FFT &fft);
private:
std::string freeMomPropName_, GFSrcName_;
std::string freeMomPropName_, GFSrcName_, prop0Name_,
propD1Name_, propD1D1Name_, propD2Name_;
std::vector<std::string> phaseName_;
ScalarField *freeMomProp_, *GFSrc_;
ScalarField *freeMomProp_, *GFSrc_, *prop0_;
std::vector<ScalarField *> phase_;
EmField *A;
};