fixing value of alpha

Reviewed-on: portelli/HadronsPresets#1

RbcUkqcd.hpp

@@ -44,10 +44,10 @@ struct RbcUkqcd
   inline static constexpr EnsembleParameters m0LCDPar{0.0006203, 0.02661, 1.8, 2., 64, 128, 12};
 
   // C1M
-  inline static constexpr EnsembleParameters c1mLCDPar  {0.005, 0.0362, 1.8, 2., 24, 64, 24};
-  inline static constexpr EnsembleParameters c1m16LCDPar{0.005, 0.0362, 1.8, 2., 16, 64, 24};
-  inline static constexpr EnsembleParameters c1m20LCDPar{0.005, 0.0362, 1.8, 2., 20, 64, 24};
-  inline static constexpr EnsembleParameters c1m32LCDPar{0.005, 0.0362, 1.8, 2., 32, 64, 24};
+  inline static constexpr EnsembleParameters c1mIRLPar{0.005, 0.0362, 1.8, 2., 24, 64, 24};
+  inline static constexpr EnsembleParameters c1m16IRLPar{0.005, 0.0362, 1.8, 2., 16, 64, 24};
+  inline static constexpr EnsembleParameters c1m20IRLPar{0.005, 0.0362, 1.8, 2., 20, 64, 24};
+  inline static constexpr EnsembleParameters c1m32IRLPar{0.005, 0.0362, 1.8, 2., 32, 64, 24};
 
   // Runtime deflation parameters
   struct DeflationParameters
@@ -56,80 +56,93 @@ struct RbcUkqcd
     unsigned int nPoly, nStop, nK, nM;
   };
 
-  inline static constexpr DeflationParameters c1mDeflPar  {5.0e-04, 5.5, 101, 100, 110, 120};
+  inline static constexpr DeflationParameters c1mDeflPar{5.0e-04, 5.5, 101, 100, 110, 120};
   inline static constexpr DeflationParameters c1m16DeflPar{5.0e-03, 5.5, 101, 200, 220, 230};
   inline static constexpr DeflationParameters c1m20DeflPar{3.0e-03, 5.5, 101, 200, 220, 230};
-  inline static constexpr DeflationParameters c1m32DeflPar{3.0e-06, 5.5, 101, 100, 110, 120};
+  inline static constexpr DeflationParameters c1m32DeflPar{5.0e-05, 5.5, 101, 100, 110, 120};
 
-  // Light solver: load deflation from disk
-  static inline void addLightLCDSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
-                                       const std::string solverName, const std::string gaugeName,
-                                       const std::string gaugeTransform,
-                                       const std::string eigenpackPath, const double residual);
-
-  static inline void addM0LightLCDSolver(Application &app, const std::string solverName,
-                                         const std::string gaugeName,
-                                         const std::string gaugeTransform,
-                                         const std::string eigenpackPath,
-                                         const double residual = 1.0e-8);
-  // Light solver: deflation at runtime
+  // Light solvers: load deflation from disk
   static inline void
-  addLightRuntimeLCDSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
-                           const RbcUkqcd::DeflationParameters &deflPar,
-                           const std::string solverName, const std::string gaugeName,
-                           const std::string gaugeTransform, const double residual);
+  addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName,
+                      const std::string gaugeTransform, const std::string eigenpackPath,
+                      const std::string boundary = "1 1 1 1", const double residual = 1.0e-8);
 
-  static inline void addC1MLightRuntimeLCDSolver(Application &app, const std::string solverName,
+  // Light solvers: deflation at runtime
+  static inline void addLightRuntimeIRLSolver(Application &app,
+                                              const RbcUkqcd::EnsembleParameters &par,
+                                              const RbcUkqcd::DeflationParameters &deflPar,
+                                              const std::string solverName,
+                                              const std::string gaugeName,
+                                              const std::string gaugeTransform,
+                                              const std::string boundary, const double residual);
+
+  static inline void addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName,
                                                  const std::string gaugeName,
                                                  const std::string gaugeTransform,
+                                                 const std::string boundary = "1 1 1 -1",
                                                  const double residual = 1.0e-8);
