more BC fix

Reviewed-on: portelli/HadronsPresets#2

RbcUkqcd.hpp

@@ -29,6 +29,7 @@
 namespace hadpresets
 {
 using namespace Grid::Hadrons;
+using namespace std::complex_literals;
 
 struct RbcUkqcd
 {
@@ -39,6 +40,23 @@ struct RbcUkqcd
     unsigned int L, T, Ls;
   };
 
+  // C0
+  inline static constexpr EnsembleParameters c0UnitaryPar{0.00078, 0.0362, 1.8, 2., 48, 96, 24};
+  inline static constexpr unsigned int c0ZMobiusLs = 10;
+  inline static constexpr std::array<std::complex<double>, c0ZMobiusLs> c0ZMobiusOmega{
+      std::complex<double>(1.458064389850479e+00, -0.000000000000000e+00),
+      std::complex<double>(1.182313183893475e+00, -0.000000000000000e+00),
+      std::complex<double>(8.309511666859551e-01, -0.000000000000000e+00),
+      std::complex<double>(5.423524091567911e-01, -0.000000000000000e+00),
+      std::complex<double>(3.419850204537295e-01, -0.000000000000000e+00),
+      std::complex<double>(2.113790261902896e-01, -0.000000000000000e+00),
+      std::complex<double>(1.260742995029118e-01, -0.000000000000000e+00),
+      std::complex<double>(9.901366519626265e-02, -0.000000000000000e+00),
+      std::complex<double>(6.863249884465925e-02, 5.506585308274019e-02),
+      std::complex<double>(6.863249884465925e-02, -5.506585308274019e-02)};
+  inline static constexpr EnsembleParameters c0LCDPar{0.00078, 0.0362, 1.8,        NAN,
+                                                      48,      96,     c0ZMobiusLs};
+
   // M0
   inline static constexpr EnsembleParameters m0UnitaryPar{0.000678, 0.02661, 1.8, 2., 64, 128, 12};
   inline static constexpr EnsembleParameters m0LCDPar{0.0006203, 0.02661, 1.8, 2., 64, 128, 12};
@@ -62,6 +80,13 @@ struct RbcUkqcd
   inline static constexpr DeflationParameters c1m32DeflPar{5.0e-05, 5.5, 101, 100, 110, 120};
 
   // Light solvers: load deflation from disk
+  static inline void addC0LightZMobiusLCDSolver(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
   addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName,
                       const std::string gaugeTransform, const std::string eigenpackPath,
@@ -73,71 +98,144 @@ struct RbcUkqcd
                                               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 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 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 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 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 std::string boundary,
-                                      const double residual);
-
+                                      const std::string boundary, const double residual);
+  static inline void addC0StrangeSolver(Application &app, const std::string solverName,
+                                        const std::string gaugeName,
+                                        const std::string boundary = "1 1 1 -1",
+                                        const double residual = 1.0e-8);
   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 addC0CharmSolver(Application &app, const std::string solverName,
+                                      const std::string gaugeName, const double mass,
+                                      const std::string boundary = "1 1 1 -1",
+                                      const double residual = 1.0e-18);
   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 std::string gaugeName, const double mass,
+                                      const std::string boundary = "1 1 1 1",
                                       const double residual = 1.0e-18);
 };
 
 // Implementations /////////////////////////////////////////////////////////////////////////////////
 
