diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index be64b5ec..a9089056 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -23,8 +23,8 @@ std::vector TChargedProp::getInput(void) std::vector TChargedProp::getOutput(void) { - std::vector out = {getName(), getName()+"_0", getName()+"_D1", - getName()+"_D1D1", getName()+"_D2"}; + std::vector out = {getName(), getName()+"_Q", + getName()+"_Sun", getName()+"_Tad"}; return out; } @@ -40,9 +40,9 @@ void TChargedProp::setup(void) } GFSrcName_ = "_" + getName() + "_DinvSrc"; prop0Name_ = getName() + "_0"; - propD1Name_ = getName() + "_D1"; - propD1D1Name_ = getName() + "_D1D1"; - propD2Name_ = getName() + "_D2"; + propQName_ = getName() + "_Q"; + propSunName_ = getName() + "_Sun"; + propTadName_ = getName() + "_Tad"; if (!env().hasRegisteredObject(freeMomPropName_)) { env().registerLattice(freeMomPropName_); @@ -63,9 +63,9 @@ void TChargedProp::setup(void) env().registerLattice(prop0Name_); } env().registerLattice(getName()); - env().registerLattice(propD1Name_); - env().registerLattice(propD1D1Name_); - env().registerLattice(propD2Name_); + env().registerLattice(propQName_); + env().registerLattice(propSunName_); + env().registerLattice(propTadName_); } // execution /////////////////////////////////////////////////////////////////// @@ -140,9 +140,9 @@ void TChargedProp::execute(void) << ", charge= " << par().charge << ")..." << std::endl; ScalarField &prop = *env().createLattice(getName()); - ScalarField &propD1 = *env().createLattice(propD1Name_); - ScalarField &propD1D1 = *env().createLattice(propD1D1Name_); - ScalarField &propD2 = *env().createLattice(propD2Name_); + ScalarField &propQ = *env().createLattice(propQName_); + ScalarField &propSun = *env().createLattice(propSunName_); + ScalarField &propTad = *env().createLattice(propTadName_); ScalarField buf(env().getGrid()); ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_; double q = par().charge; @@ -151,22 +151,22 @@ void TChargedProp::execute(void) buf = GFSrc; momD1(buf, fft); 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) buf = -buf; momD1(buf, fft); - propD1D1 = G*buf; - fft.FFT_all_dim(propD1D1, propD1D1, FFT::backward); + propSun = G*buf; + fft.FFT_all_dim(propSun, propSun, FFT::backward); // -G*momD2*G*F*Src (momD2 = F*D2*Finv) buf = GFSrc; momD2(buf, fft); buf = -G*buf; - fft.FFT_all_dim(propD2, buf, FFT::backward); + fft.FFT_all_dim(propTad, buf, FFT::backward); // 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 if (!par().output.empty()) @@ -200,29 +200,29 @@ void TChargedProp::execute(void) } write(writer, "prop_0", result); - // Write propagator D1 term - sliceSum(propD1, vecBuf, Tp); + // Write propagator O(q) term + sliceSum(propQ, vecBuf, Tp); for (unsigned int t = 0; t < vecBuf.size(); ++t) { result[t] = TensorRemove(vecBuf[t]); } - write(writer, "prop_D1", result); + write(writer, "prop_Q", result); - // Write propagator D1D1 term - sliceSum(propD1D1, vecBuf, Tp); + // Write propagator sunset term + sliceSum(propSun, vecBuf, Tp); for (unsigned int t = 0; t < vecBuf.size(); ++t) { result[t] = TensorRemove(vecBuf[t]); } - write(writer, "prop_D1D1", result); + write(writer, "prop_Sun", result); - // Write propagator D2 term - sliceSum(propD2, vecBuf, Tp); + // Write propagator tadpole term + sliceSum(propTad, vecBuf, Tp); for (unsigned int t = 0; t < vecBuf.size(); ++t) { result[t] = TensorRemove(vecBuf[t]); } - write(writer, "prop_D2", result); + write(writer, "prop_Tad", result); } } diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index aeb92179..369fff30 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -46,7 +46,7 @@ private: void momD2(ScalarField &s, FFT &fft); private: std::string freeMomPropName_, GFSrcName_, prop0Name_, - propD1Name_, propD1D1Name_, propD2Name_; + propQName_, propSunName_, propTadName_; std::vector phaseName_; ScalarField *freeMomProp_, *GFSrc_, *prop0_; std::vector phase_; diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.cc b/extras/Hadrons/Modules/MScalar/ScalarVP.cc index c91b98ae..7a3b4f9e 100644 --- a/extras/Hadrons/Modules/MScalar/ScalarVP.cc +++ b/extras/Hadrons/Modules/MScalar/ScalarVP.cc @@ -1,3 +1,4 @@ +#include #include #include @@ -16,16 +17,19 @@ TScalarVP::TScalarVP(const std::string name) // dependencies/products /////////////////////////////////////////////////////// std::vector TScalarVP::getInput(void) { - std::vector in = {par().source, par().emField}; + propQName_ = par().scalarProp + "_Q"; + propSunName_ = par().scalarProp + "_Sun"; + propTadName_ = par().scalarProp + "_Tad"; + + std::vector in = {par().emField, propQName_, propSunName_, + propTadName_}; return in; } std::vector TScalarVP::getOutput(void) { - std::vector out = {getName()+"_propQ", - getName()+"_propSun", - getName()+"_propTad"}; + std::vector out; for (unsigned int mu = 0; mu < env().getNd(); ++mu) { @@ -44,12 +48,9 @@ std::vector TScalarVP::getOutput(void) // setup /////////////////////////////////////////////////////////////////////// void TScalarVP::setup(void) { - freeMomPropName_ = FREEMOMPROP(par().mass); - GFSrcName_ = "_" + getName() + "_DinvSrc"; - prop0Name_ = getName() + "_prop0"; - propQName_ = getName() + "_propQ"; - propSunName_ = getName() + "_propSun"; - propTadName_ = getName() + "_propTad"; + freeMomPropName_ = FREEMOMPROP(static_cast(env().getModule(par().scalarProp))->par().mass); + GFSrcName_ = "_" + par().scalarProp + "_DinvSrc"; + prop0Name_ = par().scalarProp + "_0"; phaseName_.clear(); muPropQName_.clear(); @@ -74,174 +75,38 @@ void TScalarVP::setup(void) freeVpTensorName_.push_back(freeVpTensorName_mu); } - if (!env().hasRegisteredObject(freeMomPropName_)) - { - env().registerLattice(freeMomPropName_); - } - if (!env().hasRegisteredObject(phaseName_[0])) - { - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { - env().registerLattice(phaseName_[mu]); - } - } - if (!env().hasRegisteredObject(GFSrcName_)) - { - env().registerLattice(GFSrcName_); - } - if (!env().hasRegisteredObject(prop0Name_)) - { - env().registerLattice(prop0Name_); - } - env().registerLattice(propQName_); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { env().registerLattice(muPropQName_[mu]); - } - env().registerLattice(propSunName_); - env().registerLattice(propTadName_); - for (unsigned int mu = 0; mu < env().getNd(); ++mu) - { + for (unsigned int nu = 0; nu < env().getNd(); ++nu) { env().registerLattice(vpTensorName_[mu][nu]); env().registerLattice(freeVpTensorName_[mu][nu]); } - } + } } // execution /////////////////////////////////////////////////////////////////// void TScalarVP::execute(void) { - // CACHING ANALYTIC EXPRESSIONS - ScalarField &source = *env().getObject(par().source); + // Get objects cached by ChargedProp module Complex ci(0.0,1.0); FFT fft(env().getGrid()); - Real q = par().charge; + Real q = static_cast(env().getModule(par().scalarProp))->par().charge; - // cache momentum-space free scalar propagator - if (!env().hasCreatedObject(freeMomPropName_)) + freeMomProp_ = env().getObject(freeMomPropName_); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) { - LOG(Message) << "Caching momentum space free scalar propagator" - << " (mass= " << par().mass << ")..." << std::endl; - freeMomProp_ = env().createLattice(freeMomPropName_); - SIMPL::MomentumSpacePropagator(*freeMomProp_, par().mass); - } - else - { - freeMomProp_ = env().getObject(freeMomPropName_); - } - // cache phases - if (!env().hasCreatedObject(phaseName_[0])) - { - std::vector &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(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(phaseName_[mu])); - } - } - // cache G*F*src - if (!env().hasCreatedObject(GFSrcName_)) - { - GFSrc_ = env().createLattice(GFSrcName_); - fft.FFT_all_dim(*GFSrc_, source, FFT::forward); - *GFSrc_ = (*freeMomProp_)*(*GFSrc_); - } - else - { - GFSrc_ = env().getObject(GFSrcName_); - } - // cache position-space free scalar propagators - if (!env().hasCreatedObject(prop0Name_)) - { - prop0_ = env().createLattice(prop0Name_); - fft.FFT_all_dim(*prop0_, *GFSrc_, FFT::backward); - } - else - { - prop0_ = env().getObject(prop0Name_); + phase_.push_back(env().getObject(phaseName_[mu])); } + GFSrc_ = env().getObject(GFSrcName_); + prop0_ = env().getObject(prop0Name_); - // PROPAGATOR CALCULATION // Propagator from unshifted source - LOG(Message) << "Computing O(alpha) charged scalar propagator" - << " (mass= " << par().mass - << ", charge= " << q << ")..." - << std::endl; - ScalarField &propQ = *env().createLattice(propQName_); - ScalarField &propSun = *env().createLattice(propSunName_); - ScalarField &propTad = *env().createLattice(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 vecBuf; - // std::vector 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); - // } + ScalarField &propQ = *env().getObject(propQName_); + ScalarField &propSun = *env().getObject(propSunName_); + ScalarField &propTad = *env().getObject(propTadName_); // Propagators from shifted sources LOG(Message) << "Computing O(q) charged scalar propagators..." @@ -357,7 +222,7 @@ void TScalarVP::execute(void) std::vector result; write(writer, "charge", q); - write(writer, "mass", par().mass); + write(writer, "mass", static_cast(env().getModule(par().scalarProp))->par().mass); 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) { EmField &A = *env().getObject(par().emField); @@ -443,32 +280,3 @@ void TScalarVP::momD1(ScalarField &s, FFT &fft) s = result; } - -void TScalarVP::momD2(ScalarField &s, FFT &fft) -{ - EmField &A = *env().getObject(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; -} diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp index fdbad6f6..81071ca0 100644 --- a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp +++ b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp @@ -17,9 +17,7 @@ class ScalarVPPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar, std::string, emField, - std::string, source, - double, mass, - double, charge, + std::string, scalarProp, std::string, output); }; @@ -42,11 +40,7 @@ public: // execution virtual void execute(void); private: - void chargedProp(ScalarField &prop_q, ScalarField &prop_sun, - ScalarField &prop_tad, ScalarField &GFSrc, - FFT &fft); void momD1(ScalarField &s, FFT &fft); - void momD2(ScalarField &s, FFT &fft); private: std::string freeMomPropName_, GFSrcName_, prop0Name_, propQName_,