-  static inline void addC1M16LightRuntimeLCDSolver(Application &app, const std::string solverName,
+  static inline void addC1M16LightRuntimeIRLSolver(Application &app, const std::string solverName,
                                                    const std::string gaugeName,
                                                    const std::string gaugeTransform,
+                                                   const std::string boundary = "1 1 1 -1",
                                                    const double residual = 1.0e-8);
-  static inline void addC1M20LightRuntimeLCDSolver(Application &app, const std::string solverName,
+  static inline void addC1M20LightRuntimeIRLSolver(Application &app, const std::string solverName,
                                                    const std::string gaugeName,
                                                    const std::string gaugeTransform,
+                                                   const std::string boundary = "1 1 1 -1",
                                                    const double residual = 1.0e-8);
-  static inline void addC1M32LightRuntimeLCDSolver(Application &app, const std::string solverName,
+  static inline void addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName,
                                                    const std::string gaugeName,
                                                    const std::string gaugeTransform,
+                                                   const std::string boundary = "1 1 1 -1",
                                                    const double residual = 1.0e-8);
 
-  // Strange solver (undeflated)
+  // Strange solvers (undeflated)
   static inline void addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
                                       const std::string solverName, const std::string gaugeName,
-                                      const std::string gaugeTransform, const double residual);
+                                      const std::string gaugeTransform, const std::string boundary,
+                                      const double residual);
 
   static inline void addM0StrangeSolver(Application &app, const std::string solverName,
                                         const std::string gaugeName,
                                         const std::string gaugeTransform,
+                                        const std::string boundary = "1 1 1 1",
                                         const double residual = 1.0e-8);
   static inline void addC1MStrangeSolver(Application &app, const std::string solverName,
                                          const std::string gaugeName,
                                          const std::string gaugeTransform,
+                                         const std::string boundary = "1 1 1 -1",
                                          const double residual = 1.0e-8);
   static inline void addC1M16StrangeSolver(Application &app, const std::string solverName,
                                            const std::string gaugeName,
                                            const std::string gaugeTransform,
+                                           const std::string boundary = "1 1 1 -1",
                                            const double residual = 1.0e-8);
   static inline void addC1M20StrangeSolver(Application &app, const std::string solverName,
                                            const std::string gaugeName,
                                            const std::string gaugeTransform,
+                                           const std::string boundary = "1 1 1 -1",
                                            const double residual = 1.0e-8);
   static inline void addC1M32StrangeSolver(Application &app, const std::string solverName,
                                            const std::string gaugeName,
                                            const std::string gaugeTransform,
+                                           const std::string boundary = "1 1 1 -1",
                                            const double residual = 1.0e-8);
+
+  // Charm solvers, mass is a free parameter
+  static inline void addM0CharmSolver(Application &app, const std::string solverName,
+                                      const std::string gaugeName, const std::string gaugeTransform,
+                                      const double mass, const std::string boundary = "1 1 1 1",
+                                      const double residual = 1.0e-18);
 };
 
 // Implementations /////////////////////////////////////////////////////////////////////////////////
 
