/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid Source file: Copyright (C) 2015-2016 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 #include #include #include #include #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. */ NAMESPACE_END(Grid); int main(int argc, char **argv) { using namespace Grid; Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); // here make a routine to print all the relevant information on the run 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 MobiusEOFAFermionR FermionEOFAAction; typedef typename FermionAction::FermionField FermionField; typedef Grid::XmlReader Serialiser; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: IntegratorParameters MD; // typedef GenericHMCRunner HMCWrapper; // MD.name = std::string("Leap Frog"); // typedef GenericHMCRunner HMCWrapper; // MD.name = std::string("Force Gradient"); typedef GenericHMCRunner HMCWrapper; MD.name = std::string("MinimumNorm2"); MD.MDsteps = 12; MD.trajL = 1.0; HMCparameters HMCparams; HMCparams.StartTrajectory = 26; HMCparams.Trajectories = 1; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; // HMCparams.StartingType =std::string("ColdStart"); HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); // Grid from the command line arguments --grid and --mpi TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition CheckpointerParameters CPparams; CPparams.config_prefix = "ckpoint_EOFA4D_lat"; CPparams.rng_prefix = "ckpoint_EOFA4D_rng"; CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); RNGModuleParameters RNGpar; RNGpar.serial_seeds = "1 2 3 4 5"; RNGpar.parallel_seeds = "6 7 8 9 10"; TheHMC.Resources.SetRNGSeeds(RNGpar); // Momentum Dirichlet Coordinate Block({16,16,16,16}); // TheHMC.Resources.SetMomentumFilter(new DirichletFilter(Block)); TheHMC.Resources.SetMomentumFilter(new DDHMCFilter(Block,1)); // Construct observables // here there is too much indirection typedef PlaquetteMod PlaqObs; TheHMC.Resources.AddObservable(); ////////////////////////////////////////////// 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 }); auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); Coordinate latt = GridDefaultLatt(); Coordinate mpi = GridDefaultMpi(); Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd()); Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd()); auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi); auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF); auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF); auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF); IwasakiGaugeActionR GaugeAction(beta); // temporarily need a gauge field LatticeGaugeField U(GridPtr); // These lines are unecessary if BC are all periodic std::vector boundary = {1,1,1,-1}; FermionAction::ImplParams Params(boundary); FermionAction::ImplParams DirichletParams(boundary); DirichletParams.locally_periodic=true; double ActionStoppingCondition = 1e-10; double DerivativeStoppingCondition = 1e-10; double MaxCGIterations = 30000; //////////////////////////////////// // Collect actions //////////////////////////////////// ActionLevel Level1(1); ActionLevel Level2(8); ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); //////////////////////////////////// // up down action //////////////////////////////////// std::vector light_den; std::vector light_num; int n_hasenbusch = hasenbusch.size(); light_den.push_back(light_mass); for(int h=0;h PeriNumerators; std::vector PeriDenominators; std::vector DNumerators; std::vector DDenominators; std::vector DirichletNumerators; std::vector DirichletDenominators; std::vector *> BoundaryNumerators; std::vector *> BoundaryDenominators; std::vector *> Quotients; std::vector *> BoundaryQuotients; std::vector *> BoundaryNums; std::vector *> BoundaryDens; for(int h=0;h (*PeriNumerators[h], *DirichletNumerators[h], ActionCG,Block)); BoundaryDenominators.push_back (new SchurFactoredFermionOperator (*PeriDenominators[h], *DirichletDenominators[h], ActionCG,Block)); //////////////////////////////////////////////////////////////////////////// // 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, 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], ActionCG,ActionCG)); } for(int h=0;h