portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-13 01:05:36 +00:00
Code Releases Activity

weird bug...

Browse Source
This commit is contained in:
Felix Erben 2019-09-17 14:36:26 +01:00
parent f548114ff6
commit 76c93aa44e
2 changed files with 53 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
Grid/qcd/utils/BaryonUtils.h
View File

@ -53,8 +53,8 @@ public:
const PropagatorField &q3_src, const PropagatorField &q3_src,
const Gamma GammaA, const Gamma GammaA,
const Gamma GammaB, const Gamma GammaB,
const char quarks_snk[], const char * quarks_snk,
const char quarks_src[], const char * quarks_src,
const int parity, const int parity,
ComplexField &baryon_corr); ComplexField &baryon_corr);
@ -67,12 +67,15 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
const PropagatorField &q3_src, const PropagatorField &q3_src,
const Gamma GammaA, const Gamma GammaA,
const Gamma GammaB, const Gamma GammaB,
const char quarks_snk[], const char * quarks_snk,
const char quarks_src[], const char * quarks_src,
const int parity, const int parity,
ComplexField &baryon_corr) ComplexField &baryon_corr)
{ {
std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl;
std::cout << "GammaA " << (GammaA.g) << std::endl;
std::cout << "GammaB " << (GammaB.g) << std::endl;
assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
GridBase *grid = q1_src.Grid(); GridBase *grid = q1_src.Grid();
@ -88,12 +91,16 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
wick_contraction[ie]=1; wick_contraction[ie]=1;
typedef typename ComplexField::vector_object vobj; typedef typename ComplexField::vector_object vobj;
LatticeView<vobj> vbaryon_corr{ baryon_corr }; auto vbaryon_corr= baryon_corr.View();
accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto v1 = q1_src.View();
auto v2 = q2_src.View();
auto v3 = q3_src.View();
std::cout << "wick contract " << wick_contraction << std::endl;
int count=0;
// accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
thread_for(ss,grid->oSites(),{
LatticeView<pobj> v1(q1_src);
LatticeView<pobj> v2(q2_src);
LatticeView<pobj> v3(q3_src);
auto D1 = v1[ss]; auto D1 = v1[ss];
auto D2 = v2[ss]; auto D2 = v2[ss];
auto D3 = v3[ss]; auto D3 = v3[ss];
@ -103,7 +110,11 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
auto pD1 = 0.5* (gD1a + (double)parity * gD1b); auto pD1 = 0.5* (gD1a + (double)parity * gD1b);
auto gD3 = GammaB * D3; auto gD3 = GammaB * D3;
vobj result{ 0 }; vobj result=Zero();
if (count<10){
std::cout << "outside epsilon " << count << std::endl;
count++;
}
for (int ie_src=0; ie_src < 6 ; ie_src++){ for (int ie_src=0; ie_src < 6 ; ie_src++){
int a_src = epsilon[ie_src][0]; //a int a_src = epsilon[ie_src][0]; //a
@ -132,6 +143,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
}}} }}}
} }
//This is the \delta_{456}^{312} part //This is the \delta_{456}^{312} part
// for(int test=0;test<3;test++){
if (wick_contraction[2]){ if (wick_contraction[2]){
auto gD3g = gD3 * GammaB; auto gD3g = gD3 * GammaB;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){ for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
@ -149,6 +161,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,gamma_src)(c_snk,c_src)*D2()(alpha_snk,beta_src)(a_snk,b_src)*gD3g()(alpha_snk,beta_src)(b_snk,a_src); result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,gamma_src)(c_snk,c_src)*D2()(alpha_snk,beta_src)(a_snk,b_src)*gD3g()(alpha_snk,beta_src)(b_snk,a_src);
}}} }}}
} }
// }
//This is the \delta_{456}^{321} part //This is the \delta_{456}^{321} part
if (wick_contraction[4]){ if (wick_contraction[4]){
auto D2g = D2 * GammaB; auto D2g = D2 * GammaB;

32
Hadrons/Modules/MContraction/Baryon.hpp
View File

@ -34,6 +34,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Hadrons/Module.hpp> #include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp> #include <Hadrons/ModuleFactory.hpp>
#include <Grid/qcd/utils/BaryonUtils.h> #include <Grid/qcd/utils/BaryonUtils.h>
#include <Grid/qcd/utils/A2Autils.h>
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
@ -134,24 +135,49 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
<< " ) and parity " << par().parity << "." << std::endl; << " ) and parity " << par().parity << "." << std::endl;
auto &q1_src = envGet(PropagatorField1, par().q1_src); auto &q1_src = envGet(PropagatorField1, par().q1_src);
LOG(Message) << "1" << std::endl;
auto &q2_src = envGet(PropagatorField2, par().q2_src); auto &q2_src = envGet(PropagatorField2, par().q2_src);
auto &q3_src = envGet(PropagatorField3, par().q3_src); auto &q3_src = envGet(PropagatorField3, par().q3_src);
envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, c);
Result result; Result result;
int nt = env().getDim(Tp); int nt = env().getDim(Tp);
result.corr.resize(nt); result.corr.resize(nt);
LOG(Message) << "2" << std::endl;
std::vector<Gamma::Algebra> ggA = strToVec<Gamma::Algebra>(par().GammaA); std::vector<Gamma::Algebra> ggA = strToVec<Gamma::Algebra>(par().GammaA);
Gamma GammaA = ggA[0]; LOG(Message) << "3" << std::endl;
Gamma GammaA(ggA[0]);
LOG(Message) << "4" << std::endl;
std::vector<Gamma::Algebra> ggB = strToVec<Gamma::Algebra>(par().GammaB); std::vector<Gamma::Algebra> ggB = strToVec<Gamma::Algebra>(par().GammaB);
Gamma GammaB = ggB[0]; Gamma GammaB(ggB[0]);
std::vector<TComplex> buf; std::vector<TComplex> buf;
const int parity {par().parity}; const int parity {par().parity};
LOG(Message) << "5" << std::endl;
const char * quarks_snk{par().quarks_snk.c_str()}; const char * quarks_snk{par().quarks_snk.c_str()};
LOG(Message) << "6" << std::endl;
const char * quarks_src{par().quarks_src.c_str()}; const char * quarks_src{par().quarks_src.c_str()};
LOG(Message) << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl;
LOG(Message) << "GammaA " << (GammaA.g) << std::endl;
LOG(Message) << "GammaB " << (GammaB.g) << std::endl;
BaryonUtils<FIMPL>::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); GridCartesian *Ugrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
LatticePropagator q(Ugrid);
GridParallelRNG RNG4(Ugrid);
gaussian(RNG4,q);
Gamma gA = Gamma(Gamma::Algebra::Identity);
Gamma gB = Gamma(Gamma::Algebra::Gamma5);
int p=1;
char * om = "sss";
LatticeComplex c2(Ugrid);
//BaryonUtils<FIMPL>::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c);
BaryonUtils<FIMPL>::ContractBaryons(q,q,q,gA,gB,om,om,p,c);
std::vector<Gamma> GA={gA};
std::vector<Gamma> GB={gB};
//A2Autils<FIMPL>::ContractFourQuarkColourMix(q,q,GA,GB,c,c2);
LOG(Message) << "survived ContractBaryons" << std::endl;
sliceSum(c,buf,Tp); sliceSum(c,buf,Tp);
LOG(Message) << "survived sliceSum" << std::endl;
for (unsigned int t = 0; t < buf.size(); ++t) for (unsigned int t = 0; t < buf.size(); ++t)
{ {