From d5a757128f6508ad769dfc4e4a8516fde6b02ef4 Mon Sep 17 00:00:00 2001 From: Ryan Hill Date: Tue, 29 Sep 2026 08:44:45 +0100 Subject: [PATCH] Add C0L support --- RbcUkqcd.hpp | 136 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/RbcUkqcd.hpp b/RbcUkqcd.hpp index a4a5f00..85f6f42 100644 --- a/RbcUkqcd.hpp +++ b/RbcUkqcd.hpp @@ -60,6 +60,10 @@ struct RbcUkqcd 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}; @@ -101,6 +105,11 @@ struct RbcUkqcd 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, @@ -143,6 +152,10 @@ struct RbcUkqcd 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 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 boundary = "1 1 1 1", @@ -173,6 +186,10 @@ struct RbcUkqcd 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", @@ -431,6 +448,79 @@ void RbcUkqcd::addC0StrangeZMobiusSolver(Grid::Hadrons::Application &app, app.createModule(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(prefix + "_gauge_fp32", gaugeCastPar); + if (gaugeFixed) + { + MUtilities::ColourMatrixSinglePrecisionCast::Par transformCastPar; + + transformCastPar.field = gaugeTransform; + app.createModule(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(prefix + "_dwf", actionPar); + actionPar.gauge = prefix + "_gauge_fp32"; + app.createModule(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(prefix + "_epack", epPar); + + // Inner guesser + MGuesser::CoarseDeflation200F::Par iguessPar; + + iguessPar.eigenPack = prefix + "_epack"; + iguessPar.size = 4500; + app.createModule(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(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, @@ -656,6 +746,14 @@ void RbcUkqcd::addC0StrangeSolver(Application &app, const std::string solverName 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) @@ -775,6 +873,44 @@ void RbcUkqcd::addC0CharmSolver(Application &app, const std::string solverName, app.createModule(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(prefix + "_gauge_3stout", smearPar); + + // Scaled DWF action + FP32 version + MAction::ScaledDWF::Par actionPar; + + actionPar.gauge = prefix + "_gauge_3stout"; + actionPar.Ls = 24; + actionPar.M5 = 1.; + actionPar.mass = mass; + actionPar.scale = 2.; + actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; + app.createModule(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(solverName, solverPar); +} + // Charm M0 void RbcUkqcd::addM0CharmSolver(Application &app, const std::string solverName, const std::string gaugeName, const double mass,