From 143b75956c735ad349516b221822e19202fc0b41 Mon Sep 17 00:00:00 2001
From: Michael Marshall <43034299+mmphys@users.noreply.github.com>
Date: Fri, 19 Apr 2019 11:54:02 +0100
Subject: [PATCH] Stout smearing 3D fixes. Changed LapEvec to perform spatial
 smearing only

---
 Grid/qcd/smearing/StoutSmearing.h          | 28 ++++++++++++----------
 Hadrons/Modules/MDistil/LapEvec.hpp        |  8 +++----
 Hadrons/Modules/MGauge/StoutSmearing3D.hpp | 11 ++++-----
 3 files changed, 24 insertions(+), 23 deletions(-)

diff --git a/Grid/qcd/smearing/StoutSmearing.h b/Grid/qcd/smearing/StoutSmearing.h
index 65501c7e..27283d73 100644
--- a/Grid/qcd/smearing/StoutSmearing.h
+++ b/Grid/qcd/smearing/StoutSmearing.h
@@ -44,20 +44,20 @@ class Smear_Stout : public Smear<Gimpl> {
   // Smear<Gimpl>* ownership semantics:
   //    Smear<Gimpl>* passed in to constructor are owned by caller, so we don't delete them here
   //    Smear<Gimpl>* created within constructor need to be deleted as part of the destructor
-  const Smear<Gimpl>* SmearBase; // Not owned by this object, so not deleted at destruction
   const std::unique_ptr<Smear<Gimpl>> OwnedBase; // deleted at destruction
+  const Smear<Gimpl>* SmearBase; // Not owned by this object, so not deleted at destruction
 
-public:
-  INHERIT_GIMPL_TYPES(Gimpl)
-
-  // only anticipated to be used from default constructor, but might as well be visible
-  inline static std::vector<double> rho3D(double rho = 1.0, int orthogdim = -1){
+  // only anticipated to be used from default constructor
+  inline static std::vector<double> rho3D(double rho, int orthogdim){
     std::vector<double> rho3d(Nd*Nd);
     for (int mu=0; mu<Nd; mu++)
       for (int nu=0; nu<Nd; nu++)
-        rho3d[mu + Nd * nu] = (mu == orthogdim || nu == orthogdim) ? 0.0 : rho;
+        rho3d[mu + Nd * nu] = (mu == nu || mu == orthogdim || nu == orthogdim) ? 0.0 : rho;
     return rho3d;
   };
+  
+public:
+  INHERIT_GIMPL_TYPES(Gimpl)
 
   /*! Stout smearing with base explicitly specified */
   Smear_Stout(Smear<Gimpl>* base) : SmearBase{base} {
@@ -71,15 +71,17 @@ public:
     }
 
   /*! Default constructor. rho is constant in all directions, optionally except for orthogonal dimension */
-  /*Smear_Stout(double rho = 1.0, int orthogdim = -1)
-    : OwnedBase{new Smear_APE<Gimpl>(rho3D(rho,orthogdim))}, SmearBase{OwnedBase.get()} {
+  Smear_Stout(double rho, int orthogdim = -1)
+    : OwnedBase{(orthogdim<0 || orthogdim>=Nd) ? new Smear_APE<Gimpl>(rho) : new Smear_APE<Gimpl>(rho3D(rho,orthogdim))},
+      SmearBase{OwnedBase.get()} {
+    assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
+  }
+
+  /*
+  Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
     assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
   }*/
 
-  Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
-    assert(Nc == 3);//                  "Stout smearing currently implemented only for Nc==3");
-  }
-
   ~Smear_Stout() {}  // delete SmearBase...
 
   void smear(GaugeField& u_smr, const GaugeField& U) const {
diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp
index 8b9841ae..bccd873e 100644
--- a/Hadrons/Modules/MDistil/LapEvec.hpp
+++ b/Hadrons/Modules/MDistil/LapEvec.hpp
@@ -51,7 +51,7 @@ BEGIN_MODULE_NAMESPACE(MDistil)
 struct StoutParameters: Serializable {
   GRID_SERIALIZABLE_CLASS_MEMBERS(StoutParameters,
                                   int, steps,
-                                  double, parm)
+                                  double, parm) // TODO: change name of this to rho
   StoutParameters() = default;
   template <class ReaderClass> StoutParameters(Reader<ReaderClass>& Reader){read(Reader,"StoutSmearing",*this);}
 };
@@ -139,8 +139,8 @@ MODULE_REGISTER_TMP(LapEvec, TLapEvec<GIMPL>, MDistil);
 template <typename GImpl>
 TLapEvec<GImpl>::TLapEvec(const std::string name) : gridLD{nullptr}, Module<LapEvecPar>(name)
 {
-  LOG(Message) << "TLapEvec constructor : TLapEvec<GImpl>::TLapEvec(const std::string name)" << std::endl;
-  LOG(Message) << "this=" << this << ", gridLD=" << gridLD << std::endl;
+  //LOG(Message) << "TLapEvec constructor : TLapEvec<GImpl>::TLapEvec(const std::string name)" << std::endl;
+  //LOG(Message) << "this=" << this << ", gridLD=" << gridLD << std::endl;
 }
 
 // destructor /////////////////////////////////////////////////////////////////
@@ -235,7 +235,7 @@ void TLapEvec<GImpl>::execute(void)
   {
     const StoutParameters &Stout{par().Stout};
     envGetTmp(GaugeField, Umu_stout);
-    Smear_Stout<PeriodicGimplR> LS(Stout.parm);
+    Smear_Stout<PeriodicGimplR> LS(Stout.parm, Tdir);
     for (int i = 0; i < Stout.steps; i++) {
       LS.smear(Umu_stout, Umu_smear);
       Umu_smear = Umu_stout;
diff --git a/Hadrons/Modules/MGauge/StoutSmearing3D.hpp b/Hadrons/Modules/MGauge/StoutSmearing3D.hpp
index c5724943..3036c46a 100644
--- a/Hadrons/Modules/MGauge/StoutSmearing3D.hpp
+++ b/Hadrons/Modules/MGauge/StoutSmearing3D.hpp
@@ -108,13 +108,12 @@ void TStoutSmearing3D<GImpl>::setup(void)
 template <typename GImpl>
 void TStoutSmearing3D<GImpl>::execute(void)
 {
-    LOG(Message) << "Smearing '" << par().gauge << "' with " << par().steps
-                 << " step" << ((par().steps > 1) ? "s" : "") 
-                 << " of 3D-stout smearing and rho= " << par().rho
-		 << "orthogonal to dimension " << par().orthogdim << std::endl;
+    LOG(Message) << "Smearing '" << par().gauge
+                 << "' with " << par().steps << " step" << ((par().steps > 1) ? "s" : "")
+                 << " of 3D-stout smearing and rho=" << par().rho
+		 << " orthogonal to dimension " << par().orthogdim << std::endl;
 
-    //Smear_Stout<GImpl> smearer(par().rho, par().orthogdim);
-    Smear_Stout<GImpl> smearer(par().rho);
+    Smear_Stout<GImpl> smearer(par().rho, par().orthogdim);
     auto               &U    = envGet(GaugeField, par().gauge);
     auto               &Usmr = envGet(GaugeField, getName());