HadronsPresets/RbcUkqcd.hpp
2024-03-20 14:04:27 +09:00

577 lines
25 KiB
C++

/*
* Copyright (C) 2023
*
* Author: Antonin Portelli <antonin.portelli@me.com>
* Elements based on production templates from Raoul Hodgson
*
* 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>
namespace hadpresets
{
using namespace Grid::Hadrons;
using namespace std::complex_literals;
struct RbcUkqcd
{
// Ensemble parameters
struct EnsembleParameters
{
double ml, ms, M5, scale;
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<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};
// 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: 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,
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 addM0StrangeSolver(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 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 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 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 /////////////////////////////////////////////////////////////////////////////////
// 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);
}
// 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 = 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::LoadCoarseFermionEigenPack250F::Par epPar;
epPar.filestem = eigenpackPath;
epPar.multiFile = true;
epPar.redBlack = true;
epPar.sizeFine = 250;
epPar.sizeCoarse = 2000;
epPar.Ls = 12;
epPar.blockSize = "4 4 4 4 12";
epPar.orthogonalise = false;
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::MixedPrecisionRBPrecCGBatched::Par solverPar;
solverPar.innerAction = prefix + "_dwf_fp32";
solverPar.outerAction = prefix + "_dwf";
solverPar.maxInnerIteration = 400;
solverPar.maxOuterIteration = 100;
solverPar.maxPatchupIteration = 1000;
solverPar.residual = residual;
solverPar.updateResidual = true;
solverPar.innerGuesser = prefix + "_iguesser";
solverPar.outerGuesser = "";
app.createModule<MSolver::MixedPrecisionRBPrecCGBatched>(solverName, solverPar);
}
// Deflation at runtime
void RbcUkqcd::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)
{
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);
// 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";
solverPar.outerAction = prefix + "_dwf";
solverPar.maxInnerIteration = 30000;
solverPar.maxOuterIteration = 100;
solverPar.residual = residual;
solverPar.innerGuesser = "";
solverPar.outerGuesser = "";
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::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 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 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