-// Load deflation from disk
-void RbcUkqcd::addLightLCDSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
-                                 const std::string solverName, const std::string gaugeName,
-                                 const std::string gaugeTransform, const std::string eigenpackPath,
-                                 const double residual)
+// Light M0 (load deflation from disk)
+void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName,
+                                   const std::string gaugeName, const std::string gaugeTransform,
+                                   const std::string eigenpackPath, const std::string boundary,
+                                   const double residual)
 {
   const std::string prefix = solverName;
 
@@ -143,11 +156,11 @@ void RbcUkqcd::addLightLCDSolver(Application &app, const RbcUkqcd::EnsembleParam
   MAction::ScaledDWF::Par actionPar;
 
   actionPar.gauge = gaugeName;
-  actionPar.Ls = par.Ls;
-  actionPar.M5 = par.M5;
-  actionPar.mass = par.ml;
-  actionPar.scale = par.scale;
-  actionPar.boundary = "1 1 1 1";
+  actionPar.Ls = RbcUkqcd::m0LCDPar.Ls;
+  actionPar.M5 = RbcUkqcd::m0LCDPar.M5;
+  actionPar.mass = RbcUkqcd::m0LCDPar.ml;
+  actionPar.scale = RbcUkqcd::m0LCDPar.scale;
+  actionPar.boundary = boundary;
   actionPar.twist = "0. 0. 0. 0.";
   app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
   actionPar.gauge = prefix + "_gauge_fp32";
@@ -189,20 +202,12 @@ void RbcUkqcd::addLightLCDSolver(Application &app, const RbcUkqcd::EnsembleParam
   app.createModule<MSolver::MixedPrecisionRBPrecCGBatched>(solverName, solverPar);
 }
 
-// Light M0
-void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName,
-                                   const std::string gaugeName, const std::string gaugeTransform,
-                                   const std::string eigenpackPath, const double residual)
-{
-  RbcUkqcd::addLightLCDSolver(app, RbcUkqcd::m0LCDPar, solverName, gaugeName, gaugeTransform,
-                              eigenpackPath, residual);
-}
-
 // Deflation at runtime
-void RbcUkqcd::addLightRuntimeLCDSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
+void RbcUkqcd::addLightRuntimeIRLSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
                                         const RbcUkqcd::DeflationParameters &deflPar,
                                         const std::string solverName, const std::string gaugeName,
-                                        const std::string gaugeTransform, const double residual)
+                                        const std::string gaugeTransform,
+                                        const std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -220,7 +225,7 @@ void RbcUkqcd::addLightRuntimeLCDSolver(Application &app, const RbcUkqcd::Ensemb
   actionPar.M5 = par.M5;
   actionPar.mass = par.ml;
   actionPar.scale = par.scale;
-  actionPar.boundary = "1 1 1 1";
+  actionPar.boundary = boundary;
   actionPar.twist = "0. 0. 0. 0.";
   app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
   actionPar.gauge = prefix + "_gauge_fp32";
@@ -267,45 +272,50 @@ void RbcUkqcd::addLightRuntimeLCDSolver(Application &app, const RbcUkqcd::Ensemb
 }
 
 // Light C1M
-void RbcUkqcd::addC1MLightRuntimeLCDSolver(Application &app, const std::string solverName,
+void RbcUkqcd::addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName,
                                            const std::string gaugeName,
-                                           const std::string gaugeTransform, const double residual)
+                                           const std::string gaugeTransform,
+                                           const std::string boundary, const double residual)
 {
-  RbcUkqcd::addLightRuntimeLCDSolver(app, RbcUkqcd::c1mLCDPar, RbcUkqcd::c1mDeflPar, solverName,
-                                     gaugeName, gaugeTransform, residual);
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1mIRLPar, RbcUkqcd::c1mDeflPar, solverName,
+                                     gaugeName, gaugeTransform, boundary, residual);
 }
 
 // Light C1M16
-void RbcUkqcd::addC1M16LightRuntimeLCDSolver(Application &app, const std::string solverName,
-                                           const std::string gaugeName,
-                                           const std::string gaugeTransform, const double residual)
+void RbcUkqcd::addC1M16LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string gaugeTransform,
+                                             const std::string boundary, const double residual)
 {
-  RbcUkqcd::addLightRuntimeLCDSolver(app, RbcUkqcd::c1m16LCDPar, RbcUkqcd::c1m16DeflPar, solverName,
-                                     gaugeName, gaugeTransform, residual);
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m16IRLPar, RbcUkqcd::c1m16DeflPar, solverName,
+                                     gaugeName, gaugeTransform, boundary, residual);
 }
 
 // Light C1M20
-void RbcUkqcd::addC1M20LightRuntimeLCDSolver(Application &app, const std::string solverName,
-                                           const std::string gaugeName,
-                                           const std::string gaugeTransform, const double residual)
+void RbcUkqcd::addC1M20LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string gaugeTransform,
+                                             const std::string boundary, const double residual)
 {
-  RbcUkqcd::addLightRuntimeLCDSolver(app, RbcUkqcd::c1m20LCDPar, RbcUkqcd::c1m20DeflPar, solverName,
-                                     gaugeName, gaugeTransform, residual);
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m20IRLPar, RbcUkqcd::c1m20DeflPar, solverName,
+                                     gaugeName, gaugeTransform, boundary, residual);
 }
 
 // Light C1M32
-void RbcUkqcd::addC1M32LightRuntimeLCDSolver(Application &app, const std::string solverName,
-                                           const std::string gaugeName,
-                                           const std::string gaugeTransform, const double residual)
+void RbcUkqcd::addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string gaugeTransform,
+                                             const std::string boundary, const double residual)
 {
-  RbcUkqcd::addLightRuntimeLCDSolver(app, RbcUkqcd::c1m32LCDPar, RbcUkqcd::c1m32DeflPar, solverName,
-                                     gaugeName, gaugeTransform, residual);
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m32IRLPar, RbcUkqcd::c1m32DeflPar, solverName,
+                                     gaugeName, gaugeTransform, boundary, residual);
 }
 
