Adding more fermion action modules, generalised DWF

lib/qcd/hmc/GenericHMCrunner.h

@@ -36,7 +36,7 @@ namespace Grid {
 namespace QCD {
 
 
-// very ugly here but possibly resolved if we have a base Reader class
+// very ugly here but possibly resolved if we had a base Reader class
 template < class ReaderClass >
 class HMCRunnerBase {
 public:

lib/qcd/hmc/HMCModules.h

@@ -115,6 +115,7 @@ class GridFourDimModule : public GridModule {
   }
 };
 
+
 ////////////////////////////////////////////////////////////////////
 class RNGModuleParameters: Serializable {
 

lib/qcd/modules/ActionModules.h

@@ -239,7 +239,7 @@ class TwoFlavourFModule: public PseudoFermionModuleBase<Impl, TwoFlavourPseudoFe
 
    // constructor
    template <class ReaderClass>
-   TwoFlavourFModule(Reader<ReaderClass>& R): PseudoFermionModuleBase<Impl, TwoFlavourPseudoFermionAction>(R) {
+   TwoFlavourFModule(Reader<ReaderClass>& R): Base(R) {
     this->getSolverOperator(R, solver_mod, "Solver");
     this->getFermionOperator(R, fop_mod, "Operator");
    } 

lib/qcd/modules/Factory.h

@@ -35,11 +35,11 @@ namespace Grid{
 /******************************************************************************
  *                        abstract factory class                              *
  ******************************************************************************/
-template <typename T, typename ProductCreator>
+template <typename T, typename CreatorInput>
 class Factory
 {
 public:
-    typedef std::function< std::unique_ptr<T> (const ProductCreator&) > Func;
+    typedef std::function< std::unique_ptr<T> (const CreatorInput&) > Func;
 
     // constructor
     Factory(void) = default;
@@ -51,7 +51,7 @@ public:
     std::vector<std::string> getBuilderList(void) const;
     // factory
     std::unique_ptr<T> create(const std::string type,
-              								const ProductCreator& input) const;
+              								const CreatorInput& input) const;
 private:
     std::map<std::string, Func> builder_;
     virtual std::string obj_type() const = 0;
@@ -61,15 +61,15 @@ private:
  *                         template implementation                            *
  ******************************************************************************/
 // registration ////////////////////////////////////////////////////////////////
-template <typename T, typename ProductCreator>
+template <typename T, typename CreatorInput>
-void Factory<T, ProductCreator>::registerBuilder(const std::string type, const Func &f)
+void Factory<T, CreatorInput>::registerBuilder(const std::string type, const Func &f)
 {
     builder_[type] = f;
 }
 
 // get module list /////////////////////////////////////////////////////////////
-template <typename T, typename ProductCreator>
+template <typename T, typename CreatorInput>
-std::vector<std::string> Factory<T, ProductCreator>::getBuilderList(void) const
+std::vector<std::string> Factory<T, CreatorInput>::getBuilderList(void) const
 {
     std::vector<std::string> list;
     
@@ -82,9 +82,9 @@ std::vector<std::string> Factory<T, ProductCreator>::getBuilderList(void) const
 }
 
 // factory /////////////////////////////////////////////////////////////////////
-template <typename T, typename ProductCreator>
+template <typename T, typename CreatorInput>
-std::unique_ptr<T> Factory<T, ProductCreator>::create(const std::string type,
+std::unique_ptr<T> Factory<T, CreatorInput>::create(const std::string type,
-                                      								const ProductCreator& input) const
+                                      							const CreatorInput& input) const
 {
     Func func;
     

lib/qcd/modules/FermionOperatorModules.h

@@ -47,7 +47,7 @@ class FermionOperatorModule
 
 protected:
   std::unique_ptr< FOType<FermionImpl> > FOPtr;
-  std::vector< std::reference_wrapper<QCD::GridModule> >    GridRefs;
+  std::vector< QCD::GridModule* >    GridRefs;
  public:
   typedef HMCModuleBase< QCD::FermionOperator<FermionImpl> > Base;
   typedef typename Base::Product Product;
@@ -58,11 +58,21 @@ protected:
   FermionOperatorModule(Reader<ReaderClass>& Reader) : Parametrized<FOPar>(Reader){};
 
   void AddGridPair(QCD::GridModule &Mod){
-    if (GridRefs.size()>2){
+    if (GridRefs.size()>1){
       std::cout << GridLogError << "Adding too many Grids to the FermionOperatorModule" << std::endl;
       exit(1);
     }
-    GridRefs.push_back(Mod);
+    GridRefs.push_back(&Mod);
+
+    if (Ls()){
+      GridRefs.push_back(new QCD::GridModule());
+      GridRefs[1]->set_full(QCD::SpaceTimeGrid::makeFiveDimGrid(Ls(),GridRefs[0]->get_full()));
+      GridRefs[1]->set_rb(QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls(),GridRefs[0]->get_full()));
+    }
+  }
+
+  virtual unsigned int Ls(){
+    return 0;  
   }
 
   virtual void print_parameters(){
@@ -124,16 +134,86 @@ class WilsonFermionModule: public FermionOperatorModule<WilsonFermion, FermionIm
 
   // acquire resource
   virtual void initialize(){
-    auto &GridMod = this->GridRefs[0].get();
+    auto GridMod = this->GridRefs[0];
-    typename FermionImpl::GaugeField U(GridMod.get_full());
+    typename FermionImpl::GaugeField U(GridMod->get_full());
-    this->FOPtr.reset(new WilsonFermion<FermionImpl>(U, *(GridMod.get_full()), *(GridMod.get_rb()), this->Par_.mass));
+    this->FOPtr.reset(new WilsonFermion<FermionImpl>(U, *(GridMod->get_full()), *(GridMod->get_rb()), this->Par_.mass));
   }
 };
 
 
+
+class MobiusFermionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MobiusFermionParameters, 
+    RealD, mass,
+    RealD, M5,
+    RealD, b,
+    RealD, c,
+    unsigned int, Ls);
+};
+
+template <class FermionImpl >
+class MobiusFermionModule: public FermionOperatorModule<MobiusFermion, FermionImpl, MobiusFermionParameters> {
+  typedef FermionOperatorModule<MobiusFermion, FermionImpl, MobiusFermionParameters> FermBase;
+  using FermBase::FermBase; // for constructors
+
+  virtual unsigned int Ls(){
+    return this->Par_.Ls;
+  }
+
+  // acquire resource
+  virtual void initialize(){
+    auto GridMod = this->GridRefs[0];
+    auto GridMod5d = this->GridRefs[1];
+    typename FermionImpl::GaugeField U(GridMod->get_full());
+    this->FOPtr.reset(new MobiusFermion<FermionImpl>( U, *(GridMod->get_full()), *(GridMod->get_rb()), 
+                                                      *(GridMod5d->get_full()), *(GridMod5d->get_rb()),
+                                                      this->Par_.mass, this->Par_.M5, this->Par_.b, this->Par_.c));
+  }
+};
+
+
+class DomainWallFermionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(DomainWallFermionParameters, 
+    RealD, mass,
+    RealD, M5,
+    unsigned int, Ls);
+};
+
+template <class FermionImpl >
+class DomainWallFermionModule: public FermionOperatorModule<DomainWallFermion, FermionImpl, DomainWallFermionParameters> {
+  typedef FermionOperatorModule<DomainWallFermion, FermionImpl, DomainWallFermionParameters> FermBase;
+  using FermBase::FermBase; // for constructors
+
+  virtual unsigned int Ls(){
+    return this->Par_.Ls;
+  }
+
+  // acquire resource
+  virtual void initialize(){
+    auto GridMod = this->GridRefs[0];
+    auto GridMod5d = this->GridRefs[1];
+    typename FermionImpl::GaugeField U(GridMod->get_full());
+    this->FOPtr.reset(new DomainWallFermion<FermionImpl>( U, *(GridMod->get_full()), *(GridMod->get_rb()), 
+                                                      *(GridMod5d->get_full()), *(GridMod5d->get_rb()),
+                                                      this->Par_.mass, this->Par_.M5));
+  }
+};
+
+
+
 // Now a specific registration with a fermion field
 static Registrar< WilsonFermionModule<WilsonImplR>,   
                   HMC_FermionOperatorModuleFactory<fermionop_string, WilsonImplR, XmlReader> > __WilsonFOPmodXMLInit("Wilson"); 
+static Registrar< MobiusFermionModule<WilsonImplR>,   
+                  HMC_FermionOperatorModuleFactory<fermionop_string, WilsonImplR, XmlReader> > __MobiusFOPmodXMLInit("Mobius");
+static Registrar< DomainWallFermionModule<WilsonImplR>,   
+                  HMC_FermionOperatorModuleFactory<fermionop_string, WilsonImplR, XmlReader> > __DWFOPmodXMLInit("DomainWall");
+
+
+
+
 
 
 } // QCD

tests/hmc/Test_hmc_Factories.cc

@@ -37,7 +37,8 @@ typedef QCD::PeriodicGimplR ImplementationPolicy;
 typedef QCD::NoHirep RepresentationPolicy;
 typedef QCD::WilsonFermionR FermionImplementation;
 
-
+///////////////////////////////////////////////////////////////////////
+// Put all registrations in an header file
 static Registrar< HMCLeapFrog<ImplementationPolicy, RepresentationPolicy, XmlReader>      , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCLFmodXMLInit("LeapFrog");
 static Registrar< HMCMinimumNorm2<ImplementationPolicy, RepresentationPolicy, XmlReader>  , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCMN2modXMLInit("MinimumNorm2");
 static Registrar< HMCForceGradient<ImplementationPolicy, RepresentationPolicy, XmlReader> , HMCRunnerModuleFactory<hmc_string, XmlReader> > __HMCFGmodXMLInit("ForceGradient");
@@ -58,51 +59,17 @@ int main(int argc, char **argv) {
   // Reader, file should come from command line
   InputFileReader Reader("input.wilson_gauge.params.xml");
 
-
-
   // Test HMC factory (put in an external file)
   auto &HMCfactory = HMCRunnerModuleFactory<hmc_string, InputFileReader >::getInstance();
   // Simplify this step (IntergratorName field?)
   HMCparameters HMCpar(Reader);
   
   // Construct the module
-  auto myHMCmodule = HMCfactory.create(HMCpar.MD.name, Reader);
+  auto HMCmodule = HMCfactory.create(HMCpar.MD.name, Reader);
 
-  myHMCmodule->getPtr()->initialize(Reader);
+  HMCmodule->getPtr()->initialize(Reader);
-  myHMCmodule->getPtr()->Run();
+  HMCmodule->getPtr()->Run();
 
   Grid_finalize();
 
-
-/*
-// Test solver creation
-  auto &SolverFactory = HMC_SolverModuleFactory<solver_string, FermionImplementation::FermionField, XmlReader>::getInstance();
-  Reader.push("Solver");
-  std::string name;
-  read(Reader, "name",name);
-  auto SModule = SolverFactory.create(name, Reader);
-  std::cout << "Registered types " << std::endl;
-  std::cout << SolverFactory.getBuilderList() << std::endl;
-  SModule->print_parameters();
-  Reader.pop();
-*/
-
-/*
- // Test fermion operator module creation
-  auto &FOFactory = HMC_FermionOperatorModuleFactory<fermionop_string, WilsonImplR, XmlReader>::getInstance();
-  Reader.push("Operator");
-  std::string op_name;
-  Reader.readDefault("name",op_name);
-  auto FOModule = FOFactory.create(op_name, Reader);
-  std::cout << "Registered types " << std::endl;
-  std::cout << FOFactory.getBuilderList() << std::endl;
-  GridFourDimModule GMod;
-  FOModule->AddGridPair(GMod);
-  FOModule->print_parameters();
-  Reader.pop();  
-*/
-  
-
-
-
 } // main

tests/hmc/Test_hmc_WilsonFermionGauge_Binary.cc

@@ -84,11 +84,11 @@ int main(int argc, char **argv) {
   // standard
   RealD beta = 5.6 ;
   WilsonGaugeActionR Waction(beta);
-  
+    
-  // temporarily need a gauge field
   auto GridPtr = TheHMC.Resources.GetCartesian();
   auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
 
+  // temporarily need a gauge field
   LatticeGaugeField U(GridPtr);
 
   Real mass = -0.77;
@@ -101,6 +101,13 @@ int main(int argc, char **argv) {
 
   TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
 
+  // With modules
+  /*
+
+  TwoFlavourFmodule<FermionImplPolicy> TwoFMod(Reader);
+  
+  */
+
     // Set smearing (true/false), default: false
   Nf2.is_smeared = false;