Hadrons: support for twisted boundary conditions

Untested but doesn't affect answers when twists are all zero. The zero is the default behaviour
for ImplParams.

Grid/communicator/Communicator_mpi3.cc

@@ -50,8 +50,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv)
       assert(0);
   }
 
-  Grid_quiesce_nodes();
-
   // Never clean up as done once.
   MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world);
 
@@ -124,10 +122,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors,
   // split the communicator
   //////////////////////////////////////////////////////////////////////////////////////////////////////
   //  int Nparent = parent._processors ; 
-  //  std::cout << " splitting from communicator "<<parent.communicator <<std::endl;
   int Nparent;
   MPI_Comm_size(parent.communicator,&Nparent);
-  //  std::cout << " Parent size  "<<Nparent <<std::endl;
 
   int childsize=1;
   for(int d=0;d<processors.size();d++) {
@@ -136,8 +132,6 @@ CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors,
   int Nchild = Nparent/childsize;
   assert (childsize * Nchild == Nparent);
 
-  //  std::cout << " child size  "<<childsize <<std::endl;
-
   std::vector<int> ccoor(_ndimension); // coor within subcommunicator
   std::vector<int> scoor(_ndimension); // coor of split within parent
   std::vector<int> ssize(_ndimension); // coor of split within parent

Grid/communicator/SharedMemoryMPI.cc

@@ -413,7 +413,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     assert(((uint64_t)ptr&0x3F)==0);
     close(fd);
     WorldShmCommBufs[r] =ptr;
-    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+    //    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
   }
   _ShmAlloc=1;
   _ShmAllocBytes  = bytes;
@@ -455,7 +455,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     assert(((uint64_t)ptr&0x3F)==0);
     close(fd);
     WorldShmCommBufs[r] =ptr;
-    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+    //    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
   }
   _ShmAlloc=1;
   _ShmAllocBytes  = bytes;
@@ -499,7 +499,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 #endif
       void * ptr =  mmap(NULL,size, PROT_READ | PROT_WRITE, mmap_flag, fd, 0);
       
-      std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< size<< "bytes)"<<std::endl;
+      //      std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< size<< "bytes)"<<std::endl;
       if ( ptr == (void * )MAP_FAILED ) {       
 	perror("failed mmap");     
 	assert(0);    

Grid/lattice/Lattice_transfer.h

@@ -464,9 +464,11 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
   assert(orthog>=0);
 
   for(int d=0;d<nh;d++){
+    if ( d!=orthog ) {
       assert(lg->_processors[d]  == hg->_processors[d]);
       assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
     }
+  }
 
   // the above should guarantee that the operations are local
   parallel_for(int idx=0;idx<lg->lSites();idx++){
@@ -499,9 +501,11 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slic
   assert(orthog>=0);
 
   for(int d=0;d<nh;d++){
+    if ( d!=orthog ) {
       assert(lg->_processors[d]  == hg->_processors[d]);
       assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
     }
+  }
 
   // the above should guarantee that the operations are local
   parallel_for(int idx=0;idx<lg->lSites();idx++){

Grid/qcd/action/ActionParams.h

@@ -44,12 +44,15 @@ namespace QCD {
   
   struct WilsonImplParams {
     bool overlapCommsCompute;
+    std::vector<Real> twist_n_2pi_L;
     std::vector<Complex> boundary_phases;
     WilsonImplParams() : overlapCommsCompute(false) {
       boundary_phases.resize(Nd, 1.0);
+      twist_n_2pi_L.resize(Nd, 0.0);
     };
-    WilsonImplParams(const std::vector<Complex> phi)
+    WilsonImplParams(const std::vector<Complex> phi) : boundary_phases(phi), overlapCommsCompute(false) {
-      : boundary_phases(phi), overlapCommsCompute(false) {}
+      twist_n_2pi_L.resize(Nd, 0.0);
+    }
   };
 
   struct StaggeredImplParams {

Grid/qcd/action/fermion/FermionOperatorImpl.h

@@ -240,16 +240,30 @@ namespace QCD {
       GaugeLinkField tmp(GaugeGrid);
 
       Lattice<iScalar<vInteger> > coor(GaugeGrid);
+      ////////////////////////////////////////////////////
+      // apply any boundary phase or twists
+      ////////////////////////////////////////////////////
       for (int mu = 0; mu < Nd; mu++) {
 
+	////////// boundary phase /////////////
 	auto pha = Params.boundary_phases[mu];
 	scalar_type phase( real(pha),imag(pha) );
 
-        int Lmu = GaugeGrid->GlobalDimensions()[mu] - 1;
+	int L   = GaugeGrid->GlobalDimensions()[mu];
+        int Lmu = L - 1;
 
         LatticeCoordinate(coor, mu);
 
         U = PeekIndex<LorentzIndex>(Umu, mu);
+
+	// apply any twists
+	RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L;
+	if ( theta != 0.0) { 
+	  scalar_type twphase(::cos(theta),::sin(theta));
+	  U = twphase*U;
+	  std::cout << GridLogMessage << " Twist ["<<mu<<"] "<< Params.twist_n_2pi_L[mu]<< " phase"<<phase <<std::endl;
+	}
+
         tmp = where(coor == Lmu, phase * U, U);
         PokeIndex<LorentzIndex>(Uds, tmp, mu);
 

Hadrons/Modules/MAction/DWF.hpp

@@ -49,7 +49,8 @@ public:
                                     unsigned int, Ls,
                                     double      , mass,
                                     double      , M5,
-                                    std::string , boundary);
+                                    std::string , boundary,
+                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -119,8 +120,9 @@ void TDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename DomainWallFermion<FImpl>::ImplParams implParams;
-    typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
 }

Hadrons/Modules/MAction/MobiusDWF.hpp

@@ -49,7 +49,8 @@ public:
                                     double      , M5,
                                     double      , b,
                                     double      , c,
-                                    std::string , boundary);
+                                    std::string , boundary,
+                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -119,8 +120,9 @@ void TMobiusDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename MobiusFermion<FImpl>::ImplParams implParams;
-    typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                      implParams);

Hadrons/Modules/MAction/ScaledDWF.hpp

@@ -48,7 +48,8 @@ public:
                                     double      , mass,
                                     double      , M5,
                                     double      , scale,
-                                    std::string , boundary);
+                                    std::string , boundary,
+                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -118,8 +119,9 @@ void TScaledDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename ScaledShamirFermion<FImpl>::ImplParams implParams;
-    typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().scale,
                      implParams);

Hadrons/Modules/MAction/Wilson.hpp

@@ -47,7 +47,9 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar,
                                     std::string, gauge,
                                     double     , mass,
-                                    std::string, boundary);
+                                    std::string, boundary,
+                                    std::string, string,
+                                    std::string, twist);
 };
 
 template <typename FImpl>
