From 3936ed389eb63c39a965e6b7637c1799359b12e4 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 6 Jun 2023 19:17:45 +0100
Subject: [PATCH] M0 charm solver

---
 RbcUkqcd.hpp | 63 +++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 60 insertions(+), 3 deletions(-)

diff --git a/RbcUkqcd.hpp b/RbcUkqcd.hpp
index 95c05e9..f5d190c 100644
--- a/RbcUkqcd.hpp
+++ b/RbcUkqcd.hpp
@@ -61,13 +61,13 @@ struct RbcUkqcd
   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};
 
-  // Light solver: load deflation from disk
+  // Light solvers: load deflation from disk
   static inline void
   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);
 
-  // Light solver: deflation at runtime
+  // Light solvers: deflation at runtime
   static inline void addLightRuntimeIRLSolver(Application &app,
                                               const RbcUkqcd::EnsembleParameters &par,
                                               const RbcUkqcd::DeflationParameters &deflPar,
@@ -97,7 +97,7 @@ struct RbcUkqcd
                                                    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 std::string boundary,
@@ -128,6 +128,12 @@ struct RbcUkqcd
                                            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 /////////////////////////////////////////////////////////////////////////////////
@@ -386,4 +392,55 @@ void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverN
                              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