+// Light C0 (load deflation from disk)
+void RbcUkqcd::addC0LightZMobiusLCDSolver(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;
+  const bool gaugeFixed = !gaugeTransform.empty();
+
+  // Gauge field & transform FP32 cast
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+
+  gaugeCastPar.field = gaugeName;
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+  if (gaugeFixed)
+  {
+    MUtilities::ColourMatrixSinglePrecisionCast::Par transformCastPar;
+
+    transformCastPar.field = gaugeTransform;
+    app.createModule<MUtilities::ColourMatrixSinglePrecisionCast>(prefix + "_gaugeTransform_fp32",
+                                                                  transformCastPar);
+  }
+
+  // Scaled DWF action + FP32 version
+  MAction::ZMobiusDWF::Par actionPar;
+
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = RbcUkqcd::c0LCDPar.Ls;
+  actionPar.M5 = RbcUkqcd::c0LCDPar.M5;
+  actionPar.mass = RbcUkqcd::c0LCDPar.ml;
+  actionPar.b = 1.;
+  actionPar.c = 0.;
+  actionPar.omega = std::vector<std::complex<double>>(RbcUkqcd::c0ZMobiusOmega.begin(),
+                                                      RbcUkqcd::c0ZMobiusOmega.end());
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ZMobiusDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ZMobiusDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  // Compressed eigenpack
+  MIO::LoadCoarseFermionEigenPack250F::Par epPar;
+
+  epPar.filestem = eigenpackPath;
+  epPar.multiFile = true;
+  epPar.redBlack = true;
+  epPar.sizeFine = 250;
+  epPar.sizeCoarse = 2000;
+  epPar.Ls = 10;
+  epPar.blockSize = "4 4 4 3 10";
+  epPar.orthogonalise = gaugeFixed;
+  epPar.gaugeXform = gaugeFixed ? (prefix + "_gaugeTransform_fp32") : "";
+  app.createModule<MIO::LoadCoarseFermionEigenPack250F>(prefix + "_epack", epPar);
+
+  // Inner guesser
+  MGuesser::CoarseDeflation250F::Par iguessPar;
+
+  iguessPar.eigenPack = prefix + "_epack";
+  iguessPar.size = 2000;
+  app.createModule<MGuesser::CoarseDeflation250F>(prefix + "_iguesser", iguessPar);
+
+  // Batched mixed-precision red-black preconditionned CG
+  MSolver::ZMixedPrecisionRBPrecCGBatched::Par solverPar;
+
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 300;
+  solverPar.maxOuterIteration = 100;
+  solverPar.maxPatchupIteration = 1000;
+  solverPar.residual = residual;
+  solverPar.updateResidual = true;
+  solverPar.innerGuesser = prefix + "_iguesser";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::ZMixedPrecisionRBPrecCGBatched>(solverName, solverPar);
+}
+
 // Light M0 (load deflation from disk)
 void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName,
                                    const std::string gaugeName, const std::string gaugeTransform,
@@ -145,12 +243,21 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam
                                    const double residual)
 {
   const std::string prefix = solverName;
+  const bool gaugeFixed = !gaugeTransform.empty();
 
   // Gauge field FP32 cast
   MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
 
   gaugeCastPar.field = gaugeName;
   app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+  if (gaugeFixed)
+  {
+    MUtilities::ColourMatrixSinglePrecisionCast::Par transformCastPar;
+
+    transformCastPar.field = gaugeTransform;
+    app.createModule<MUtilities::ColourMatrixSinglePrecisionCast>(prefix + "_gaugeTransform_fp32",
+                                                                  transformCastPar);
+  }
 
   // Scaled DWF action + FP32 version
   MAction::ScaledDWF::Par actionPar;
@@ -177,7 +284,7 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam
   epPar.Ls = 12;
   epPar.blockSize = "4 4 4 4 12";
   epPar.orthogonalise = false;
-  epPar.gaugeXform = gaugeTransform;
+  epPar.gaugeXform = gaugeFixed ? (prefix + "_gaugeTransform_fp32") : "";
   app.createModule<MIO::LoadCoarseFermionEigenPack250F>(prefix + "_epack", epPar);
 
   // Inner guesser
@@ -206,7 +313,6 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam
 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 std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
@@ -273,49 +379,44 @@ void RbcUkqcd::addLightRuntimeIRLSolver(Application &app, const RbcUkqcd::Ensemb
 
 // Light C1M
 void RbcUkqcd::addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName,
-                                           const std::string gaugeName,
-                                           const std::string gaugeTransform,
-                                           const std::string boundary, const double residual)
+                                           const std::string gaugeName, const std::string boundary,
+                                           const double residual)
 {
   RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1mIRLPar, RbcUkqcd::c1mDeflPar, solverName,
-                                     gaugeName, gaugeTransform, boundary, residual);
+                                     gaugeName, boundary, residual);
 }
 
 // Light C1M16
 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::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m16IRLPar, RbcUkqcd::c1m16DeflPar, solverName,
-                                     gaugeName, gaugeTransform, boundary, residual);
+                                     gaugeName, boundary, residual);
 }
 
 // Light C1M20
 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::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m20IRLPar, RbcUkqcd::c1m20DeflPar, solverName,
-                                     gaugeName, gaugeTransform, boundary, residual);
+                                     gaugeName, boundary, residual);
 }
 
 // Light C1M32
 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::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m32IRLPar, RbcUkqcd::c1m32DeflPar, solverName,
-                                     gaugeName, gaugeTransform, boundary, residual);
+                                     gaugeName, boundary, residual);
 }
 
 // Strange
 void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
                                 const std::string solverName, const std::string gaugeName,
-                                const std::string gaugeTransform, const std::string boundary,
-                                const double residual)
+                                const std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -352,51 +453,54 @@ void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParame
   app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
 }
 