@@ -113,8 +115,9 @@ void TWilson<FImpl>::setup(void)
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename WilsonFermion<FImpl>::ImplParams implParams;
-    typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
 }

Hadrons/Modules/MAction/WilsonClover.hpp

@@ -51,7 +51,8 @@ public:
 				                    double     , csw_r,
 				                    double     , csw_t,
 				                    WilsonAnisotropyCoefficients ,clover_anisotropy,
-                                    std::string, boundary
+                                    std::string, boundary,
+                                    std::string, twist
 				    );
 };
 
@@ -119,8 +120,9 @@ void TWilsonClover<FImpl>::setup(void)
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename WilsonCloverFermion<FImpl>::ImplParams implParams;
-    typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                      gridRb, par().mass, par().csw_r, par().csw_t, 
                      par().clover_anisotropy, implParams); 

Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -50,7 +50,8 @@ public:
                                     double                           , b,
                                     double                           , c,
                                     std::vector<std::complex<double>>, omega,
-                                    std::string                      , boundary);
+                                    std::string                      , boundary,
+                                    std::string                      , twist);
 };
 
 template <typename FImpl>
@@ -127,8 +128,9 @@ void TZMobiusDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename ZMobiusFermion<FImpl>::ImplParams implParams;
-    typename ZMobiusFermion<FImpl>::ImplParams implParams(boundary);
+    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, omega,
                      par().b, par().c, implParams);

configure.ac

@@ -485,6 +485,7 @@ DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg])
 
 ############### Ouput
 cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
+GRID_CXX="$CXX"
 GRID_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
 GRID_LDFLAGS="$AM_LDFLAGS $LDFLAGS"
 GRID_LIBS=$LIBS
@@ -497,6 +498,7 @@ AM_LDFLAGS="-L${cwd}/Grid $AM_LDFLAGS"
 AC_SUBST([AM_CFLAGS])
 AC_SUBST([AM_CXXFLAGS])
 AC_SUBST([AM_LDFLAGS])
+AC_SUBST([GRID_CXX])
 AC_SUBST([GRID_CXXFLAGS])
 AC_SUBST([GRID_LDFLAGS])
 AC_SUBST([GRID_LIBS])

grid-config.in

@@ -61,6 +61,10 @@ while test $# -gt 0; do
       echo @GRID_CXXFLAGS@
     ;;
     
