Merge pull request 'Feature/iblep workflow presets' (#3) from RChHill/HadronsPresets:feature/iblep-workflow-presets into main

Reviewed-on: portelli/HadronsPresets#3

Jlqcd.hpp

@@ -0,0 +1,435 @@
+/*
+ * Copyright (C) 2026
+ *
+ * Author: Antonin Portelli <antonin.portelli@me.com>
+ * Author: Ryan Hill <Ryan.Hill@ed.ac.uk>
+ * Elements based on workflows from Teseo San Jose Perez
+ *
+ * Hadrons is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * Hadrons is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Hadrons.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * See the full license in the file "LICENSE" in the top level distribution
+ * directory.
+ */
+
+#pragma once
+
+#include <Hadrons/Application.hpp>
+#include <Hadrons/Modules.hpp>
+
+#include <array>
+#include <stdexcept>
+#include <string>
+#include <string_view>
+
+namespace hadpresets
+{
+using namespace Grid::Hadrons;
+
+struct Jlqcd
+{
+  struct EnsembleParameters
+  {
+    double ml, ms, M5, scale;
+    unsigned int L, T, Ls;
+  };
+
+  inline static constexpr EnsembleParameters fUd3SaPar {0.003,  0.015, 1., 2., 64, 128, 8};
+  inline static constexpr EnsembleParameters mUd3SaPar {0.0042, 0.025, 1., 2., 48, 96,  8};
+  inline static constexpr EnsembleParameters mUd3SbPar {0.0042, 0.018, 1., 2., 48, 96,  8};
+  inline static constexpr EnsembleParameters cUd2SaPar {0.0035, 0.040, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters cUd2SaLPar{0.0035, 0.040, 1., 2., 48, 96,  12};
+  inline static constexpr EnsembleParameters cUd3SaPar {0.007,  0.040, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters cUd3SbPar {0.007,  0.030, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters testPar   {0.010,  0.030, 1., 2.,  8, 16,   8};
+
+  static inline const EnsembleParameters &ensemble(const std::string& name);
+
+  static inline void addMixedPrecisionSolver(Application &app,
+                                             const EnsembleParameters &par,
+                                             const std::string &solverName,
+                                             const std::string &gaugeName,
+                                             const double mass,
+                                             const std::string &twist,
+                                             const std::string &boundary = "1 1 1 -1",
+                                             const double residual = 1.0e-8);
+
+  static inline void addNoFailSolver(Application &app, const EnsembleParameters &par,
+                                     const std::string &solverName,
+                                     const std::string &gaugeName, const double mass,
+                                     const std::string &twist,
+                                     const std::string &boundary = "1 1 1 -1",
+                                     const double residual = 1.0e-20);
+
+  // Light solvers (no deflation)
+  static inline void addLightSolver(Application &app, const EnsembleParameters &par,
+                                    const std::string &solverName, const std::string &gaugeName,
+                                    const double mass, const std::string &twist,
+                                    const std::string &boundary = "1 1 1 -1",
+                                    const double residual = 1.0e-8);
+
+  static inline void addLightSolver(Application &app, const EnsembleParameters &par,
+                                    const std::string &solverName, const std::string &gaugeName,
+                                    const std::string &twist = "0. 0. 0. 0.",
+                                    const std::string &boundary = "1 1 1 -1",
+                                    const double residual = 1.0e-8);
+
+  static inline void addfUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SbLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SbLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+
+  // Strange solvers
+  static inline void addStrangeSolver(Application &app, const EnsembleParameters &par,
+                                      const std::string &solverName,
+                                      const std::string &gaugeName, const double mass,
+                                      const std::string &twist,
+                                      const std::string &boundary = "1 1 1 -1",
+                                      const double residual = 1.0e-8);
+
+  static inline void addStrangeSolver(Application &app, const EnsembleParameters &par,
+                                      const std::string &solverName,
+                                      const std::string &gaugeName,
+                                      const std::string &twist = "0. 0. 0. 0.",
+                                      const std::string &boundary = "1 1 1 -1",
+                                      const double residual = 1.0e-8);
+
+  static inline void addfUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+
+  // Heavy solvers
+  static inline void addHeavySolver(Application &app, const EnsembleParameters &par,
+                                    const std::string &solverName, const std::string &gaugeName,
+                                    const double mass, const std::string &twist,
+                                    const std::string &boundary = "1 1 1 -1",
+                                    const double residual = 1.0e-20);
+
+  static inline void addfUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addmUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addmUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd2SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd2SaLHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+
+};
+
+inline const Jlqcd::EnsembleParameters &Jlqcd::ensemble(const std::string& name)
+{
+  if (name == "fUd3Sa" ) return fUd3SaPar;
+  if (name == "mUd3Sa" ) return mUd3SaPar;
+  if (name == "mUd3Sb" ) return mUd3SbPar;
+  if (name == "cUd2Sa" ) return cUd2SaPar;
+  if (name == "cUd2SaL") return cUd2SaLPar;
+  if (name == "cUd3Sa" ) return cUd3SaPar;
+  if (name == "cUd3Sb" ) return cUd3SbPar;
+  if (name == "test"   ) return testPar;
+
+  throw std::invalid_argument("unknown JLQCD ensemble preset '" + std::string(name) + "'");
+}
+
+inline void Jlqcd::addMixedPrecisionSolver(Application &app, const EnsembleParameters &par,
+                                           const std::string &solverName, const std::string &gaugeName,
+                                           const double mass, const std::string &twist,
+                                           const std::string &boundary, const double residual)
+{
+  const std::string prefix = solverName;
+
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+  gaugeCastPar.field = gaugeName;
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+
+  MAction::ScaledDWF::Par actionPar;
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = par.Ls;
+  actionPar.M5 = par.M5;
+  actionPar.mass = mass;
+  actionPar.scale = par.scale;
+  actionPar.boundary = boundary;
+  actionPar.twist = twist;
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+
+  actionPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  MSolver::MixedPrecisionRBPrecCG::Par solverPar;
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 10000;
+  solverPar.maxOuterIteration = 100;
+  solverPar.residual = residual;
+  solverPar.innerGuesser = "";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
+}
+
+inline void Jlqcd::addNoFailSolver(Application &app, const EnsembleParameters &par,
+                                   const std::string &solverName, const std::string &gaugeName,
+                                   const double mass, const std::string &twist,
+                                   const std::string &boundary, const double residual)
+{
+  const std::string prefix = solverName;
+
+  MAction::ScaledDWF::Par actionPar;
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = par.Ls;
+  actionPar.M5 = par.M5;
+  actionPar.mass = mass;
+  actionPar.scale = par.scale;
+  actionPar.boundary = boundary;
+  actionPar.twist = twist;
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+
+  MSolver::RBPrecCGNoFail::Par solverPar;
+  solverPar.action = prefix + "_dwf";
+  solverPar.maxIteration = 30000;
+  solverPar.residual = residual;
+  solverPar.guesser = "";
+  app.createModule<MSolver::RBPrecCGNoFail>(solverName, solverPar);
+}
+
+// Light solvers (no deflation)
+inline void Jlqcd::addLightSolver(Application &app, const EnsembleParameters &par,
+                                  const std::string &solverName, const std::string &gaugeName,
+                                  const double mass, const std::string &twist,
+                                  const std::string &boundary, const double residual)
+{
+  addMixedPrecisionSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
+}
+
+inline void Jlqcd::addLightSolver(Application &app, const EnsembleParameters &par,
+                                  const std::string &solverName, const std::string &gaugeName,
+                                  const std::string &twist, const std::string &boundary,
+                                  const double residual)
+{
+  addLightSolver(app, par, solverName, gaugeName, par.ml, twist, boundary, residual);
+}
+
+inline void Jlqcd::addfUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, fUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, mUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SbLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, mUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd2SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd2SaLPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SbLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+// Strange solvers
+inline void Jlqcd::addStrangeSolver(Application &app, const EnsembleParameters &par,
+                                    const std::string &solverName, const std::string &gaugeName,
+                                    const double mass, const std::string &twist,
+                                    const std::string &boundary, const double residual)
+{
+    addMixedPrecisionSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
+}
+
+inline void Jlqcd::addStrangeSolver(Application &app, const EnsembleParameters &par,
+                                    const std::string &solverName, const std::string &gaugeName,
+                                    const std::string &twist, const std::string &boundary,
+                                    const double residual)
+{
+    addStrangeSolver(app, par, solverName, gaugeName, par.ms, twist, boundary, residual);
+}
+
+inline void Jlqcd::addfUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, fUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, mUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, mUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd2SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLStrangeSolver(Application &app, const std::string &solverName,
+                                           const std::string &gaugeName, const std::string &twist,
+                                           const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd2SaLPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+// Heavy solvers
+inline void Jlqcd::addHeavySolver(Application &app, const EnsembleParameters &par,
+                                  const std::string &solverName, const std::string &gaugeName,
+                                  const double mass, const std::string &twist,
+                                  const std::string &boundary, const double residual)
+{
+  addNoFailSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
+}
+
+inline void Jlqcd::addfUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, fUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addmUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, mUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addmUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, mUd3SbPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd2SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd2SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd2SaLHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                         const double mass, const std::string &twist,
+                                         const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd2SaLPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd3SbPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+
+} // namespace hadpresets

RbcUkqcd.hpp

@@ -29,6 +29,7 @@
 namespace hadpresets
 {
 using namespace Grid::Hadrons;
+using namespace std::complex_literals;
 
 struct RbcUkqcd
 {
@@ -39,25 +40,168 @@ struct RbcUkqcd
     unsigned int L, T, Ls;
   };
 
-  static constexpr EnsembleParameters m0UnitaryPar{0.000678, 0.02661, 1.8, 2., 64, 128, 12};
+  inline static constexpr hadpresets::RbcUkqcd::EnsembleParameters testPar{0.005, 0.0362, 1., 2.,
-  static constexpr EnsembleParameters m0LCDPar{0.0006203, 0.02661, 1.8, 2., 64, 128, 12};
+                                                                           8,     8,      4};
 
-  // Solvers (LCD: Local Coherence Deflation)
+  // C0
-  static inline void addM0LightLCDSolver(Application &app, const std::string solverName,
+  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};
+
+  // C0L
+  inline static constexpr EnsembleParameters c0LUnitaryPar{0.00078, 0.0362, 1.8, 2., 64, 128, 24};
+  inline static constexpr EnsembleParameters c0LLCDPar{0.00078, 0.0362, 1.8, 2., 64, 128, 24};
+  
+  // 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};
+
+  // C1M
+  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
+  {
+    double alpha, beta;
+    unsigned int nPoly, nStop, nK, nM;
+  };
+
+  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{5.0e-05, 5.5, 101, 100, 110, 120};
+
+  // Light solvers: no deflation
+  static inline void addLightSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
+                                    const std::string solverName, const std::string gaugeName,
+                                    const std::string boundary, const double residual);
+
+  // 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
+  addC0LightMADWFLCDSolver(Grid::Hadrons::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
+  addC0LLightLCDSolver(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,
+                      const std::string boundary = "1 1 1 1", const double residual = 1.0e-8);
+
+  // 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 boundary, const double residual);
+
+  static inline void addC1MLightRuntimeIRLSolver(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 addC1M16LightRuntimeIRLSolver(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 addC1M20LightRuntimeIRLSolver(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 addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                                   const std::string gaugeName,
+                                                   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 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 addC0StrangeZMobiusSolver(Grid::Hadrons::Application &app,
+                                               const std::string solverName,
+                                               const std::string gaugeName,
+                                               const std::string boundary, const double residual);
+  static inline void addC0LStrangeSolver(Application &app, const std::string solverName,
                                          const std::string gaugeName,
-                                         const std::string gaugeTransform,
+                                         const std::string boundary = "1 1 1 -1",
-                                         const std::string eigenpackPath,
                                          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 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 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 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 boundary = "1 1 1 -1",
+                                           const double residual = 1.0e-8);
+
+  // Charm solvers, mass is a free parameter
+  static inline void addUnsmearedCharmSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName, const double mass,
+                                             const int Ls, const double M5,
+                                             const std::string boundary, const double residual);
+  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 addC0LCharmSolver(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 double mass,
+                                      const std::string boundary = "1 1 1 1",
+                                      const double residual = 1.0e-18);
 };
 
 // Implementations /////////////////////////////////////////////////////////////////////////////////
-void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName,
+
-                                   const std::string gaugeName, const std::string gaugeTransform,
+// Light
-                                   const std::string eigenpackPath, const double residual)
+void RbcUkqcd::addLightSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
+                              const std::string solverName, const std::string gaugeName,
+                              const std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -70,12 +214,345 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam
   // Scaled DWF action + FP32 version
   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 = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_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);
+}
+
+// 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);
+}
+
+void RbcUkqcd::addC0LightMADWFLCDSolver(Grid::Hadrons::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)
+{
+  using namespace Grid;
+  using namespace Hadrons;
+  using namespace hadpresets;
+
+  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
+  MAction::ScaledDWF::Par actionPar;
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = RbcUkqcd::c0UnitaryPar.Ls;
+  actionPar.M5 = RbcUkqcd::c0UnitaryPar.M5;
+  actionPar.mass = RbcUkqcd::c0UnitaryPar.ml;
+  actionPar.scale = 2.;
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+
+  // Single-prec ZMobius Action
+  MAction::ZMobiusDWFF::Par actionFPar;
+  actionFPar.Ls = RbcUkqcd::c0LCDPar.Ls;
+  actionFPar.M5 = RbcUkqcd::c0LCDPar.M5;
+  actionFPar.mass = RbcUkqcd::c0LCDPar.ml;
+  actionFPar.b = 1.;
+  actionFPar.c = 0.;
+  actionFPar.omega = std::vector<std::complex<double>>(RbcUkqcd::c0ZMobiusOmega.begin(),
+                                                       RbcUkqcd::c0ZMobiusOmega.end());
+  actionFPar.boundary = boundary;
+  actionFPar.twist = "0. 0. 0. 0.";
+  actionFPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ZMobiusDWFF>(prefix + "_dwf_fp32", actionFPar);
+
+  // 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);
+
+  // ZMobius-MADWF red-black preconditionned CG
+  MSolver::ZMADWFMixedPrecCG::Par solverPar;
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 10000;
+  solverPar.maxPVIteration = 1000;
+  solverPar.maxOuterIteration = 100;
+  solverPar.innerResidual = 1e-5;
+  solverPar.PVResidual = 1e-5;
+  solverPar.outerResidual = residual;
+  solverPar.guesser = prefix + "_iguesser";
+  app.createModule<MSolver::ZMADWFMixedPrecCG>(solverName, solverPar);
+}
+
+void RbcUkqcd::addC0StrangeZMobiusSolver(Grid::Hadrons::Application &app,
+                                         const std::string solverName, const std::string gaugeName,
+                                         const std::string boundary, const double residual)
+{
+  using namespace Grid;
+  using namespace Hadrons;
+  using namespace hadpresets;
+
+  const std::string prefix = solverName;
+
+  // Gauge field & transform FP32 cast
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+
+  gaugeCastPar.field = gaugeName;
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+
+  // 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.ms;
+  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);
+
+  // Batched mixed-precision red-black preconditionned CG
+  MSolver::ZMixedPrecisionRBPrecCG::Par solverPar;
+
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 10000;
+  solverPar.maxOuterIteration = 100;
+  solverPar.residual = residual;
+  solverPar.innerGuesser = "";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::ZMixedPrecisionRBPrecCG>(solverName, solverPar);
+}
+
+// Light C0L (load deflation from disk)
+void RbcUkqcd::addC0LLightLCDSolver(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 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;
+
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = RbcUkqcd::c0LLCDPar.Ls;
+  actionPar.M5 = RbcUkqcd::c0LLCDPar.M5;
+  actionPar.mass = RbcUkqcd::c0LLCDPar.ml;
+  actionPar.scale = RbcUkqcd::c0LLCDPar.scale;
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  // Compressed eigenpack
+  MIO::LoadCoarseFermionEigenPack200F::Par epPar;
+
+  epPar.filestem = eigenpackPath;
+  epPar.multiFile = true;
+  epPar.redBlack = true;
+  epPar.sizeFine = 200;
+  epPar.sizeCoarse = 4500;
+  epPar.Ls = 24;
+  epPar.blockSize = "4 4 4 4 24";
+  epPar.orthogonalise = false;
+  epPar.gaugeXform = gaugeFixed ? (prefix + "_gaugeTransform_fp32") : "";
+  app.createModule<MIO::LoadCoarseFermionEigenPack200F>(prefix + "_epack", epPar);
+
+  // Inner guesser
+  MGuesser::CoarseDeflation200F::Par iguessPar;
+
+  iguessPar.eigenPack = prefix + "_epack";
+  iguessPar.size = 4500;
+  app.createModule<MGuesser::CoarseDeflation200F>(prefix + "_iguesser", iguessPar);
+
+  // Batched mixed-precision red-black preconditionned CG
+  MSolver::MixedPrecisionRBPrecCGBatched::Par solverPar;
+
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 250;
+  solverPar.maxOuterIteration = 100;
+  solverPar.maxPatchupIteration = 10000;
+  solverPar.residual = residual;
+  solverPar.updateResidual = true;
+  solverPar.innerGuesser = prefix + "_iguesser";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::MixedPrecisionRBPrecCGBatched>(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,
+                                   const std::string eigenpackPath, const std::string boundary,
+                                   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;
+
   actionPar.gauge = gaugeName;
   actionPar.Ls = RbcUkqcd::m0LCDPar.Ls;
   actionPar.M5 = RbcUkqcd::m0LCDPar.M5;
   actionPar.mass = RbcUkqcd::m0LCDPar.ml;
   actionPar.scale = RbcUkqcd::m0LCDPar.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";
@@ -92,7 +569,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
@@ -117,9 +594,11 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam
   app.createModule<MSolver::MixedPrecisionRBPrecCGBatched>(solverName, solverPar);
 }
 
-void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName,
+// Deflation at runtime
-                                  const std::string gaugeName, const std::string gaugeTransform,
+void RbcUkqcd::addLightRuntimeIRLSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
-                                  const double residual)
+                                        const RbcUkqcd::DeflationParameters &deflPar,
+                                        const std::string solverName, const std::string gaugeName,
+                                        const std::string boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -133,17 +612,120 @@ void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName
   MAction::ScaledDWF::Par actionPar;
 
   actionPar.gauge = gaugeName;
-  actionPar.Ls = RbcUkqcd::m0LCDPar.Ls;
+  actionPar.Ls = par.Ls;
-  actionPar.M5 = RbcUkqcd::m0LCDPar.M5;
+  actionPar.M5 = par.M5;
-  actionPar.mass = RbcUkqcd::m0LCDPar.ms;
+  actionPar.mass = par.ml;
-  actionPar.scale = RbcUkqcd::m0LCDPar.scale;
+  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";
   app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
 
-  // Batched mixed-precision red-black preconditionned CG
+  // Guesser
+  MFermion::Operators::Par opPar;
+  opPar.action = prefix + "_dwf";
+  app.createModule<MFermion::Operators>(prefix + "_dwf_op", opPar);
+
+  MSolver::FermionImplicitlyRestartedLanczosIo32::Par lanPar;
+  lanPar.op = prefix + "_dwf_op_schur";
+  lanPar.multiFile = false;
+  lanPar.redBlack = true;
+  lanPar.lanczosParams.Cheby.alpha = deflPar.alpha;
+  lanPar.lanczosParams.Cheby.beta = deflPar.beta;
+  lanPar.lanczosParams.Cheby.Npoly = deflPar.nPoly;
+  lanPar.lanczosParams.Nstop = deflPar.nStop;
+  lanPar.lanczosParams.Nk = deflPar.nK;
+  lanPar.lanczosParams.Nm = deflPar.nM;
+  lanPar.lanczosParams.resid = 3e-6;
+  lanPar.lanczosParams.MaxIt = 10000;
+  lanPar.lanczosParams.betastp = 0;
+  lanPar.lanczosParams.MinRes = 0;
+  lanPar.output = "";
+  app.createModule<MSolver::FermionImplicitlyRestartedLanczosIo32>(prefix + "_epack", lanPar);
+
+  MGuesser::ExactDeflation::Par guessPar;
+  guessPar.eigenPack = prefix + "_epack";
+  guessPar.size = deflPar.nStop;
+  app.createModule<MGuesser::ExactDeflation>(prefix + "_defl", guessPar);
+
+  // 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 = prefix + "_defl";
+  app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
+}
+
+// Light C1M
+void RbcUkqcd::addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                           const std::string gaugeName, const std::string boundary,
+                                           const double residual)
+{
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1mIRLPar, RbcUkqcd::c1mDeflPar, solverName,
+                                     gaugeName, boundary, residual);
+}
+
+// Light C1M16
+void RbcUkqcd::addC1M16LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string boundary, const double residual)
+{
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m16IRLPar, RbcUkqcd::c1m16DeflPar, solverName,
+                                     gaugeName, boundary, residual);
+}
+
+// Light C1M20
+void RbcUkqcd::addC1M20LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string boundary, const double residual)
+{
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m20IRLPar, RbcUkqcd::c1m20DeflPar, solverName,
+                                     gaugeName, boundary, residual);
+}
+
+// Light C1M32
+void RbcUkqcd::addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName,
+                                             const std::string gaugeName,
+                                             const std::string boundary, const double residual)
+{
+  RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m32IRLPar, RbcUkqcd::c1m32DeflPar, solverName,
+                                     gaugeName, boundary, residual);
+}
+
+// Strange
+void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par,
+                                const std::string solverName, const std::string gaugeName,
+                                const std::string boundary, const double residual)
+{
+  const std::string prefix = solverName;
+
+  // Gauge field FP32 cast
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+
+  gaugeCastPar.field = gaugeName;
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+
+  // Scaled DWF action + FP32 version
+  MAction::ScaledDWF::Par actionPar;
+
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = par.Ls;
+  actionPar.M5 = par.M5;
+  actionPar.mass = par.ms;
+  actionPar.scale = par.scale;
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  // Mixed-precision red-black preconditionned CG
   MSolver::MixedPrecisionRBPrecCG::Par solverPar;
 
   solverPar.innerAction = prefix + "_dwf_fp32";
@@ -156,4 +738,215 @@ void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName
   app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
 }
 
