Domain wall even-odd 2f HMC with wilson gauge and PV 2f ratio now running and giving small dH.

Azusa is working hard on the rectangle term and we'll hopefully start reproducing plaquettes
from RBC-UKQCD parameters soon !

My new laptop is pretty warm and is starting to groan ;)

TODO

@@ -7,19 +7,19 @@ Done:
   - TwoFlavourEvenOdd        
   - TwoFlavourRatio
   - TwoFlavourRatioEvenOdd
-  => generalise to non-const EE
 
-fill in:
-  - OneFlavourRationalEvenOddConstEE
+Done:
   - OneFlavourRationalEvenOdd
-  - OneFlavourRationalRatioEvenOddConstEE
   - OneFlavourRationalRatioEvenOdd
   - OneFlavourRationalRatio
 
-- Clean up HMC
-- Fix a threading bug that has been introduced and prevents HMC running hybrid OMP mode
+TODO:
+=> generalise to non-const EE
+=> Test DWF HMC
+=> Clean up HMC
+   - Fix a threading bug that has been introduced and prevents HMC running hybrid OMP mode
 
-- Integrators
+=> Integrators
   - Force Gradient
   - Multi-timescale looks broken and operating on single timescale for now.
     Fix/debug/rewrite this 
@@ -31,8 +31,6 @@ fill in:
 
 - Link smearing/boundary conds; Policy class based implementation
 
-- Test DWF HMC
-
 - Rectangle gauge actions.
   Iwasaki,
   Symanzik,
@@ -102,9 +100,6 @@ Insert/Extract
    // localMaxAbs
    // Fourier transform equivalent.]
 
-================================================================
-*** New Functionality
-================================================================
 
 * CovariantShift support                             -----Use a class to store gauge field? (parallel transport?)
 

lib/lattice/Lattice_rng.h

@@ -5,6 +5,37 @@
 
 namespace Grid {
 
+
+  //////////////////////////////////////////////////////////////
+  // Allow the RNG state to be less dense than the fine grid
+  //////////////////////////////////////////////////////////////
+  inline int RNGfillable(GridBase *coarse,GridBase *fine)
+  {
+
+    int rngdims = coarse->_ndimension;
+
+    // trivially extended in higher dims, with locality guaranteeing RNG state is local to node
+    int lowerdims   = fine->_ndimension - coarse->_ndimension;
+    assert(lowerdims >= 0);
+    for(int d=0;d<lowerdims;d++){
+      assert(fine->_simd_layout[d]==1);
+      assert(fine->_processors[d]==1);
+    }
+
+    // local and global volumes subdivide cleanly after SIMDization
+    int multiplicity=1;
+    for(int d=0;d<rngdims;d++){
+      int fd= d+lowerdims;
+      assert(coarse->_processors[d]  == fine->_processors[fd]);
+      assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
+      assert((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[fd]); 
+
+      multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d]; 
+    }
+
+    return multiplicity;
+  }
+
   // Wrap seed_seq to give common interface with random_device
   class fixedSeed {
   public:
@@ -226,7 +257,7 @@ namespace Grid {
       typedef typename vobj::scalar_type scalar_type;
       typedef typename vobj::vector_type vector_type;
       
-      conformable(_grid,l._grid);
+      int multiplicity = RNGfillable(_grid,l._grid);
 
       int     Nsimd =_grid->Nsimd();
       int     osites=_grid->oSites();
@@ -236,17 +267,22 @@ namespace Grid {
 
 PARALLEL_FOR_LOOP
       for(int ss=0;ss<osites;ss++){
-	for(int si=0;si<Nsimd;si++){
 
-	  int gdx = generator_idx(ss,si); // index of generator state
-	  scalar_type *pointer = (scalar_type *)&buf[si];
-	  for(int idx=0;idx<words;idx++){
-	    fillScalar(pointer[idx],dist,_generators[gdx]);
+	for(int m=0;m<multiplicity;m++) {// Draw from same generator multiplicity times
+
+	  int sm=multiplicity*ss+m;      // Maps the generator site to the fine site
+
+	  for(int si=0;si<Nsimd;si++){
+	    int gdx = generator_idx(ss,si); // index of generator state
+	    scalar_type *pointer = (scalar_type *)&buf[si];
+	    for(int idx=0;idx<words;idx++){
+	      fillScalar(pointer[idx],dist,_generators[gdx]);
+	    }
 	  }
 
+	  // merge into SIMD lanes
+	  merge(l._odata[sm],buf);
 	}
-	// merge into SIMD lanes
-	merge(l._odata[ss],buf);
       }
     };
 

lib/lattice/Lattice_transfer.h

@@ -17,6 +17,7 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
   }
 }
 
+
   ////////////////////////////////////////////////////////////////////////////////////////////
   // remove and insert a half checkerboard
   ////////////////////////////////////////////////////////////////////////////////////////////
@@ -298,5 +299,42 @@ inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
 }
 
 
+template<class vobj>
+void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
+{
+  typedef typename vobj::scalar_object sobj;
+
+  GridBase *cg = coarse._grid;
+  GridBase *fg =   fine._grid;
+
+  int nd = cg->_ndimension;
+
+  subdivides(cg,fg); 
+
+  assert(cg->_ndimension==fg->_ndimension);
+
+  std::vector<int> ratio(cg->_ndimension);
+
+  for(int d=0;d<cg->_ndimension;d++){
+    ratio[d] = fg->_fdimensions[d]/cg->_fdimensions[d];
+  }
+
+  std::vector<int> fcoor(nd);
+  std::vector<int> ccoor(nd);
+  for(int g=0;g<fg->gSites();g++){
+
+    fg->GlobalIndexToGlobalCoor(g,fcoor);
+    for(int d=0;d<nd;d++){
+      ccoor[d] = fcoor[d]%cg->_gdimensions[d];
+    }
+    
+    sobj tmp;
+    peekSite(tmp,coarse,ccoor);
+    pokeSite(tmp,fine,fcoor);
+  }
+
+}
+
+
 }
 #endif

scripts/hmc.sh

@@ -3,4 +3,4 @@
 LOG=$1
 SWEEPS=$2
 grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10} END { print S/NR} ' 
