From 21165ed48961bca438136980ecceb86724be8ef1 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:08:23 +0200 Subject: [PATCH] Better logging and moving momentumfilter --- Grid/qcd/hmc/HMC.h | 80 ++++++++++++++----- Grid/qcd/hmc/integrators/Integrator.h | 1 - Grid/qcd/hmc/integrators/MomentumFilter.h | 94 ----------------------- 3 files changed, 59 insertions(+), 116 deletions(-) delete mode 100644 Grid/qcd/hmc/integrators/MomentumFilter.h diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index b05e0561..2e391caa 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -116,22 +116,17 @@ private: random(sRNG, rn_test); - std::cout << GridLogMessage - << "--------------------------------------------------\n"; - std::cout << GridLogMessage << "exp(-dH) = " << prob - << " Random = " << rn_test << "\n"; - std::cout << GridLogMessage - << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "exp(-dH) = " << prob << " Random = " << rn_test << "\n"; + std::cout << GridLogMessage << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n"; if ((prob > 1.0) || (rn_test <= prob)) { // accepted std::cout << GridLogMessage << "Metropolis_test -- ACCEPTED\n"; - std::cout << GridLogMessage - << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; return true; } else { // rejected std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n"; - std::cout << GridLogMessage - << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; return false; } } @@ -140,18 +135,63 @@ private: // Evolution ///////////////////////////////////////////////////////// RealD evolve_hmc_step(Field &U) { - TheIntegrator.refresh(U, sRNG, pRNG); // set U and initialize P and phi's - RealD H0 = TheIntegrator.S(U); // initial state action + GridBase *Grid = U.Grid(); + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // Mainly for DDHMC perform a random translation of U modulo volume + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Random shifting gauge field by ["; + for(int d=0;dNd();d++) { + + int L = Grid->GlobalDimensions()[d]; + + RealD rn_uniform; random(sRNG, rn_uniform); + + int shift = (int) (rn_uniform*L); + + std::cout << shift; + if(dNd()-1) std::cout <<","; + else std::cout <<"]\n"; + + U = Cshift(U,d,shift); + } + std::cout << GridLogMessage << "--------------------------------------------------\n"; + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // set U and initialize P and phi's + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Refresh momenta and pseudofermions"; + TheIntegrator.refresh(U, sRNG, pRNG); + std::cout << GridLogMessage << "--------------------------------------------------\n"; + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // initial state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Compute initial action"; + RealD H0 = TheIntegrator.S(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; std::streamsize current_precision = std::cout.precision(); std::cout.precision(15); std::cout << GridLogMessage << "Total H before trajectory = " << H0 << "\n"; std::cout.precision(current_precision); + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << " Molecular Dynamics evolution "; TheIntegrator.integrate(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; - RealD H1 = TheIntegrator.S(U); // updated state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // updated state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Compute final action"; + RealD H1 = TheIntegrator.S(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; /////////////////////////////////////////////////////////// if(0){ @@ -164,17 +204,14 @@ private: } /////////////////////////////////////////////////////////// - std::cout.precision(15); - std::cout << GridLogMessage << "Total H after trajectory = " << H1 - << " dH = " << H1 - H0 << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Total H after trajectory = " << H1 << " dH = " << H1 - H0 << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; std::cout.precision(current_precision); return (H1 - H0); } - - - public: ///////////////////////////////////////// @@ -196,8 +233,11 @@ public: // Actual updates (evolve a copy Ucopy then copy back eventually) unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory; + for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) { + std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n"; + if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) { std::cout << GridLogMessage << "-- Thermalization" << std::endl; } @@ -217,8 +257,6 @@ public: if (accept) Ucur = Ucopy; - - double t1=usecond(); std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl; diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index aa28c6c8..6b0e3caf 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -33,7 +33,6 @@ directory #define INTEGRATOR_INCLUDED #include -#include "MomentumFilter.h" NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/hmc/integrators/MomentumFilter.h b/Grid/qcd/hmc/integrators/MomentumFilter.h deleted file mode 100644 index 2a15d80c..00000000 --- a/Grid/qcd/hmc/integrators/MomentumFilter.h +++ /dev/null @@ -1,94 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h - -Copyright (C) 2015 - -Author: Christopher Kelly -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 */ -//-------------------------------------------------------------------- -#ifndef MOMENTUM_FILTER -#define MOMENTUM_FILTER - -NAMESPACE_BEGIN(Grid); - -//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh - -template -struct MomentumFilterBase{ - virtual void applyFilter(MomentaField &P) const; -}; - -//Do nothing -template -struct MomentumFilterNone: public MomentumFilterBase{ - void applyFilter(MomentaField &P) const override{} -}; - -//Multiply each site/direction by a Lorentz vector complex number field -//Can be used to implement a mask, zeroing out sites -template -struct MomentumFilterApplyPhase: public MomentumFilterBase{ - typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type - typedef typename MomentaField::scalar_type scalar_type; //scalar complex type - typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction - typedef Lattice LatticeLorentzScalarType; - - LatticeLorentzScalarType phase; - - MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){} - - //Default to uniform field of (1,0) - MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){ - LorentzScalarType one; - for(int mu=0;mu