Feature/iblep workflow presets #3

Open
RChHill wants to merge 2 commits from RChHill/HadronsPresets:feature/iblep-workflow-presets into main
2 changed files with 446 additions and 0 deletions
+206
View File
@@ -0,0 +1,206 @@
/*
* Copyright (C) 2026
*
* Hadrons presets for the JLQCD scaled-domain-wall ensembles used by iblep.
*/
#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);
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 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 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);
};
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);
}
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::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::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);
}
} // namespace hadpresets
+240
View File
@@ -40,6 +40,9 @@ struct RbcUkqcd
unsigned int L, T, Ls; unsigned int L, T, Ls;
}; };
inline static constexpr hadpresets::RbcUkqcd::EnsembleParameters testPar{0.005, 0.0362, 1., 2.,
8, 8, 4};
// C0 // C0
inline static constexpr EnsembleParameters c0UnitaryPar{0.00078, 0.0362, 1.8, 2., 48, 96, 24}; 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 unsigned int c0ZMobiusLs = 10;
@@ -79,6 +82,11 @@ struct RbcUkqcd
inline static constexpr DeflationParameters c1m20DeflPar{3.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}; 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 // Light solvers: load deflation from disk
static inline void addC0LightZMobiusLCDSolver(Application &app, const std::string solverName, static inline void addC0LightZMobiusLCDSolver(Application &app, const std::string solverName,
const std::string gaugeName, const std::string gaugeName,
@@ -87,6 +95,12 @@ struct RbcUkqcd
const std::string boundary = "1 1 1 -1", const std::string boundary = "1 1 1 -1",
const double residual = 1.0e-8); 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 static inline void
addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName, addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName,
const std::string gaugeTransform, const std::string eigenpackPath, const std::string gaugeTransform, const std::string eigenpackPath,
@@ -125,6 +139,10 @@ struct RbcUkqcd
const std::string gaugeName, const std::string gaugeName,
const std::string boundary = "1 1 1 -1", const std::string boundary = "1 1 1 -1",
const double residual = 1.0e-8); 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 addM0StrangeSolver(Application &app, const std::string solverName, static inline void addM0StrangeSolver(Application &app, const std::string solverName,
const std::string gaugeName, const std::string gaugeName,
const std::string boundary = "1 1 1 1", const std::string boundary = "1 1 1 1",
@@ -147,6 +165,10 @@ struct RbcUkqcd
const double residual = 1.0e-8); const double residual = 1.0e-8);
// Charm solvers, mass is a free parameter // 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, static inline void addC0CharmSolver(Application &app, const std::string solverName,
const std::string gaugeName, const double mass, const std::string gaugeName, const double mass,
const std::string boundary = "1 1 1 -1", const std::string boundary = "1 1 1 -1",
@@ -159,6 +181,46 @@ struct RbcUkqcd
// Implementations ///////////////////////////////////////////////////////////////////////////////// // Implementations /////////////////////////////////////////////////////////////////////////////////
// Light
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;
// 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.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) // Light C0 (load deflation from disk)
void RbcUkqcd::addC0LightZMobiusLCDSolver(Application &app, const std::string solverName, void RbcUkqcd::addC0LightZMobiusLCDSolver(Application &app, const std::string solverName,
const std::string gaugeName, const std::string gaugeName,
@@ -236,6 +298,139 @@ void RbcUkqcd::addC0LightZMobiusLCDSolver(Application &app, const std::string so
app.createModule<MSolver::ZMixedPrecisionRBPrecCGBatched>(solverName, solverPar); 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 M0 (load deflation from disk) // Light M0 (load deflation from disk)
void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName, void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName,
const std::string gaugeName, const std::string gaugeTransform, const std::string gaugeName, const std::string gaugeTransform,
@@ -497,6 +692,51 @@ void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverN
RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, boundary, 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 // Charm C0
void RbcUkqcd::addC0CharmSolver(Application &app, const std::string solverName, void RbcUkqcd::addC0CharmSolver(Application &app, const std::string solverName,
const std::string gaugeName, const double mass, const std::string gaugeName, const double mass,