-void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName,
-                                  const std::string gaugeName, const std::string gaugeTransform,
-                                  const std::string boundary, const double residual)
+void RbcUkqcd::addC0StrangeSolver(Application &app, const std::string solverName,
+                                  const std::string gaugeName, const std::string boundary,
+                                  const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0LCDPar, solverName, gaugeName, gaugeTransform,
-                             boundary, residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c0UnitaryPar, solverName, gaugeName, boundary,
+                             residual);
+}
+
+void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName,
+                                  const std::string gaugeName, const std::string boundary,
+                                  const double residual)
+{
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0UnitaryPar, solverName, gaugeName, boundary,
+                             residual);
 }
 
 void RbcUkqcd::addC1MStrangeSolver(Application &app, const std::string solverName,
-                                   const std::string gaugeName, const std::string gaugeTransform,
-                                   const std::string boundary, const double residual)
+                                   const std::string gaugeName, const std::string boundary,
+                                   const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1mIRLPar, solverName, gaugeName, gaugeTransform,
-                             boundary, residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1mIRLPar, solverName, gaugeName, boundary, residual);
 }
 
 void RbcUkqcd::addC1M16StrangeSolver(Application &app, const std::string solverName,
-                                     const std::string gaugeName, const std::string gaugeTransform,
-                                     const std::string boundary, const double residual)
+                                     const std::string gaugeName, const std::string boundary,
+                                     const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m16IRLPar, solverName, gaugeName, gaugeTransform,
-                             boundary, residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m16IRLPar, solverName, gaugeName, boundary, residual);
 }
 
 void RbcUkqcd::addC1M20StrangeSolver(Application &app, const std::string solverName,
-                                     const std::string gaugeName, const std::string gaugeTransform,
-                                     const std::string boundary, const double residual)
+                                     const std::string gaugeName, const std::string boundary,
+                                     const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m20IRLPar, solverName, gaugeName, gaugeTransform,
-                             boundary, residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m20IRLPar, solverName, gaugeName, boundary, residual);
 }
 
 void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverName,
-                                     const std::string gaugeName, const std::string gaugeTransform,
-                                     const std::string boundary, const double residual)
+                                     const std::string gaugeName, const std::string boundary,
+                                     const double residual)
 {
-  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, gaugeTransform,
-                             boundary, residual);
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, 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)
+// Charm C0
+void RbcUkqcd::addC0CharmSolver(Application &app, const std::string solverName,
+                                const std::string gaugeName, const double mass,
+                                const std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -409,38 +513,64 @@ void RbcUkqcd::addM0CharmSolver(Application &app, const std::string solverName,
   smearPar.orthogDim = "";
   app.createModule<MGauge::StoutSmearing>(prefix + "_gauge_3stout", smearPar);
 
-  // Gauge field FP32 cast
-  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+  // Scaled DWF action + FP32 version
+  MAction::ScaledDWF::Par actionPar;
 
-  gaugeCastPar.field = prefix + "_gauge_3stout";
-  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_3stout_fp32",
-                                                         gaugeCastPar);
+  actionPar.gauge = prefix + "_gauge_3stout";
+  actionPar.Ls = 12;
+  actionPar.M5 = 1.;
+  actionPar.mass = mass;
+  actionPar.scale = 2.;
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+
+  // Red-black preconditionned CG
+  MSolver::RBPrecCG::Par solverPar;
+
+  solverPar.action = prefix + "_dwf";
+  solverPar.maxIteration = 30000;
+  solverPar.residual = residual;
+  solverPar.guesser = "";
+  app.createModule<MSolver::RBPrecCGNoFail>(solverName, solverPar);
+}
+
+// Charm M0
+void RbcUkqcd::addM0CharmSolver(Application &app, const std::string solverName,
+                                const std::string gaugeName, 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);
 
   // Scaled DWF action + FP32 version
   MAction::ScaledDWF::Par actionPar;
 
   actionPar.gauge = prefix + "_gauge_3stout";
-  actionPar.Ls = m0UnitaryPar.Ls;
+  actionPar.Ls = 12;
   actionPar.M5 = 1.;
   actionPar.mass = mass;
-  actionPar.scale = m0UnitaryPar.scale;
+  actionPar.scale = 2.;
   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;
+  // Red-black preconditionned CG
+  MSolver::RBPrecCG::Par solverPar;
 
-  solverPar.innerAction = prefix + "_dwf_fp32";
-  solverPar.outerAction = prefix + "_dwf";
-  solverPar.maxInnerIteration = 30000;
-  solverPar.maxOuterIteration = 100;
+  solverPar.action = prefix + "_dwf";
+  solverPar.maxIteration = 30000;
   solverPar.residual = residual;
-  solverPar.innerGuesser = "";
-  solverPar.outerGuesser = "";
-  app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
+  solverPar.guesser = "";
+  app.createModule<MSolver::RBPrecCGNoFail>(solverName, solverPar);
 }
 
 } // namespace hadpresets