diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 2aa0d65a..de8f6199 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -126,7 +126,8 @@ template void TDistilVectors::execute(void) { - auto &noise = envGet(std::vector>>, par().noise); + //auto &noise = envGet(std::vector>>, par().noise); + auto &noise = envGet(std::vector, par().noise); auto &perambulator = envGet(Perambulator, getName() + "_perambulator_light"); auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); auto &rho = envGet(std::vector, getName() + "_rho"); @@ -170,7 +171,8 @@ void TDistilVectors::execute(void) for (int ik = dk; ik < nvec; ik += LI){ for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced ExtractSliceLocal(evec3d,epack.evec[ik],0,it,3); - tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector + tmp3d_nospin = evec3d * noise[inoise + nnoise*(it + Nt*(ik+nvec*is))]; + //tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector tmp3d=zero; pokeSpin(tmp3d,tmp3d_nospin,is); tmp2=zero; diff --git a/Hadrons/Modules/MDistil/PerambLight.hpp b/Hadrons/Modules/MDistil/PerambLight.hpp index 774882d4..dbac54ca 100644 --- a/Hadrons/Modules/MDistil/PerambLight.hpp +++ b/Hadrons/Modules/MDistil/PerambLight.hpp @@ -24,7 +24,6 @@ class PerambLightPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(PerambLightPar, - std::string, noise, std::string, eigenPack, bool, multiFile, int, tsrc, @@ -77,7 +76,7 @@ std::vector TPerambLight::getInput(void) { std::vector in; - in.push_back(par().noise); + //in.push_back(par().noise); in.push_back(par().eigenPack); return in; @@ -86,7 +85,7 @@ std::vector TPerambLight::getInput(void) template std::vector TPerambLight::getOutput(void) { - std::vector out = {getName() + "_perambulator_light"}; + std::vector out = {getName() + "_perambulator_light",getName() + "_noise"}; return out; } @@ -96,13 +95,19 @@ template void TPerambLight::setup(void) { - auto &noise = envGet(std::vector>>, par().noise); + // auto &noise = envGet(std::vector>>, par().noise); - int nvec = 6; - int Nt=64; - - envCreate(Perambulator, getName() + "_perambulator_light", 1, - noise.size() *nvec*Nt); + int LI=par().LI; + int SI=par().SI; + int TI=par().TI; + int nnoise=par().nnoise; + int Nt=par().Nt; + int nvec=par().nvec; + + envCreate(Perambulator, getName() + "_perambulator_light", 1, + LI*SI*TI*nnoise*nvec*Nt); + envCreate(std::vector, getName() + "_noise", 1, + nvec*Ns*Nt*nnoise); GridCartesian * grid4d = env().getGrid(); std::vector latt_size = GridDefaultLatt(); @@ -133,7 +138,8 @@ template void TPerambLight::execute(void) { - auto &noise = envGet(std::vector>>, par().noise); + //auto &noise = envGet(std::vector>>, par().noise); + auto &noise = envGet(std::vector, getName() + "_noise"); auto &perambulator = envGet(Perambulator, getName() + "_perambulator_light"); auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); @@ -180,6 +186,30 @@ void TPerambLight::execute(void) int nvec=par().nvec; bool full_tdil=(TI==Nt); + bool exact_distillation = (full_tdil && LI==nvec); + + //Create Noises + //std::cout << pszGaugeConfigFile << std::endl; + //GridSerialRNG sRNG; sRNG.SeedUniqueString(std::string(pszGaugeConfigFile)); + GridSerialRNG sRNG; sRNG.SeedUniqueString("unique_string"); + Real rn; + + for (int inoise=0;inoise 0) - (rn-0.5 < 0); //TODO: This could be 0 if rn==0.5!! + //noises[inoise][t][ivec]()(is)() = (rn-0.5 > 0) - (rn-0.5 < 0); //TODO: This could be 0 if rn==0.5!! + } + } + } + } + } Real mass=par().mass; // TODO Infile Real M5 =par().M5; // TODO Infile @@ -218,7 +248,8 @@ void TPerambLight::execute(void) for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced std::cout << "LapH source vector from noise " << it << " and dilution component (d_k,d_t,d_alpha) : (" << ik << ","<< is << ")" << std::endl; ExtractSliceLocal(evec3d,epack.evec[ik],0,it,3); - tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector + tmp3d_nospin = evec3d * noise[inoise + nnoise*(it + Nt*(ik+nvec*is))]; + //tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector tmp3d=zero; pokeSpin(tmp3d,tmp3d_nospin,is); tmp2=zero; diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc index 3e411a58..17044563 100644 --- a/tests/hadrons/Test_hadrons_distil.cc +++ b/tests/hadrons/Test_hadrons_distil.cc @@ -274,24 +274,43 @@ void test_DistilVectors(Application &application) globalPar.trajCounter.step = 20; globalPar.runId = "test"; application.setPar(globalPar); - // Module parameters - MDistil::DistilVectors::Par DistilPar; - DistilPar.noise="noise"; - DistilPar.perambulator="perambulator"; - DistilPar.eigenPack="ePack"; - DistilPar.tsrc = 0; - DistilPar.nnoise = 1; - DistilPar.LI=6; - DistilPar.SI=4; - DistilPar.TI=64; - DistilPar.nvec=6; - DistilPar.Ns=4; - DistilPar.Nt=64; - DistilPar.Nt_inv=1; + // PerambLight parameters + MDistil::PerambLight::Par PerambPar; + PerambPar.eigenPack="ePack"; + PerambPar.tsrc = 0; + PerambPar.nnoise = 1; + PerambPar.LI=6; + PerambPar.SI=4; + PerambPar.TI=64; + PerambPar.nvec=6; + PerambPar.Ns=4; + PerambPar.Nt=64; + PerambPar.Nt_inv=1; + PerambPar.mass=0.005; + PerambPar.M5=1.8; + PerambPar.Ls=16; + PerambPar.CGPrecision=1e-8; + PerambPar.MaxIterations=10000; + // DistilVectors parameters + MDistil::DistilVectors::Par DistilVecPar; + DistilVecPar.noise="noise"; + DistilVecPar.perambulator="perambulator"; + DistilVecPar.eigenPack="ePack"; + DistilVecPar.tsrc = 0; + DistilVecPar.nnoise = 1; + DistilVecPar.LI=6; + DistilVecPar.SI=4; + DistilVecPar.TI=64; + DistilVecPar.nvec=6; + DistilVecPar.Ns=4; + DistilVecPar.Nt=64; + DistilVecPar.Nt_inv=1; // gauge field application.createModule("gauge"); - // Now make an instance of the LapEvec object - application.createModule("DistilVectorsInstance",DistilPar); + // Now make an instance of the Perambulator object + application.createModule("PerambulatorsInstance",PerambPar); + // Now make an instance of the DistilVectors object + application.createModule("DistilVectorsInstance",DistilVecPar); } bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )