diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc b/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc index 059a0f20..5572d11f 100644 --- a/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc +++ b/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc @@ -164,11 +164,6 @@ int main(int argc, char **argv) { typedef MobiusEOFAFermionF FermionEOFAActionF; typedef typename FermionActionF::FermionField FermionFieldF; - typedef WilsonImplD2 FermionImplPolicyD2; - typedef MobiusFermionD2 FermionActionD2; - typedef MobiusEOFAFermionD2 FermionEOFAActionD2; - typedef typename FermionActionD2::FermionField FermionFieldD2; - typedef Grid::XmlReader Serialiser; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -268,10 +263,8 @@ int main(int argc, char **argv) { typedef SchurDiagMooeeOperator LinearOperatorF; typedef SchurDiagMooeeOperator LinearOperatorD; - typedef SchurDiagMooeeOperator LinearOperatorD2; typedef SchurDiagMooeeOperator LinearOperatorEOFAF; typedef SchurDiagMooeeOperator LinearOperatorEOFAD; - typedef SchurDiagMooeeOperator LinearOperatorEOFAD2; typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG_EOFA; @@ -324,7 +317,6 @@ int main(int argc, char **argv) { // temporarily need a gauge field LatticeGaugeFieldD U(GridPtr); U=Zero(); LatticeGaugeFieldF UF(GridPtrF); UF=Zero(); - LatticeGaugeFieldD2 UD2(GridPtrF); UD2=Zero(); std::cout << GridLogMessage << " Running the HMC "<< std::endl; TheHMC.ReadCommandLine(argc,argv); // params on CML or from param file @@ -428,7 +420,7 @@ int main(int argc, char **argv) { ActionCGL, ActionCGR, DerivativeCGL, DerivativeCGR, SFRp, true); - // Level2.push_back(&EOFA); + Level2.push_back(&EOFA); //////////////////////////////////// // up down action @@ -453,15 +445,13 @@ int main(int argc, char **argv) { std::vector Denominators; std::vector NumeratorsF; std::vector DenominatorsF; - std::vector NumeratorsD2; - std::vector DenominatorsD2; std::vector *> Quotients; std::vector ActionMPCG; std::vector MPCG; #define MIXED_PRECISION #ifdef MIXED_PRECISION - std::vector *> Bdys; + std::vector *> Bdys; #else std::vector *> Bdys; #endif @@ -536,27 +526,15 @@ int main(int argc, char **argv) { Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],CG)); } else { #ifdef MIXED_PRECISION - // Use the D2 data types and make them use same grid as single - FermionActionD2::ImplParams ParamsDenD2(boundary); - FermionActionD2::ImplParams ParamsNumD2(boundary); - - ParamsDenD2.dirichlet = ParamsDen.dirichlet; - ParamsDenD2.partialDirichlet = ParamsDen.partialDirichlet; - DenominatorsD2.push_back(new FermionActionD2(UD2,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_den[h],M5,b,c, ParamsDenD2)); - - ParamsNumD2.dirichlet = ParamsNum.dirichlet; - ParamsNumD2.partialDirichlet = ParamsNum.partialDirichlet; - NumeratorsD2.push_back (new FermionActionD2(UD2,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_num[h],M5,b,c, ParamsNumD2)); - - Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction( + Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction( *Numerators[h],*Denominators[h], *NumeratorsF[h],*DenominatorsF[h], - *NumeratorsD2[h],*DenominatorsD2[h], + *Numerators[h],*Denominators[h], OFRp, SP_iters) ); - Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction( + Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction( *Numerators[h],*Denominators[h], *NumeratorsF[h],*DenominatorsF[h], - *NumeratorsD2[h],*DenominatorsD2[h], + *Numerators[h],*Denominators[h], OFRp, SP_iters) ); #else Bdys.push_back( new GeneralEvenOddRatioRationalPseudoFermionAction(*Numerators[h],*Denominators[h],OFRp)); diff --git a/TODO b/TODO index a6d0f2ac..750deb55 100644 --- a/TODO +++ b/TODO @@ -2,6 +2,8 @@ - - Also faster non-atomic reduction - - Remaining PRs - - DDHMC + - - MixedPrec is the action eval, high precision + - - MixedPrecCleanup is the force eval, low precision ================= ================= diff --git a/tests/forces/Test_bdy.cc b/tests/forces/Test_bdy.cc new file mode 100644 index 00000000..c2c97d0d --- /dev/null +++ b/tests/forces/Test_bdy.cc @@ -0,0 +1,305 @@ +/* + + 2f Full det MdagM 10^6 force ~ 1.3e7 +rid : Message : 1767.283471 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Grid : Message : 1767.283476 s : S1 : 1.52885e+09 +Grid : Message : 1767.283480 s : S2 : 1.52886e+09 +Grid : Message : 1767.283482 s : dS : 8877.34 +Grid : Message : 1767.283483 s : dSpred : 8877.7 +Grid : Message : 1767.283484 s : diff : -0.360484 +Grid : Message : 1767.283485 s : ********************************************************* + + 2f Full det MpcdagMpc 10^6 force ~ 1.8e6 +Grid : Message : 2399.576962 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Grid : Message : 2399.576968 s : S1 : 1.52885e+09 +Grid : Message : 2399.576972 s : S2 : 1.52886e+09 +Grid : Message : 2399.576974 s : dS : 9728.49 +Grid : Message : 2399.576975 s : dSpred : 9726.58 +Grid : Message : 2399.576976 s : diff : 1.90683 +Grid : Message : 2399.576977 s : ********************************************************* + + 2f bdy MdagM 1500 force Force ~ 2800 +Grid : Message : 4622.385061 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Grid : Message : 4622.385067 s : S1 : 1.52885e+09 +Grid : Message : 4622.385071 s : S2 : 1.52885e+09 +Grid : Message : 4622.385072 s : dS : 25.4944 +Grid : Message : 4622.385073 s : dSpred : 25.4672 +Grid : Message : 4622.385074 s : diff : 0.0271414 +Grid : Message : 4622.385075 s : ********************************************************* + + 2f bdy MpcdagMpc 10^6 force Force ~ 2200 +Grid : Message : 4622.385061 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Grid : Message : 4622.385067 s : S1 : 1.52885e+09 +Grid : Message : 4622.385071 s : S2 : 1.52885e+09 +Grid : Message : 4622.385072 s : dS : 25.4944 +Grid : Message : 4622.385073 s : dSpred : 25.4672 +Grid : Message : 4622.385074 s : diff : 0.0271414 +Grid : Message : 4622.385075 s : ********************************************************* + + 1f Bdy Det + Optimisation log: looser rational AND MD tolerances sloppy +MobiusForce.221179 -- same as HMC. dS is mispredicted Forece ~2.8 +Grid : Message : 6582.258991 s : dS : 0.024478 +Grid : Message : 6582.258992 s : dSpred : 0.00791876 +Grid : Message : 6582.258994 s : diff : 0.0165592 + +MobiusForce.221193 -- tight rational AND MD tolerances to 1e-8 ~ 2.8 same +Grid : Message : 1964.939209 s : S1 : 7.64404e+08 +Grid : Message : 1964.939213 s : S2 : 7.64404e+08 +Grid : Message : 1964.939215 s : dS : -0.00775838 <--- too loose even on action +Grid : Message : 1964.939216 s : dSpred : -0.00416793 +Grid : Message : 1964.939217 s : diff : -0.00359045 + +MobiusForce.221394 -- looser rational, MD tol 1e-8 ~ 2.8 same +Grid : Message : 1198.346720 s : S1 : 764404649.48886 +Grid : Message : 1198.346760 s : S2 : 764404649.5133 +Grid : Message : 1198.346780 s : dS : 0.024440884590149 +Grid : Message : 1198.346800 s : dSpred : 0.0079145154465184 +Grid : Message : 1198.346810 s : diff : 0.016526369143631 + +MobiusForce.221394 -- tight rational, MD tol sloppy Force ~ 2.8 +Grid : Message : 2376.921950 s : S1 : 764404436.44069 +Grid : Message : 2376.921954 s : S2 : 764404436.43299 +Grid : Message : 2376.921956 s : dS : -0.0076971054077148 +Grid : Message : 2376.921958 s : dSpred : -0.0041610472282526 +Grid : Message : 2376.921959 s : diff : -0.0035360581794623 + +*/ + +// +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_double_ratio.cc + + Copyright (C) 2022 + +Author: Peter Boyle + + This program 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. + + This program 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 this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; + +typedef MobiusFermionD FermionAction; +typedef WilsonImplD FimplD; +typedef WilsonImplD FermionImplPolicy; + +template +void ForceTest(Action &action,LatticeGaugeField & U,MomentumFilterBase &Filter) +{ + GridBase *UGrid = U.Grid(); + + std::vector seeds({1,2,3,5}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + + LatticeColourMatrix Pmu(UGrid); + LatticeGaugeField P(UGrid); + LatticeGaugeField UdSdU(UGrid); + + std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); + Pmu= PeekIndex(P,mu); + dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; + } + ComplexD dSpred = sum(dS); + RealD diff = S2-S1-dSpred.real(); + + std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<1 ? 1 : 0; + + Coordinate NonDirichlet(Nd+1,0); + Coordinate Dirichlet(Nd+1,0); + Dirichlet[1] = CommDim[0]*latt_size[0]/mpi_layout[0] * shm[0]; + Dirichlet[2] = CommDim[1]*latt_size[1]/mpi_layout[1] * shm[1]; + Dirichlet[3] = CommDim[2]*latt_size[2]/mpi_layout[2] * shm[2]; + Dirichlet[4] = CommDim[3]*latt_size[3]/mpi_layout[3] * shm[3]; + + Coordinate Block4(Nd); + Block4[0] = Dirichlet[1]; + Block4[1] = Dirichlet[2]; + Block4[2] = Dirichlet[3]; + Block4[3] = Dirichlet[4]; + + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + FermionAction::ImplParams ParamsDir(boundary); + Params.dirichlet=NonDirichlet; + ParamsDir.dirichlet=Dirichlet; + ParamsDir.partialDirichlet=1; + + ///////////////////// Gauge Field and Gauge Forces //////////////////////////// + LatticeGaugeField U(UGrid); + + RealD beta=6.0; + WilsonGaugeActionR PlaqAction(beta); + IwasakiGaugeActionR RectAction(beta); + + MomentumFilterNone FilterNone; + ForceTest(PlaqAction,U,FilterNone); + ForceTest(RectAction,U,FilterNone); + + //////////////////////////////////// + // Action + //////////////////////////////////// + RealD mass=0.00078; + RealD pvmass=1.0; + RealD M5=1.8; + RealD b=1.5; + RealD c=0.5; + + // Double versions + FermionAction DdwfPeriodic(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,Params); + FermionAction PVPeriodic (U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pvmass,M5,b,c,Params); + FermionAction DdwfDirichlet(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,ParamsDir); + + double StoppingCondition = 1.0e-8; + double MaxCGIterations = 50000; + ConjugateGradient CG(StoppingCondition,MaxCGIterations); + + //////////////////// Two Flavour Determinant Ratio /////////////////////////////// + TwoFlavourRatioPseudoFermionAction Nf2(PVPeriodic, DdwfPeriodic,CG,CG); + // ForceTest(Nf2,U,FilterNone); + + //////////////////// Two Flavour Determinant force test Even Odd /////////////////////////////// + TwoFlavourEvenOddRatioPseudoFermionAction Nf2eo(PVPeriodic, DdwfPeriodic,CG,CG); + // ForceTest(Nf2eo,U,FilterNone); + + //////////////////// Domain forces //////////////////// + int Width=4; + DDHMCFilter DDHMCFilter(Block4,Width); + + //////////////////// Two flavour boundary det //////////////////// + TwoFlavourRatioPseudoFermionAction BdyNf2(DdwfDirichlet, DdwfPeriodic,CG,CG); + // ForceTest(BdyNf2,U,DDHMCFilter); + + //////////////////// Two flavour eo boundary det //////////////////// + TwoFlavourEvenOddRatioPseudoFermionAction BdyNf2eo(DdwfDirichlet, DdwfPeriodic,CG,CG); + // ForceTest(BdyNf2eo,U,DDHMCFilter); + + //////////////////// One flavour boundary det //////////////////// + OneFlavourRationalParams OFRp; // Up/down + OFRp.lo = 4.0e-5; + OFRp.hi = 90.0; + OFRp.MaxIter = 60000; + OFRp.tolerance= 1.0e-8; + OFRp.mdtolerance= 1.0e-6; + OFRp.degree = 18; + OFRp.precision= 80; + OFRp.BoundsCheckFreq=0; + std::vector ActionTolByPole({ + 1.0e-8,1.0e-8,1.0e-8,1.0e-8, + 1.0e-8,1.0e-8,1.0e-8,1.0e-8, + 1.0e-8,1.0e-8,1.0e-8,1.0e-8 + }); + std::vector MDTolByPole({ + 1.0e-6,3.0e-7,1.0e-7,1.0e-7, // Orig sloppy + // 1.0e-8,1.0e-8,1.0e-8,1.0e-8, + 1.0e-8,1.0e-8,1.0e-8,1.0e-8, + 1.0e-8,1.0e-8,1.0e-8,1.0e-8 + }); + OneFlavourEvenOddRatioRationalPseudoFermionAction BdySqrt(DdwfDirichlet,DdwfPeriodic,OFRp); + ForceTest(BdySqrt,U,DDHMCFilter); + + Grid_finalize(); +}