From 97169d9f2fd5cc3bf421f600ffd71edcb497ce3a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 8 Mar 2024 14:21:15 +0900 Subject: [PATCH] added C0 presets --- RbcUkqcd.hpp | 154 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 150 insertions(+), 4 deletions(-) diff --git a/RbcUkqcd.hpp b/RbcUkqcd.hpp index 60b7212..92b15ec 100644 --- a/RbcUkqcd.hpp +++ b/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, c0ZMobiusLs> c0ZMobiusOmega{ + std::complex(1.458064389850479e+00, -0.000000000000000e+00), + std::complex(1.182313183893475e+00, -0.000000000000000e+00), + std::complex(8.309511666859551e-01, -0.000000000000000e+00), + std::complex(5.423524091567911e-01, -0.000000000000000e+00), + std::complex(3.419850204537295e-01, -0.000000000000000e+00), + std::complex(2.113790261902896e-01, -0.000000000000000e+00), + std::complex(1.260742995029118e-01, -0.000000000000000e+00), + std::complex(9.901366519626265e-02, -0.000000000000000e+00), + std::complex(6.863249884465925e-02, 5.506585308274019e-02), + std::complex(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, @@ -102,7 +127,11 @@ struct RbcUkqcd const std::string solverName, const std::string gaugeName, const std::string gaugeTransform, const std::string boundary, const double residual); - + static inline void addC0StrangeSolver(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 addM0StrangeSolver(Application &app, const std::string solverName, const std::string gaugeName, const std::string gaugeTransform, @@ -130,6 +159,10 @@ struct RbcUkqcd 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 std::string gaugeTransform, + 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", @@ -138,6 +171,72 @@ struct RbcUkqcd // 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; + + // Gauge field FP32 cast + MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar; + + gaugeCastPar.field = gaugeName; + app.createModule(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.ml; + actionPar.omega = std::vector>(RbcUkqcd::c0ZMobiusOmega.begin(), + RbcUkqcd::c0ZMobiusOmega.end()); + 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::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 = false; + epPar.gaugeXform = gaugeTransform; + app.createModule(prefix + "_epack", epPar); + + // Inner guesser + MGuesser::CoarseDeflation250F::Par iguessPar; + + iguessPar.eigenPack = prefix + "_epack"; + iguessPar.size = 2000; + app.createModule(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(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, @@ -352,11 +451,19 @@ void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParame app.createModule(solverName, solverPar); } +void RbcUkqcd::addC0StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c0UnitaryPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName, const std::string gaugeName, const std::string gaugeTransform, const std::string boundary, const double residual) { - RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0LCDPar, solverName, gaugeName, gaugeTransform, + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0UnitaryPar, solverName, gaugeName, gaugeTransform, boundary, residual); } @@ -392,6 +499,45 @@ void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverN boundary, residual); } +// Charm C0 +void RbcUkqcd::addC0CharmSolver(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(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(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 std::string gaugeTransform, @@ -413,10 +559,10 @@ void RbcUkqcd::addM0CharmSolver(Application &app, const std::string solverName, 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(prefix + "_dwf", actionPar);