diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 075d6f41..a250d869 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 970e8679..2efd10fa 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -85,7 +85,7 @@ std::vector TDistilVectors::getInput(void) template std::vector TDistilVectors::getOutput(void) { - std::vector out = {getName() + "_rho", getName() + "_phi", getName() + "_rho_all_tsrc"}; + std::vector out = {getName() + "_rho", getName() + "_phi"}; return out; } @@ -106,9 +106,6 @@ void TDistilVectors::setup(void) nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); envCreate(std::vector, getName() + "_phi", 1, nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); - envCreate(std::vector, getName() + "_rho_all_tsrc", 1, - nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField)); - GridCartesian * grid4d = env().getGrid(); std::vector latt_size = GridDefaultLatt(); @@ -140,8 +137,6 @@ void TDistilVectors::execute(void) auto &epack = envGet(Grid::Hadrons::EigenPack, par().eigenPack); auto &rho = envGet(std::vector, getName() + "_rho"); auto &phi = envGet(std::vector, getName() + "_phi"); - auto &rho_all = envGet(std::vector, getName() + "_rho_all_tsrc"); - envGetTmp(LatticeSpinColourVector, tmp2); envGetTmp(LatticeColourVector, tmp_nospin); @@ -197,34 +192,6 @@ void TDistilVectors::execute(void) } } - for (int inoise = 0; inoise < nnoise; inoise++) { - for (int dk = 0; dk < LI; dk++) { - for (int dt = 0; dt < Nt_inv; dt++) { - for (int ds = 0; ds < Ns; ds++) { - vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * Ns*dt; - rho_all[vecindex] = zero; - tmp3d_nospin = zero; - for (int it = 0; it < Nt; it++){ - t_inv = it; - if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) { - 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,t_inv,3); - tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))]; - tmp3d=zero; - pokeSpin(tmp3d,tmp3d_nospin,is); - tmp2=zero; - InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir); - rho_all[vecindex] += tmp2; - } - } - } - } - } - } - } - } - for (int inoise = 0; inoise < nnoise; inoise++) { for (int dk = 0; dk < LI; dk++) { for (int dt = 0; dt < Nt_inv; dt++) { diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index d994d41e..eb4c6d88 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -47,6 +47,7 @@ modules_cc =\ Modules/MDistil/Baryon2pt.cc \ Modules/MDistil/BC2.cc \ Modules/MDistil/PerambLight.cc \ + Modules/MDistil/DistilSink.cc \ Modules/MSource/Momentum.cc \ Modules/MSource/Z2.cc \ Modules/MSource/Point.cc \ @@ -130,6 +131,7 @@ modules_hpp =\ Modules/MDistil/LapEvec.hpp \ Modules/MDistil/BContraction.hpp \ Modules/MDistil/DistilVectors.hpp \ + Modules/MDistil/DistilSink.hpp \ Modules/MSource/SeqConserved.hpp \ Modules/MSource/Z2.hpp \ Modules/MSource/Wall.hpp \ diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc index 5852b854..cb97d727 100644 --- a/tests/hadrons/Test_hadrons_distil.cc +++ b/tests/hadrons/Test_hadrons_distil.cc @@ -178,7 +178,7 @@ void test_PerambulatorsS(Application &application) PerambPar.Distil.Ns=4; PerambPar.Distil.Nt=8; PerambPar.Distil.Nt_inv=1; - PerambPar.Solver.mass=0.04; //strange mass??? + PerambPar.Solver.mass=0.005; //strange mass??? PerambPar.Solver.M5=1.8; PerambPar.Ls=16; PerambPar.Solver.CGPrecision=1e-8; @@ -200,10 +200,10 @@ void test_DistilVectorsS(Application &application) DistilVecPar.nnoise = 1; DistilVecPar.LI=3; DistilVecPar.SI=4; - DistilVecPar.TI=8; + DistilVecPar.TI=32; DistilVecPar.nvec=3; DistilVecPar.Ns=4; - DistilVecPar.Nt=8; + DistilVecPar.Nt=32; DistilVecPar.Nt_inv=1; application.createModule("DistilVecsS",DistilVecPar); } @@ -427,12 +427,69 @@ void test_AslashSeq(Application &application) { // DistilVectors parameters MSolver::A2AAslashVectors::Par A2AAslashVectorsPar; - A2AAslashVectorsPar.vector="Peramb_unsmeared_sink"; + A2AAslashVectorsPar.vector="PerambS_unsmeared_sink"; A2AAslashVectorsPar.emField="Em"; A2AAslashVectorsPar.solver="CG_s"; A2AAslashVectorsPar.output="Aslash_seq"; application.createModule("Aslash_seq",A2AAslashVectorsPar); } +///////////////////////////////////////////////////////////// +// DistilVectors +///////////////////////////////////////////////////////////// + +void test_DistilVectorsAslashSeq(Application &application) +{ + // DistilVectors parameters + MDistil::DistilSink::Par DistilSinkPar; + DistilSinkPar.perambulator="PerambAslashS"; + DistilSinkPar.eigenPack="LapEvec"; + DistilSinkPar.tsrc = 0; + DistilSinkPar.nnoise = 1; + DistilSinkPar.LI=3; + DistilSinkPar.SI=4; + DistilSinkPar.TI=8; + DistilSinkPar.nvec=3; + DistilSinkPar.Ns=4; + DistilSinkPar.Nt=8; + DistilSinkPar.Nt_inv=1; + application.createModule("DistilVecsAslashSeq",DistilSinkPar); +} +void test_PerambulatorsSolve(Application &application) +{ + // PerambLight parameters + MDistil::PerambFromSolve::Par PerambFromSolvePar; + PerambFromSolvePar.eigenPack="LapEvec"; + PerambFromSolvePar.solve="Aslash_seq"; + PerambFromSolvePar.PerambFileName="perambAslashS.bin"; + PerambFromSolvePar.Distil.tsrc = 0; + PerambFromSolvePar.Distil.nnoise = 1; + PerambFromSolvePar.Distil.LI=3; + PerambFromSolvePar.Distil.SI=4; + PerambFromSolvePar.Distil.TI=8; + PerambFromSolvePar.nvec=3; + PerambFromSolvePar.Distil.Ns=4; + PerambFromSolvePar.Distil.Nt=8; + PerambFromSolvePar.Distil.Nt_inv=1; + application.createModule("PerambAslashS",PerambFromSolvePar); +} +///////////////////////////////////////////////////////////// +// MesonFields - aslaaaash +///////////////////////////////////////////////////////////// + +void test_MesonFieldAslashSeq(Application &application) +{ + // DistilVectors parameters + MContraction::A2AMesonField::Par A2AMesonFieldPar; + A2AMesonFieldPar.left="DistilVecsAslashSeq"; + //A2AMesonFieldPar.right="DistilVecs_rho"; + A2AMesonFieldPar.right="DistilVecsAslashSeq"; + A2AMesonFieldPar.output="MesonSinksAslashSeq"; + A2AMesonFieldPar.gammas="all"; + A2AMesonFieldPar.mom={"0 0 0"}; + A2AMesonFieldPar.cacheBlock=2; + A2AMesonFieldPar.block=4; + application.createModule("DistilMesonFieldAslashSeq",A2AMesonFieldPar); +} bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true ) { @@ -917,6 +974,17 @@ int main(int argc, char *argv[]) test_BaryonFieldPhi2( application ); test_BaryonFieldRho2( application ); break; + case 12: // 3 + test_Global( application ); + test_SolverS( application ); + test_LapEvec( application ); + test_PerambulatorsS( application ); + test_em( application ); + test_AslashSeq( application ); + test_PerambulatorsSolve( application ); + test_DistilVectorsAslashSeq( application ); + test_MesonFieldAslashSeq( application ); + break; } LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;