-grep dH log_cb | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '
+grep dH $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_gpwilson_even_odd Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -98,6 +98,10 @@ Test_gpwilson_even_odd_SOURCES=Test_gpwilson_even_odd.cc
 Test_gpwilson_even_odd_LDADD=-lGrid
 
 
+Test_hmc_EODWFRatio_SOURCES=Test_hmc_EODWFRatio.cc
+Test_hmc_EODWFRatio_LDADD=-lGrid
+
+
 Test_hmc_EOWilsonFermionGauge_SOURCES=Test_hmc_EOWilsonFermionGauge.cc
 Test_hmc_EOWilsonFermionGauge_LDADD=-lGrid
 

tests/Test_gparity.cc

@@ -6,41 +6,6 @@ using namespace Grid::QCD;
 
 typedef typename GparityDomainWallFermionR::FermionField FermionField;
 
-template<class vobj>
-void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
-{
-  typedef typename vobj::scalar_object sobj;
-
-  GridBase *cg = coarse._grid;
-  GridBase *fg =   fine._grid;
-
-  int nd = cg->_ndimension;
-
-  subdivides(cg,fg); 
-
-  assert(cg->_ndimension==fg->_ndimension);
-
-  std::vector<int> ratio(cg->_ndimension);
-
-  for(int d=0;d<cg->_ndimension;d++){
-    ratio[d] = fg->_fdimensions[d]/cg->_fdimensions[d];
-  }
-
-  std::vector<int> fcoor(nd);
-  std::vector<int> ccoor(nd);
-  for(int g=0;g<fg->gSites();g++){
-
-    fg->GlobalIndexToGlobalCoor(g,fcoor);
-    for(int d=0;d<nd;d++){
-      ccoor[d] = fcoor[d]%cg->_gdimensions[d];
-    }
-    
-    sobj tmp;
-    peekSite(tmp,coarse,ccoor);
-    pokeSite(tmp,fine,fcoor);
-  }
-
-}
 
 int main (int argc, char ** argv)
 {

tests/Test_hmc_EODWFRatio.cc

@@ -0,0 +1,60 @@
+#include "Grid.h"
+
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls = 8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  GridParallelRNG  pRNG(UGrid);
+  pRNG.SeedRandomDevice();
+
+  LatticeLorentzColourMatrix     U(UGrid);
+
+  SU3::HotConfiguration(pRNG, U);
+
+  // simplify template declaration? Strip the lorentz from the second template
+  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
+
+  Real mass=0.04;
+  Real pv  =1.0;
+  RealD M5=1.5;
+
+  DomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  DomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
+  
+  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
+  TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
+  
+  //Collect actions
+  ActionLevel Level1;
+  Level1.push_back(&Nf2);
+  Level1.push_back(&Waction);
+  ActionSet FullSet;
+  FullSet.push_back(Level1);
+
+  // Create integrator
+  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
+  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  IntegratorParameters MDpar(12,20,1.0);
+  Integrator<IntegratorAlgorithm> MDynamics(UGrid,MDpar, FullSet);
+
+  // Create HMC
+  HMCparameters HMCpar;
+  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
+
+  // Run it
+  HMC.evolve(U);
+
+}