Added src&snk Gaussian smearing + APE smearing

projects/generate_confs/parser.ypp

@@ -237,6 +237,10 @@
 %token TK_GAUSS_NITER_SNK
 %type <double_numb> gauss_kappa
 %type <int_numb> gauss_niter_src gauss_niter_snk
+%token TK_APE_ALPHA
+%token TK_APE_NITER
+%type <double_numb> ape_alpha
+%type <int_numb> ape_niter
 //putpourri
 %type <fermionic_putpourri_meas> fermionic_putpourri_meas
 %token TK_MEAS_PUTPOURRI
@@ -609,6 +613,8 @@ meson_corr_meas: TK_MEAS_MESON_CORRS {$$=new meson_corr_meas_pars_t();}
                | meson_corr_meas gauss_kappa {$$->gauss_kappa=$2;}
                | meson_corr_meas gauss_niter_src {$$->gauss_niter_src=$2;}
                | meson_corr_meas gauss_niter_snk {$$->gauss_niter_snk=$2;}
+               | meson_corr_meas ape_alpha {$$->ape_alpha=$2;}
+               | meson_corr_meas ape_niter {$$->ape_niter=$2;}
 ;
 
 dir: TK_DIR '=' int_numb {$$=$3;}
@@ -623,6 +629,12 @@ gauss_niter_src: TK_GAUSS_NITER_SRC '=' int_numb {$$=$3;}
 gauss_niter_snk: TK_GAUSS_NITER_SNK '=' int_numb {$$=$3;}
 ;
 