+void RbcUkqcd::addC0StrangeSolver(Application &app, const std::string solverName,
+                                  const std::string gaugeName, const std::string boundary,
+                                  const double residual)
+{
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c0UnitaryPar, solverName, gaugeName, boundary,
+                             residual);
+}
+
+void RbcUkqcd::addC0LStrangeSolver(Application &app, const std::string solverName,
+                                   const std::string gaugeName, const std::string boundary,
+                                   const double residual)
+{
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c0LUnitaryPar, 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 boundary,
+                                   const double 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 boundary,
+                                     const double 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 boundary,
+                                     const double 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 boundary,
+                                     const double residual)
+{
+  RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, boundary, residual);
+}
+
+// Charm
+void RbcUkqcd::addUnsmearedCharmSolver(Grid::Hadrons::Application &app,
+                                       const std::string solverName, const std::string gaugeName,
+                                       const double mass, const int Ls, const double M5,
+                                       const std::string boundary, const double residual)
+{
+  using namespace Grid;
+  using namespace Hadrons;
+
+  const std::string prefix = solverName;
+
+  // Gauge field & transform FP32 cast
+  MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar;
+
+  gaugeCastPar.field = gaugeName;
+  app.createModule<MUtilities::GaugeSinglePrecisionCast>(prefix + "_gauge_fp32", gaugeCastPar);
+
+  // Scaled DWF action + FP32 version
+  MAction::ScaledDWF::Par actionPar;
+
+  actionPar.gauge = gaugeName;
+  actionPar.Ls = Ls;
+  actionPar.M5 = M5;
+  actionPar.mass = mass;
+  actionPar.scale = 2.;
+
+  actionPar.boundary = boundary;
+  actionPar.twist = "0. 0. 0. 0.";
+  app.createModule<MAction::ScaledDWF>(prefix + "_dwf", actionPar);
+  actionPar.gauge = prefix + "_gauge_fp32";
+  app.createModule<MAction::ScaledDWFF>(prefix + "_dwf_fp32", actionPar);
+
+  // Batched mixed-precision red-black preconditionned CG
+  MSolver::MixedPrecisionRBPrecCG::Par solverPar;
+
+  solverPar.innerAction = prefix + "_dwf_fp32";
+  solverPar.outerAction = prefix + "_dwf";
+  solverPar.maxInnerIteration = 10000;
+  solverPar.maxOuterIteration = 100;
+  solverPar.residual = residual;
+  solverPar.innerGuesser = "";
+  solverPar.outerGuesser = "";
+  app.createModule<MSolver::MixedPrecisionRBPrecCG>(solverName, solverPar);
+}
+
+// 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;
+
+  // 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 = 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 C0L
+void RbcUkqcd::addC0LCharmSolver(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 = 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 = 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);
+}
+
 } // namespace hadpresets