+    --cxx)
+      echo @GRID_CXX@
+    ;;
+    
     --ldflags)
       echo @GRID_LDFLAGS@
     ;;

tests/debug/Test_split_laplacian.cc

@@ -0,0 +1,104 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_mrhs_cg.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    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 <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+int main (int argc, char ** argv)
+{
+  typedef LatticeComplex ComplexField; 
+
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  int nd   = latt_size.size();
+  int ndm1 = nd-1;
+
+  std::vector<int> simd_layout = GridDefaultSimd(nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  std::vector<int> mpi_split (mpi_layout.size(),1);
+
+  std::cout << " Full " << GridCmdVectorIntToString(latt_size)  << " subgrid"         <<std::endl;
+  std::cout << " Full " << GridCmdVectorIntToString(mpi_layout) << " sub communicator"<<std::endl;
+  std::cout << " Full " << GridCmdVectorIntToString(simd_layout)<< " simd layout "    <<std::endl;
+
+  GridCartesian         * GridN = new GridCartesian(latt_size,
+						    simd_layout,
+						    mpi_layout);
+
+  std::vector<int> latt_m  = latt_size;   latt_m[nd-1] = 1;
+  std::vector<int> mpi_m   = mpi_layout;  mpi_m [nd-1] = 1;
+  std::vector<int> simd_m  = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1);
+
+
+  std::cout << " Requesting " << GridCmdVectorIntToString(latt_m)<< " subgrid"         <<std::endl;
+  std::cout << " Requesting " << GridCmdVectorIntToString(mpi_m) << " sub communicator"<<std::endl;
+  std::cout << " Requesting " << GridCmdVectorIntToString(simd_m)<< " simd layout "    <<std::endl;
+  GridCartesian         * Grid_m = new GridCartesian(latt_m,
+						     simd_m,
+						     mpi_m,
+						     *GridN); 
+
+  Complex C(1.0);
+  Complex tmp;
+
+  ComplexField Full(GridN); Full = C;
+  ComplexField Full_cpy(GridN);
+  ComplexField Split(Grid_m);Split= C;
+
+  std::cout << GridLogMessage<< " Full  volume "<< norm2(Full) <<std::endl;
+  std::cout << GridLogMessage<< " Split volume "<< norm2(Split) <<std::endl;
+
+  tmp=C;
+  GridN->GlobalSum(tmp);
+  std::cout << GridLogMessage<< " Full  nodes "<< tmp <<std::endl;
+
+  tmp=C;
+  Grid_m->GlobalSum(tmp);
+  std::cout << GridLogMessage<< " Split nodes "<< tmp <<std::endl;
+  GridN->Barrier();
+
+  auto local_latt = GridN->LocalDimensions();
+
+  Full_cpy = zero;
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          RNG(GridN);  RNG.SeedFixedIntegers(seeds);
+
+  random(RNG,Full);
+  for(int t=0;t<local_latt[nd-1];t++){
+    ExtractSliceLocal(Split,Full,0,t,Tp);
+    InsertSliceLocal (Split,Full_cpy,0,t,Tp);
+  }
+  Full_cpy = Full_cpy - Full;
+  std::cout << " NormFull " << norm2(Full)<<std::endl;
+  std::cout << " NormDiff " << norm2(Full_cpy)<<std::endl;
+  Grid_finalize();
+}

tests/hadrons/Test_QED.cc

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
+    std::string twist    = "0. 0. 0. 0.";
 
     //stochastic photon field
     MGauge::StochEm::Par photonPar;
@@ -90,6 +91,7 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
+        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
 
         

tests/hadrons/Test_hadrons.hpp

@@ -126,6 +126,7 @@ inline void makeWilsonAction(Application &application, std::string actionName,
         actionPar.gauge = gaugeField;
         actionPar.mass  = mass;
         actionPar.boundary = boundary;
+        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::Wilson>(actionName, actionPar);
     }
 }
@@ -154,6 +155,7 @@ inline void makeDWFAction(Application &application, std::string actionName,
         actionPar.M5    = M5;
         actionPar.mass  = mass;
         actionPar.boundary = boundary;
+        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::DWF>(actionName, actionPar);
     }
 }

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -66,6 +66,7 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
+    std::string twist = "0. 0. 0. 0.";
 
     // sink
     MSink::Point::Par sinkPar;
@@ -80,6 +81,7 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
+        actionPar.twist = twist;
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
         
         // solvers

tests/hadrons/Test_hadrons_spectrum.cc

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
+    std::string twist = "0. 0. 0. 0.";
 
     for (unsigned int i = 0; i < flavour.size(); ++i)
     {
@@ -82,6 +83,7 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
+        actionPar.twist = twist;
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
         
         // solvers