+ape_alpha: TK_APE_ALPHA '=' double_numb {$$=$3;}
+;
+
+ape_niter: TK_APE_NITER '=' int_numb {$$=$3;}
+;
+
 operators: TK_OPERATORS '=' int_pair_list {
     for(size_t i=0;i<$3->size();i++)
 	{

projects/generate_confs/tokenizer.lpp

@@ -134,6 +134,10 @@ GaussKappa DEBUG_PRINTF("Found GaussKappa\n");return TK_GAUSS_KAPPA;
 GaussNiterSrc DEBUG_PRINTF("Found GaussNiterSrc\n");return TK_GAUSS_NITER_SRC;
 GaussNiterSnk DEBUG_PRINTF("Found GaussNiterSnk\n");return TK_GAUSS_NITER_SNK;
 
+ /* ape smearing (meson corrs) */
+ApeAlpha DEBUG_PRINTF("Found ApeAlpha\n");return TK_APE_ALPHA;
+ApeNiter DEBUG_PRINTF("Found ApeNiter\n");return TK_APE_NITER;
+
  /* rendens pars */
 MaxOrder DEBUG_PRINTF("Found MaxOrder\n");return TK_MAX_ORDER;
 

src/measures/fermions/mesons.cpp

@@ -16,6 +16,7 @@
 #include "routines/mpi_routines.hpp"
 #include "stag.hpp"
 #include "operations/smearing/gaussian.hpp"
+#include "operations/smearing/APE.hpp"
 
 namespace nissa
 {
@@ -26,9 +27,9 @@ namespace nissa
     int nop;
     int ncombo;
     int nflavs;
-    
-    /// Compute the index where to store
-    inline int icombo(const int& iflav_so,
+
+  /// Compute the index where to store
+  inline int icombo(const int& iflav_so,
 		      const int& iflav_si,
 		      const int& iop_so,
 		      const int& iop_si,
@@ -43,26 +44,30 @@ namespace nissa
     }
   }
 
-  // Reuse gaussian_smearing using lx as buffer
+
+
+  // reuse gaussian_smearing using lx as buffer
   template <typename T>
-  void gaussian_smear_eo(EoField<T>& out,
-			 const EoField<T>& in,
-			 LxField<T>& lx_buf,
-			 LxField<T>& lx_temp,
-			 LxField<T>& lx_H,
-			 const LxField<quad_su3>& conf_lx,
-			 const double gauss_kappa,
-			 const int gauss_niter)
-  {
-    if(gauss_niter<=0)
-      out=in;
-    else
-      {    
-	paste_eo_parts_into_lx_vector(lx_buf,in);
-	gaussian_smearing(lx_buf,lx_buf,conf_lx,gauss_kappa,gauss_niter,&lx_temp,&lx_H);
-	split_lx_vector_into_eo_parts(out,lx_buf);
-      }
-  }
+	void gaussian_smear_eo(EoField<T>& out,
+                       const EoField<T>& in,
+                       LxField<T>& lx_buf,
+                       LxField<T>& lx_temp,
+                       LxField<T>& lx_H,
+                       const LxField<quad_su3>& conf_lx,
+                       const double gauss_kappa,
+                       const int gauss_niter,
+					   const int corr_dir)
+    {
+      if(gauss_niter<=0)
+		out = in; 
+	  else
+		{
+	  paste_eo_parts_into_lx_vector(lx_buf,in);
+      gaussian_smearing(lx_buf,lx_buf,conf_lx,gauss_kappa,gauss_niter,corr_dir,&lx_temp,&lx_H);
+	  split_lx_vector_into_eo_parts(out,lx_buf);
+    	}
+	}
+
 
   void compute_meson_corr(complex* corr,
 			  const EoField<quad_su3>& conf,
@@ -71,9 +76,12 @@ namespace nissa
   {
     const int& dir=
       meas_pars.dir;
-    const double gauss_kappa=meas_pars.gauss_kappa;
-    const int gauss_niter_src=meas_pars.gauss_niter_src;
-    const int gauss_niter_snk=meas_pars.gauss_niter_snk;
+	const double gauss_kappa=meas_pars.gauss_kappa;
+	const int gauss_niter_src=meas_pars.gauss_niter_src;
+	const int gauss_niter_snk=meas_pars.gauss_niter_snk;
+	const double ape_alpha=meas_pars.ape_alpha;
+	const int ape_niter=meas_pars.ape_niter;
+
 
     //allocate
     EoField<color> ori_source("ori_source",WITH_HALO);
@@ -83,21 +91,21 @@ namespace nissa
     std::vector<EoField<color>> quark(nflavs*nop,{"quark",WITH_HALO});
     std::vector<EoField<color>> quark0s(nflavs*nop,{"quark0s",WITH_HALO});
 
-    LxField<quad_su3> conf_lx("conf_lx",WITH_HALO);
-    paste_eo_parts_into_lx_vector(conf_lx,conf);
-
+	LxField<quad_su3> conf_lx("conf_lx",WITH_HALO_EDGES);
+	LxField<quad_su3> conf_lx_ape("conf_lx_ape",WITH_HALO_EDGES);
+	paste_eo_parts_into_lx_vector(conf_lx,conf);
+	ape_orth_dir_smear_conf(conf_lx_ape, conf_lx, ape_alpha, ape_niter, dir);
+
     LxField<color> gauss_lx_buf("gauss_lx_buf",WITH_HALO);
     LxField<color> gauss_lx_temp("gauss_lx_temp",WITH_HALO);
     LxField<color> gauss_lx_H("gauss_lx_H",WITH_HALO);
-
-    auto gSmear=
-      [&](EoField<color>& out,
-	  const int& nIter,
-	  const EoField<color>& in)
+
+	auto gSmear =
+      [&](EoField<color>& out, const int& nIter, const EoField<color>& in)
       {
-	gaussian_smear_eo(out,in,gauss_lx_buf,gauss_lx_temp,gauss_lx_H,conf_lx,gauss_kappa,nIter);
+        gaussian_smear_eo(out, in, gauss_lx_buf, gauss_lx_temp, gauss_lx_H, conf_lx_ape, gauss_kappa, nIter, dir);
       };
-    
+
     //form the masks
     int mask[nop],shift[nop];
     for(int iop=0;iop<nop;iop++)
@@ -128,24 +136,23 @@ namespace nissa
 	generate_fully_undiluted_eo_source(ori_source,meas_pars.rnd_type,source_coord,dir);
 
 	EoField<color> smeared_ori_source("smeared_ori_source",WITH_HALO);
-	gSmear(smeared_ori_source,gauss_niter_src,ori_source);
-	
+    gSmear(smeared_ori_source, gauss_niter_src, ori_source);
+
 	for(int iflav=0;iflav<nflavs;iflav++)
 	  {
 	    for(int iop=0;iop<nop;iop++)
 	      {
-		const int idx=iflav*nop+iop;
-		apply_shift_op(source,temp[0],temp[1],conf,tp.backfield[iflav],shift[iop],smeared_ori_source);
-		put_stag_phases(source,mask[iop]);
-		
-		mult_Minv(quark[idx],conf,tp,iflav,meas_pars.residue,source);
-		
-		//smear also the sink prop before contraction
-		gSmear(quark[idx],gauss_niter_snk,quark[idx]);
+			const int idx=iflav*nop+iop;
+			apply_shift_op(source,temp[0],temp[1],conf,tp.backfield[iflav],shift[iop],smeared_ori_source);
+			put_stag_phases(source,mask[iop]);
+
+			mult_Minv(quark[idx],conf,tp,iflav,meas_pars.residue,source);
+
+			//smear also the sink prop before contraction
+			gSmear(quark[idx], gauss_niter_snk, quark[idx]);
 	      }
 	  }
-
-	/// Sink
+	    /// Sink
 	for(int iflav=0;iflav<nflavs;iflav++)
 	  for(int iop=0;iop<nop;iop++)
 	    {

src/measures/fermions/mesons.hpp

@@ -12,6 +12,10 @@ namespace nissa
     std::vector<std::pair<int,int>> mesons;
 
     int dir;
+
+    double ape_alpha;
+
+    int ape_niter;
 
     double gauss_kappa;
 
@@ -24,7 +28,18 @@ namespace nissa
     {
       return 0;
     }
-
+
+    //ape smearing default disabled
+    double def_ape_alpha() const 
+    {
+       return 0.5;
+    }
+
+    int def_ape_niter() const 
+    { 
+      return 0;
+    }
+
     // Smearing defaults: disabled
     double def_gauss_kappa() const
     {
@@ -72,15 +87,22 @@ namespace nissa
 	    }
 	  os<<"}\n";
 	}
+
+    if(ape_alpha!=def_ape_alpha() or full)
+  os<<" ApeAlpha\t=\t"<<ape_alpha<<"\n";
+
+    if(ape_niter!=def_ape_niter() or full)
+  os<<" ApeNiter\t=\t"<<ape_niter<<"\n";
+
+    if(gauss_kappa!=def_gauss_kappa() or full)
+  os<<" GaussKappa\t=\t"<<gauss_kappa<<"\n";
+
+    if(gauss_niter_src!=def_gauss_niter_src() or full)
+  os<<" GaussNiterSrc\t=\t"<<gauss_niter_src<<"\n";
+
+    if(gauss_niter_snk!=def_gauss_niter_snk() or full)
+  os<<" GaussNiterSnk\t=\t"<<gauss_niter_snk<<"\n";
 
-      if(gauss_kappa!=def_gauss_kappa() or full)
-	os<<" GaussKappa\t=\t"<<gauss_kappa<<"\n";
-
-      if(gauss_niter_src!=def_gauss_niter_src() or full)
-	os<<" GaussNiterSrc\t=\t"<<gauss_niter_src<<"\n";
-
-      if(gauss_niter_snk!=def_gauss_niter_snk() or full)
-	os<<" GaussNiterSnk\t=\t"<<gauss_niter_snk<<"\n";
 
       return os.str();
     }
@@ -91,15 +113,19 @@ namespace nissa
 	base_fermionic_meas_t::is_nonstandard() or
 	mesons.size() or
 	dir!=def_dir() or
-	gauss_kappa!=def_gauss_kappa() or
-	gauss_niter_src!=def_gauss_niter_src() or
-	gauss_niter_snk!=def_gauss_niter_snk() or
+  gauss_kappa!=def_gauss_kappa() or
+  gauss_niter_src!=def_gauss_niter_src() or
+  gauss_niter_snk!=def_gauss_niter_snk() or
+  ape_alpha!=def_ape_alpha() or
+  ape_niter!=def_ape_niter() or
 	path!=def_path();
     }
 
     meson_corr_meas_pars_t() :
       base_fermionic_meas_t(),
       dir(def_dir()),
+      ape_alpha(def_ape_alpha()),
+      ape_niter(def_ape_niter()),
       gauss_kappa(def_gauss_kappa()),
       gauss_niter_src(def_gauss_niter_src()),
       gauss_niter_snk(def_gauss_niter_snk())

src/operations/smearing/APE.cpp

@@ -11,7 +11,8 @@ namespace nissa
 		      const double& alpha,
 		      const int& nSteps,
 		      const WhichDirs& dirs,
-		      const int& minStapleDir)
+		      const int& minStapleDir,
+			  const int& skipStapleDir)
   {
     char listed_dirs[21]="";
     for(int mu=0;mu<NDIM;mu++)
@@ -32,7 +33,7 @@ namespace nissa
 	conf.updateEdges();
 
 	PAR(0,locVol,
-	    CAPTURE(alpha,dirs,minStapleDir,
+	    CAPTURE(alpha,dirs,minStapleDir,skipStapleDir,
 		    TO_READ(conf),
 		    TO_WRITE(smearConf)),
 	    ivol,
@@ -44,10 +45,10 @@ namespace nissa
 		    su3 stap,temp1,temp2;
 		    su3_put_to_zero(stap);
 
-		    for(int nu=minStapleDir;nu<NDIM;nu++)     //  E---F---C
-		      if(nu!=mu)                              //  |   |   | mu
-			{                                     //  D---A---B
-			  const int A=ivol;                   //   nu
+		    for(int nu=minStapleDir;nu<NDIM;nu++)     							//  E---F---C
+		      if(nu!=mu and (skipStapleDir<0 || nu!=skipStapleDir))             //  |   |   | mu
+			{                                     							//  D---A---B
+			  const int A=ivol;                  							//   nu
 			  const int B=loclxNeighup[A][nu];
 			  const int F=loclxNeighup[A][mu];
 			  unsafe_su3_prod_su3(temp1,conf[A][nu],conf[B][mu]);

src/operations/smearing/APE.hpp

@@ -67,7 +67,8 @@ namespace nissa
 		      const double& alpha,
 		      const int& nSteps,
 		      const WhichDirs& dirs=allDirs,
-		      const int& minStapleDir=0);
+		      const int& minStapleDir=0,
+          const int& skipStapleDir=-1);
 
   inline void ape_single_dir_smear_conf(LxField<quad_su3>& smear_conf,
 					const LxField<quad_su3>& origi_conf,
@@ -88,6 +89,17 @@ namespace nissa
   {
     ape_smear_conf(smear_conf,origi_conf,alpha,nstep,allOtherDirs[mu]);
   }
+
+  // - smear only mu != corr_dir
+  // - staples use only nu != corr_dir
+  inline void ape_orth_dir_smear_conf(LxField<quad_su3>& smear_conf,
+                                      const LxField<quad_su3>& origi_conf,
+                                      const double& alpha,
+                                      const int& nstep,
+                                      const int& corr_dir)
+  {
+    ape_smear_conf(smear_conf,origi_conf,alpha,nstep,allOtherDirs[corr_dir],0,corr_dir);
+  }
 
   inline void ape_temporal_smear_conf(LxField<quad_su3>& smear_conf,
 				      const LxField<quad_su3>& origi_conf,