Hadrons: sequential gamma source

extras/Hadrons/Modules.hpp

@@ -7,5 +7,6 @@
 #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
 #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Grid/Hadrons/Modules/MSource/Point.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
 #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
 #include <Grid/Hadrons/Modules/Quark.hpp>

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -0,0 +1,137 @@
+#ifndef Hadrons_SeqGamma_hpp_
+#define Hadrons_SeqGamma_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Sequential source
+ -----------------------------
+ * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * gamma * exp(i x.mom)
+ 
+ * options:
+ - q: input propagator (string)
+ - tA: begin timeslice (integer)
+ - tB: end timesilce (integer)
+ - gamma: gamma product to insert (integer)
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ 
+ */
+
+/******************************************************************************
+ *                         SeqGamma                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class SeqGammaPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
+                                    std::string,  q,
+                                    unsigned int, tA,
+                                    unsigned int, tB,
+                                    unsigned int, gamma,
+                                    std::string,  mom);
+};
+
+template <typename FImpl>
+class TSeqGamma: public Module<SeqGammaPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TSeqGamma(const std::string name);
+    // destructor
+    virtual ~TSeqGamma(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+typedef TSeqGamma<FIMPL> SeqGamma;
+
+/******************************************************************************
+ *                 TSeqGamma implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSeqGamma<FImpl>::TSeqGamma(const std::string name)
+: Module<SeqGammaPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqGamma<FImpl>::execute(void)
+{
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source at t= " << par().tA << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating gamma_" << par().gamma
+                     << " sequential source for "
+                     << par().tA << " <= t <= " << par().tB << std::endl;
+    }
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    SpinMatrix                 g;
+    std::vector<Real>          p;
+    Complex                    i(0.0,1.0);
+    
+    g  = makeGammaProd(par().gamma);
+    p  = strToVec<Real>(par().mom);
+    ph = zero;
+    for(unsigned int mu = 0; mu < Nd; mu++)
+    {
+        LatticeCoordinate(coor, mu);
+        ph = ph + p[mu]*coor;
+    }
+    ph = exp(i*ph);
+    LatticeCoordinate(t, Tp);
+    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+}
+
+END_MODULE_NAMESPACE
+
+MODULE_REGISTER_NS(SeqGamma, MSource);
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SeqGamma_hpp_

extras/Hadrons/modules.inc

@@ -13,6 +13,7 @@ modules_hpp =\
   Modules/MGauge/Unit.hpp \
   Modules/MSolver/RBPrecCG.hpp \
   Modules/MSource/Point.hpp \
+  Modules/MSource/SeqGamma.hpp \
   Modules/MSource/Z2.hpp \
   Modules/Quark.hpp
 

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -0,0 +1,171 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_meson_3pt.cc
+ 
+ Copyright (C) 2015
+ 
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ 
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ 
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+    
+    // run setup ///////////////////////////////////////////////////////////////
+    Application                   application;
+
+    std::vector<std::string>              flavour = {"l", "s", "c1", "c2", "c3"};
+    std::vector<double>                   mass    = {.01, .04, .2  , .25 , .3  };
+    unsigned int                          nt      = GridDefaultLatt()[Tp];
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start = 1500;
+    globalPar.trajCounter.end   = 1520;
+    globalPar.trajCounter.step  = 20;
+    globalPar.seed              = "1 2 3 4";
+    globalPar.genetic.maxGen    = 5000;
+    globalPar.genetic.maxCstGen = 100;
+    globalPar.genetic.popSize   = 50;
+    globalPar.genetic.mutationRate = .25;
+    application.setPar(globalPar);
+    
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 12;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+    }
+    for (unsigned int t = 0; t < nt; t += 8)
+    {
+        std::string                           srcName;
+        std::vector<std::string>              qName;
+        std::vector<std::vector<std::string>> seqName;
+        
+        // Z2 source
+        MSource::Z2::Par z2Par;
+        z2Par.tA = t;
+        z2Par.tB = t;
+        srcName  = "z2_" + std::to_string(t);
+        application.createModule<MSource::Z2>(srcName, z2Par);
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        {
+            // sequential sources
+            MSource::SeqGamma::Par seqPar;
+            qName.push_back("QZ2_" + flavour[i] + srcName);
+            seqPar.q   = qName[i];
+            seqPar.tA  = (t + nt/4) % nt;
+            seqPar.tB  = (t + nt/4) % nt;
+            seqPar.mom = "1. 0. 0. 0.";
+            seqName.push_back(std::vector<std::string>(Nd));
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                seqPar.gamma   = 0x1 << mu;
+                seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
+                                 + "_" + std::to_string(seqPar.tA) + "-"
+                                 + qName[i];
+                application.createModule<MSource::SeqGamma>(seqName[i][mu], seqPar);
+            }
+            
+            // propagators
+            Quark::Par quarkPar;
+            quarkPar.solver = "CG_" + flavour[i];
+            quarkPar.source = srcName;
+            application.createModule<Quark>(qName[i], quarkPar);
+            for (unsigned int mu = 0; mu < Nd; ++mu)
+            {
+                quarkPar.source = seqName[i][mu];
+                seqName[i][mu]  = qName[i] + "-" + seqName[i][mu];
+                application.createModule<Quark>(seqName[i][mu], quarkPar);
+            }
+        }
+        
+        // contractions
+        MContraction::Meson::Par mesPar;
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = i; j < flavour.size(); ++j)
+        {
+            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.q1     = qName[i];
+            mesPar.q2     = qName[j];
+            application.createModule<MContraction::Meson>("meson_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j],
+                                                          mesPar);
+        }
+        for (unsigned int i = 0; i < flavour.size(); ++i)
+        for (unsigned int j = 0; j < flavour.size(); ++j)
+        for (unsigned int mu = 0; mu < Nd; ++mu)
+        {
+            MContraction::Meson::Par mesPar;
+            
+            mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_"
+                            + std::to_string(mu);
+            mesPar.q1     = qName[i];
+            mesPar.q2     = seqName[j][mu];
+            application.createModule<MContraction::Meson>("3pt_Z2_"
+                                                          + std::to_string(t)
+                                                          + "_"
+                                                          + flavour[i]
+                                                          + flavour[j]
+                                                          + "_"
+                                                          + std::to_string(mu),
+                                                          mesPar);
+        }
+    }
+    
+    Environment::getInstance().printContent();
+    // execution
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}