-// Strange 
+// Strange
 void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
                                 const std::string solverName, const std::string gaugeName,
-                                const std::string gaugeTransform, const double residual)
+                                const std::string gaugeTransform, const std::string boundary,
+                                const double residual)
 {
   const std::string prefix = solverName;
 
@@ -323,7 +333,7 @@ void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParame
   actionPar.M5 = par.M5;
   actionPar.mass = par.ms;
   actionPar.scale = par.scale;
-  actionPar.boundary = "1 1 1 1";
+  actionPar.boundary = boundary;
   actionPar.twist = "0. 0. 0. 0.";
   app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
   actionPar.gauge = prefix + "_gauge_fp32";
@@ -344,42 +354,93 @@ void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParame
 
 void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName,
                                   const std::string gaugeName, const std::string gaugeTransform,
-                                  const double residual)
+                                  const std::string boundary, const double residual)
 {
   RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0LCDPar, solverName, gaugeName, gaugeTransform,
-                             residual);
+                             boundary, residual);
 }
 
 void RbcUkqcd::addC1MStrangeSolver(Application &app, const std::string solverName,
                                    const std::string gaugeName, const std::string gaugeTransform,
-                                   const double residual)
+                                   const std::string boundary, const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1mLCDPar, solverName, gaugeName, gaugeTransform,
-                             residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1mIRLPar, solverName, gaugeName, gaugeTransform,
+                             boundary, residual);
 }
 
 void RbcUkqcd::addC1M16StrangeSolver(Application &app, const std::string solverName,
                                      const std::string gaugeName, const std::string gaugeTransform,
-                                     const double residual)
+                                     const std::string boundary, const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m16LCDPar, solverName, gaugeName, gaugeTransform,
-                             residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m16IRLPar, solverName, gaugeName, gaugeTransform,
+                             boundary, residual);
 }
 
 void RbcUkqcd::addC1M20StrangeSolver(Application &app, const std::string solverName,
                                      const std::string gaugeName, const std::string gaugeTransform,
-                                     const double residual)
+                                     const std::string boundary, const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m20LCDPar, solverName, gaugeName, gaugeTransform,
-                             residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m20IRLPar, solverName, gaugeName, gaugeTransform,
+                             boundary, residual);
 }
 
 void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverName,
                                      const std::string gaugeName, const std::string gaugeTransform,
-                                     const double residual)
+                                     const std::string boundary, const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32LCDPar, solverName, gaugeName, gaugeTransform,
-                             residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, gaugeTransform,
+                             boundary, residual);
+}
+
+// Charm M0
+void RbcUkqcd::addM0CharmSolver(Application &app, const std::string solverName,
+                                const std::string gaugeName, const std::string gaugeTransform,
+                                const double mass, const std::string boundary,
+                                const double residual)
+{
+  const std::string prefix = solverName;
+
+  // stout smearing
+  MGauge::StoutSmearing::Par smearPar;
+
+  smearPar.gauge = gaugeName;
+  smearPar.steps = 3;
+  smearPar.rho = 0.1;
+  smearPar.orthogDim = "";
+  app.createModule<MGauge::StoutSmearing>(prefix + "_gauge_3stout", smearPar);
+
+  // Gauge field FP32 cast
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+
+  gaugeCastPar.field = prefix + "_gauge_3stout";
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_3stout_fp32",
+                                                         gaugeCastPar);
+
+  // Scaled DWF action + FP32 version
+  MAction::ScaledDWF::Par actionPar;
+
+  actionPar.gauge = prefix + "_gauge_3stout";
+  actionPar.Ls = m0UnitaryPar.Ls;
+  actionPar.M5 = 1.;
+  actionPar.mass = mass;
+  actionPar.scale = m0UnitaryPar.scale;
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_3stout_fp32";
+  app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  // Mixed-precision red-black preconditionned CG
+  MSolver::MixedPrecisionRBPrecCG::Par solverPar;
+
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 30000;
+  solverPar.maxOuterIteration = 100;
+  solverPar.residual = residual;
+  solverPar.innerGuesser = "";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
 }
 
 } // namespace hadpresets