From 4522f1e831a2710fe733398ab63624abea36422f Mon Sep 17 00:00:00 2001 From: ferben Date: Tue, 19 Feb 2019 17:22:30 +0000 Subject: [PATCH] separated final 2pt Contraction --- Hadrons/Modules/MDistil/BContraction.hpp | 23 ++- Hadrons/Modules/MDistil/Baryon2pt.cc | 7 + Hadrons/Modules/MDistil/Baryon2pt.hpp | 187 +++++++++++++++++++++++ Hadrons/Modules/MDistil/Distil.hpp | 8 + Hadrons/modules.inc | 2 + tests/hadrons/Test_hadrons_distil.cc | 17 +++ 6 files changed, 241 insertions(+), 3 deletions(-) create mode 100644 Hadrons/Modules/MDistil/Baryon2pt.cc create mode 100644 Hadrons/Modules/MDistil/Baryon2pt.hpp diff --git a/Hadrons/Modules/MDistil/BContraction.hpp b/Hadrons/Modules/MDistil/BContraction.hpp index 75c505cc..3e8c523f 100644 --- a/Hadrons/Modules/MDistil/BContraction.hpp +++ b/Hadrons/Modules/MDistil/BContraction.hpp @@ -72,6 +72,14 @@ protected: GridCartesian * grid3d; }; +/*class BFieldIO: Serializable{ +public: + using BaryonTensorSet = Eigen::Tensor; + GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, + BaryonTensorSet, BField + ); +};*/ + MODULE_REGISTER_TMP(BContraction, TBContraction, MDistil); /****************************************************************************** @@ -120,6 +128,7 @@ void TBContraction::execute(void) int N_3 = three.size(); int parity = par().parity; + const std::string &output{par().output}; LOG(Message) << "Computing distillation baryon fields" << std::endl; LOG(Message) << "One: '" << par().one << "' Two: '" << par().two << "' Three: '" << par().three << "'" << std::endl; @@ -133,7 +142,7 @@ void TBContraction::execute(void) grid3d = MakeLowerDimGrid(grid4d); int Nmom=1; int Nt=64; - std::vector BField(Nmom*Nt*N_1*N_2*N_3); + // std::vector BField(Nmom*Nt*N_1*N_2*N_3); int Bindex; int Nc=3; //Num colours @@ -223,7 +232,15 @@ void TBContraction::execute(void) } } - int Npairs = 0; + BFieldIO BField_save; + BField_save.BField = BField3; + + std::string filename ="./" + output + ".h5"; + std::cout << "Writing to file " << filename << std::endl; + Hdf5Writer writer(filename); + write(writer,"BaryonField",BField_save.BField); + + /* int Npairs = 0; char left[] = "uud"; char right[] = "uud"; std::vector pairs(6); @@ -261,7 +278,7 @@ void TBContraction::execute(void) for (int t=0 ; t < Nt ; t++){ std::cout << "C2(is=" << is << ",t=" << t << ") = " << corr(0,is,t) << std::endl; } - } + }*/ } END_MODULE_NAMESPACE diff --git a/Hadrons/Modules/MDistil/Baryon2pt.cc b/Hadrons/Modules/MDistil/Baryon2pt.cc new file mode 100644 index 00000000..58de7d33 --- /dev/null +++ b/Hadrons/Modules/MDistil/Baryon2pt.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MDistil; + +template class Grid::Hadrons::MDistil::TBaryon2pt; diff --git a/Hadrons/Modules/MDistil/Baryon2pt.hpp b/Hadrons/Modules/MDistil/Baryon2pt.hpp new file mode 100644 index 00000000..e97dab16 --- /dev/null +++ b/Hadrons/Modules/MDistil/Baryon2pt.hpp @@ -0,0 +1,187 @@ +#ifndef Hadrons_MDistil_Baryon2pt_hpp_ +#define Hadrons_MDistil_Baryon2pt_hpp_ + +#include +#include +#include +#include +#include +#include +#include + +// These are members of Distillation +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Baryon2pt * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MDistil) + +class Baryon2ptPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Baryon2ptPar, + std::string, inputL, + std::string, inputR, + std::string, output + ); +}; + +template +class TBaryon2pt: public Module +{ +public: + // constructor + TBaryon2pt(const std::string name); + // destructor + virtual ~TBaryon2pt(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +class C2IO: Serializable{ +public: +using C2Set = Eigen::Tensor; + GRID_SERIALIZABLE_CLASS_MEMBERS(C2IO, + C2Set, C2 + ); +}; + +MODULE_REGISTER_TMP(Baryon2pt, TBaryon2pt, MDistil); + +/****************************************************************************** + * TBaryon2pt implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBaryon2pt::TBaryon2pt(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBaryon2pt::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TBaryon2pt::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TBaryon2pt::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TBaryon2pt::execute(void) +{ + + const std::string &inputL{par().inputL}; + const std::string &inputR{par().inputR}; + const std::string &output{par().output}; + + int Nmom=1; + int Nt=64; + int Nc=3; //Num colours + std::vector> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + std::vector epsilon_sgn = {1,1,1,-1,-1,-1}; + int Ngamma=3; + + int N_1=20; + int N_2=20; + int N_3=20; + + // using BaryonTensorSet = Eigen::Tensor; + + BFieldIO BFieldL; + BFieldL.BField.resize(Nmom,4*Ngamma,Nt,N_1,N_2,N_3); + + std::string filenameL ="./" + inputL + ".h5"; + std::cout << "Reading from file " << filenameL << std::endl; + Hdf5Reader readerL(filenameL); + read(readerL,"BaryonField",BFieldL.BField); + + BFieldIO BFieldR; + BFieldR.BField.resize(Nmom,4*Ngamma,Nt,N_1,N_2,N_3); + + std::string filenameR ="./" + inputR + ".h5"; + std::cout << "Reading from file " << filenameR << std::endl; + Hdf5Reader readerR(filenameR); + read(readerR,"BaryonField",BFieldR.BField); + + Eigen::Tensor corr(Nmom,Nt); + + int Npairs = 0; + char left[] = "uud"; + char right[] = "uud"; + std::vector pairs(6); + for (int ie=0, i=0 ; ie < 6 ; ie++){ + if (left[0] == right[epsilon[ie][0]] && left[1] == right[epsilon[ie][1]] && left[2] == right[epsilon[ie][2]]){ + pairs[i] = ie; + i++; + Npairs++; + } + } + pairs.resize(Npairs); + std::cout << Npairs << " pairs: " << pairs << std::endl; + for (int imom=0 ; imom < Nmom ; imom++){ + for (int t=0 ; t < Nt ; t++){ + corr(imom,t) = 0.; + } + } + + int tsrc=0; + + for (int ipair=0 ; ipair < Npairs ; ipair++){ + Eigen::array, 3> product_dims = { Eigen::IndexPair(0,epsilon[pairs[ipair]][0]),Eigen::IndexPair(1,epsilon[pairs[ipair]][1]) ,Eigen::IndexPair(2,epsilon[pairs[ipair]][2]) }; + for (int imom=0 ; imom < Nmom ; imom++){ + std::cout << imom << std::endl; + Eigen::Tensor B5L = BFieldL.BField.chip(imom,0); + Eigen::Tensor B5R = BFieldR.BField.chip(imom,0); + for (int is=0 ; is < 4 ; is++){ + Eigen::Tensor B4L = B5L.chip(is,0); + Eigen::Tensor B4R = B5R.chip(is,0); + for (int t=0 ; t < Nt ; t++){ + Eigen::Tensor B3L = B4L.chip(t,0); + Eigen::Tensor B3R = B4R.chip(tsrc,0); + Eigen::Tensor C2 = B3L.contract(B3R,product_dims); + corr(imom,t) += (double)epsilon_sgn[pairs[ipair]]*C2(0); + } + } + } + } + for (int t=0 ; t < Nt ; t++){ + std::cout << "C2(t=" << t << ") = " << corr(0,t) << std::endl; + } + + C2IO C2_save; + C2_save.C2 = corr; + + std::string filename ="./" + output + ".h5"; + std::cout << "Writing to file " << filename << std::endl; + Hdf5Writer writer(filename); + write(writer,"C2",C2_save.C2); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MDistil_Baryon2pt_hpp_ diff --git a/Hadrons/Modules/MDistil/Distil.hpp b/Hadrons/Modules/MDistil/Distil.hpp index f2ada50f..d6673dd0 100644 --- a/Hadrons/Modules/MDistil/Distil.hpp +++ b/Hadrons/Modules/MDistil/Distil.hpp @@ -200,6 +200,14 @@ public: }; +class BFieldIO: Serializable{ +public: + using BaryonTensorSet = Eigen::Tensor; + GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, + BaryonTensorSet, BField + ); +}; + END_MODULE_NAMESPACE // Grid /****************************************************************************** diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 6d5ce763..441df51b 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -10,6 +10,7 @@ modules_cc =\ Modules/MSolver/RBPrecCG.cc \ Modules/MDistil/BContraction.cc \ Modules/MDistil/LapEvec.cc \ + Modules/MDistil/Baryon2pt.cc \ Modules/MDistil/PerambLight.cc \ Modules/MDistil/DistilVectors.cc \ Modules/MContraction/WeakHamiltonianEye.cc \ @@ -88,6 +89,7 @@ modules_hpp =\ Modules/MDistil/LapEvec.hpp \ Modules/MDistil/Distil.hpp \ Modules/MDistil/DistilVectors.hpp \ + Modules/MDistil/Baryon2pt.hpp \ Modules/MDistil/BContraction.hpp \ Modules/MDistil/PerambLight.hpp \ Modules/MContraction/WeakHamiltonian.hpp \ diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc index 8d865d20..bfa8c305 100644 --- a/tests/hadrons/Test_hadrons_distil.cc +++ b/tests/hadrons/Test_hadrons_distil.cc @@ -282,6 +282,19 @@ void test_BaryonFieldRho(Application &application) BContractionPar.mom={"0 0 0"}; application.createModule("BaryonFieldRho",BContractionPar); } +///////////////////////////////////////////////////////////// +// BaryonContraction +///////////////////////////////////////////////////////////// + +void test_Baryon2pt(Application &application) +{ + // DistilVectors parameters + MDistil::Baryon2pt::Par Baryon2ptPar; + Baryon2ptPar.inputL="BaryonFieldPhi"; + Baryon2ptPar.inputR="BaryonFieldRho"; + Baryon2ptPar.output="C2_baryon"; + application.createModule("C2_b",Baryon2ptPar); +} bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true ) { @@ -743,6 +756,10 @@ int main(int argc, char *argv[]) test_MesonField( application ); test_MesonFieldRho( application ); break; + case 9: // 3 + test_Global( application ); + test_Baryon2pt( application ); + break; } LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;