From 7d8d2503893369c6ceb31c23c22e4b943725106b Mon Sep 17 00:00:00 2001 From: Quadro Date: Wed, 9 Jun 2021 16:30:39 -0400 Subject: [PATCH] Complete ? --- HMC/Mobius2f_DDHMC.cc | 208 ++++++++++++++++++++++++++---------------- 1 file changed, 129 insertions(+), 79 deletions(-) diff --git a/HMC/Mobius2f_DDHMC.cc b/HMC/Mobius2f_DDHMC.cc index c44f3ef1..d3648813 100644 --- a/HMC/Mobius2f_DDHMC.cc +++ b/HMC/Mobius2f_DDHMC.cc @@ -36,22 +36,16 @@ directory #include #include -//#include -//#include - +#include NAMESPACE_BEGIN(Grid); - /* - * Need a plan for gauge field update for mixed precision in HMC (2x speed up) - * -- Store the single prec action operator. - * -- Clone the gauge field from the operator function argument. - * -- Build the mixed precision operator dynamically from the passed operator and single prec clone. - */ +#define MIXED_PRECISION NAMESPACE_END(Grid); -int main(int argc, char **argv) { +int main(int argc, char **argv) +{ using namespace Grid; Grid_init(&argc, &argv); @@ -60,14 +54,24 @@ int main(int argc, char **argv) { std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; // Typedefs to simplify notation - typedef WilsonImplR FermionImplPolicy; - typedef MobiusFermionR FermionAction; - typedef typename FermionAction::Impl_t Fimpl; - typedef DirichletFermionOperator DirichletFermion; + typedef WilsonImplR FimplD; + typedef WilsonImplF FimplF; + typedef FermionOperator FermionOperatorF; + typedef FermionOperator FermionOperatorD; + typedef MobiusFermionR FermionActionD; + typedef MobiusFermionF FermionActionF; + typedef DirichletFermionOperator DirichletFermionD; + typedef DirichletFermionOperator DirichletFermionF; typedef MobiusEOFAFermionR FermionEOFAAction; - typedef typename FermionAction::FermionField FermionField; - + typedef typename FermionActionD::FermionField FermionFieldD; + typedef typename FermionActionF::FermionField FermionFieldF; + + typedef SchurDiagMooeeOperator,FermionFieldF> LinearOperatorF; + typedef SchurDiagMooeeOperator,FermionFieldD> LinearOperatorD; + typedef SchurDiagMooeeDagOperator,FermionFieldF> LinearOperatorDagF; + typedef SchurDiagMooeeDagOperator,FermionFieldD> LinearOperatorDagD; + typedef Grid::XmlReader Serialiser; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -78,7 +82,7 @@ int main(int argc, char **argv) { // MD.name = std::string("Force Gradient"); typedef GenericHMCRunner HMCWrapper; MD.name = std::string("MinimumNorm2"); - MD.MDsteps = 12; + MD.MDsteps = 6; MD.trajL = 1.0; HMCparameters HMCparams; @@ -107,10 +111,9 @@ int main(int argc, char **argv) { TheHMC.Resources.SetRNGSeeds(RNGpar); // Momentum Dirichlet - Coordinate Block({16,16,16,16}); + Coordinate Block({0,0,0,16}); // Temporal only decomposition - // TheHMC.Resources.SetMomentumFilter(new DirichletFilter(Block)); - TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block,1)); + TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block)); // Construct observables // here there is too much indirection typedef PlaquetteMod PlaqObs; @@ -120,14 +123,13 @@ int main(int argc, char **argv) { const int Ls = 16; Real beta = 2.13; Real light_mass = 0.01; - Real strange_mass = 0.04; Real pv_mass = 1.0; RealD M5 = 1.8; RealD b = 1.0; RealD c = 0.0; - std::vector hasenbusch({ 0.04, 0.2 }); - + std::vector hasenbusch({ 0.15, 0.4, 0.7 }); + auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); @@ -146,25 +148,28 @@ int main(int argc, char **argv) { // temporarily need a gauge field LatticeGaugeField U(GridPtr); + LatticeGaugeFieldF UF(GridPtrF); // These lines are unecessary if BC are all periodic std::vector boundary = {1,1,1,-1}; - FermionAction::ImplParams Params(boundary); - FermionAction::ImplParams DirichletParams(boundary); + FermionActionD::ImplParams Params(boundary); + FermionActionD::ImplParams DirichletParams(boundary); DirichletParams.locally_periodic=true; double ActionStoppingCondition = 1e-10; - double DerivativeStoppingCondition = 1e-10; + double DerivativeStoppingCondition = 1e-8; + double BoundaryDerivativeStoppingCondition = 1e-6; double MaxCGIterations = 30000; //////////////////////////////////// // Collect actions //////////////////////////////////// ActionLevel Level1(1); - ActionLevel Level2(8); + ActionLevel Level2(1); + ActionLevel Level3(8); - ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); - ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); + ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); + ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); //////////////////////////////////// // up down action @@ -186,78 +191,123 @@ int main(int argc, char **argv) { // into single op for any operator pair. // Same issue prevents using MxPCG in the Heatbath step ////////////////////////////////////////////////////////////// - std::vector PeriNumerators; - std::vector PeriDenominators; - std::vector DNumerators; - std::vector DDenominators; - std::vector DirichletNumerators; - std::vector DirichletDenominators; - std::vector *> BoundaryNumerators; - std::vector *> BoundaryDenominators; + ///////////////////////////////////////////////// + // These are consumed/owned by the Dirichlet wrappers + ///////////////////////////////////////////////// + std::vector DNumeratorsD; + std::vector DNumeratorsF; + std::vector DDenominatorsD; + std::vector DDenominatorsF; + ///////////////////////////////////////////////// + // Dirichlet wrappers + ///////////////////////////////////////////////// + std::vector DirichletNumeratorsD; + std::vector DirichletNumeratorsF; + std::vector DirichletDenominatorsD; + std::vector DirichletDenominatorsF; + + std::vector *> Quotients; + std::vector *> BoundaryQuotients; + + typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; + std::vector ActionMPCG; + std::vector MPCG; + std::vector LinOpD; + std::vector LinOpF; + + const int MX_inner = 1000; + const RealD MX_tol = 1.0e-6; - std::vector *> Quotients; - std::vector *> BoundaryQuotients; - std::vector *> BoundaryNums; - std::vector *> BoundaryDens; for(int h=0;h - (*PeriNumerators[h], - *DirichletNumerators[h], - ActionCG,ActionCG, - ActionCG,ActionCG, - Block)); + DirichletNumeratorsD.push_back (new DirichletFermionD(*DNumeratorsD[h],Block)); + DirichletNumeratorsF.push_back (new DirichletFermionF(*DNumeratorsF[h],Block)); + DirichletDenominatorsD.push_back(new DirichletFermionD(*DDenominatorsD[h],Block)); + DirichletDenominatorsF.push_back(new DirichletFermionF(*DDenominatorsF[h],Block)); - BoundaryDenominators.push_back (new SchurFactoredFermionOperator - (*PeriDenominators[h], - *DirichletDenominators[h], - ActionCG,ActionCG, - ActionCG,ActionCG, - Block)); + // Dirichlet Schur even odd MpsDagMpc operators on local domains + LinOpD.push_back(new LinearOperatorD(*DirichletDenominatorsD[h])); + LinOpF.push_back(new LinearOperatorF(*DirichletDenominatorsF[h])); + + // Derivative + MPCG.push_back(new MxPCG(DerivativeStoppingCondition, MX_tol, + MX_inner, + MaxCGIterations, + FrbGridF, + *DirichletDenominatorsF[h],*DirichletDenominatorsD[h], + *LinOpF[h], *LinOpD[h]) ); + + // Action + ActionMPCG.push_back(new MxPCG(ActionStoppingCondition, MX_tol, + MX_inner, + MaxCGIterations, + FrbGridF, + *DirichletDenominatorsF[h],*DirichletDenominatorsD[h], + *LinOpF[h], *LinOpD[h]) ); //////////////////////////////////////////////////////////////////////////// // Standard CG for 2f force //////////////////////////////////////////////////////////////////////////// - std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl; - Quotients.push_back (new - TwoFlavourEvenOddRatioPseudoFermionAction - (*DirichletNumerators[h], - *DirichletDenominators[h], - ActionCG, + TwoFlavourEvenOddRatioPseudoFermionAction + (*DirichletNumeratorsD[h], + *DirichletDenominatorsD[h], + *ActionMPCG[h], + *ActionMPCG[h], ActionCG)); - std::cout << GridLogMessage << " 2f Boundary quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl; - - BoundaryQuotients.push_back(new - DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion - (*BoundaryNumerators[h], - *BoundaryDenominators[h])); - + Level2.push_back(Quotients[h]); } - for(int h=0;h BoundaryNumerator(PeriNumeratorD,PeriNumeratorF, + DirichletNumeratorD,DirichletNumeratorF, + Block); + + SchurFactoredFermionOperator BoundaryDenominator(PeriDenominatorD,PeriDenominatorF, + DirichletDenominatorD,DirichletDenominatorF, + Block); + Level1.push_back(new + DomainDecomposedBoundaryTwoFlavourRatioPseudoFermion + (BoundaryNumerator, + BoundaryDenominator, + BoundaryDerivativeStoppingCondition,ActionStoppingCondition,MX_tol)); ///////////////////////////////////////////////////////////// // Gauge action ///////////////////////////////////////////////////////////// - Level2.push_back(&GaugeAction); + Level3.push_back(&GaugeAction); TheHMC.TheAction.push_back(Level1); TheHMC.TheAction.push_back(Level2); + TheHMC.TheAction.push_back(Level3); std::cout << GridLogMessage << " Action complete "<< std::endl; /////////////////////////////////////////////////////////////