rpg/System.h

#ifndef RPG_SYSTEM_H
#define RPG_SYSTEM_H
 
/*
* PSCF - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
// Header file includes
#include <util/param/ParamComposite.h>     // base class
#include <rpg/solvers/Mixture.h>           // member
#include <rpg/field/Domain.h>              // member
#include <rpg/field/FieldIo.h>             // member
#include <rpg/field/WFieldContainer.h>     // member
#include <rpg/field/CFieldContainer.h>     // member
#include <rpg/field/Mask.h>                // member
#include <prdc/cuda/RField.h>              // member (tmpFieldsRGrid_)
#include <prdc/cuda/RFieldDft.h>           // member (tmpFieldsKGrid_)
#include <pscf/homogeneous/Mixture.h>      // member
#include <util/misc/FileMaster.h>          // member
#include <util/containers/DArray.h>        // member (tmpFields ...)
//#include <util/containers/FSArray.h>       // ????
 
// Forward references
namespace Util {
   template <typename T, int N> class FSArray;
}
namespace Pscf {
   class Interaction;
   template <typename Data> class DeviceArray;
   namespace Prdc {
      template <int D> class UnitCell;
   }
   namespace Rpg {
      template <int D> class Iterator;
      template <int D> class IteratorFactory;
      template <int D> class Sweep;
      template <int D> class SweepFactory;
      template <int D> class Simulator;
      template <int D> class SimulatorFactory;
   }
}
 
 
namespace Pscf {
namespace Rpg {
 
   // Namespaces that are implicitly available, without qualification
   using namespace Util;
   using namespace Pscf::Prdc;
   using namespace Pscf::Prdc::Cuda;
 
   template <int D>
   class System : public ParamComposite
   {
 
   public:
 
 
      System();
 
      ~System();
 
 
      void setOptions(int argc, char **argv);
 
      virtual void readParam(std::istream& in);
 
      void readParam();
 
      virtual void readParameters(std::istream& in);
 
      void readCommands(std::istream& in);
 
      void readCommands();
 
 
      void readWBasis(const std::string & filename);
 
      void readWRGrid(const std::string & filename);
 
      void setWBasis(DArray< DArray<double> > const & fields);
 
      void setWRGrid(DArray< RField<D> > const & fields);
 
      void setWRGrid(DeviceArray<cudaReal> & fields);
 
      void symmetrizeWFields();
 
      void estimateWfromC(const std::string& filename);
 
 
      void setUnitCell(UnitCell<D> const & unitCell);
 
      void setUnitCell(typename UnitCell<D>::LatticeSystem lattice,
                       FSArray<double, 6> const & parameters);
 
      void setUnitCell(FSArray<double, 6> const & parameters);
 
 
      void compute(bool needStress = false);
 
      int iterate(bool isContinuation = false);
 
      void sweep();
 
      void simulate(int nStep);
 
 
      void computeFreeEnergy();
 
      double fHelmholtz() const;
 
      double pressure() const;
 
 
      void writeParamNoSweep(std::ostream& out) const;
 
      void writeThermo(std::ostream& out);
 
      void writeStress(std::ostream& out);
 
 
      void writeWBasis(const std::string & filename);
 
      void writeWRGrid(const std::string & filename) const;
 
      void writeCBasis(const std::string & filename);
 
      void writeCRGrid(const std::string & filename) const;
 
      void writeBlockCRGrid(const std::string & filename) const;
 
 
      void writeQSlice(std::string const & filename,
                       int polymerId, int blockId,
                       int directionId, int segmentId)  const;
 
      void writeQTail(std::string const & filename, int polymerId,
                      int blockId, int directionId)  const;
 
      void writeQ(std::string const & filename, int polymerId,
                  int blockId, int directionId)  const;
 
      void writeQAll(std::string const & basename);
 
 
      void writeStars(const std::string & filename) const;
 
      void writeWaves(const std::string & filename) const;
 
      void writeGroup(std::string const & filename) const;
 
 
      void basisToRGrid(const std::string & inFileName,
                        const std::string & outFileName);
 
      void rGridToBasis(const std::string & inFileName,
                        const std::string & outFileName);
 
      void kGridToRGrid(const std::string& inFileName,
                        const std::string& outFileName);
 
      void rGridToKGrid(const std::string & inFileName,
                        const std::string & outFileName);
 
      void kGridToBasis(const std::string& inFileName,
                        const std::string& outFileName);
 
      void basisToKGrid(const std::string & inFileName,
                        const std::string & outFileName);
 
      void compare(const DArray< DArray<double> > field1,
                   const DArray< DArray<double> > field2);
 
      void compare(const DArray< RField<D> > field1,
                   const DArray< RField<D> > field2);
 
      bool checkRGridFieldSymmetry(const std::string & inFileName,
                                   double epsilon = 1.0E-8);
 
      void scaleFieldsBasis(const std::string & inFileName,
                            const std::string & outFileName,
                            double factor);
 
      void scaleFieldsRGrid(const std::string & inFileName,
                            const std::string & outFileName,
                            double factor) const;
 
      void expandRGridDimension(const std::string & inFileName,
                                const std::string & outFileName,
                                int d,
                                DArray<int> newGridDimensions);
      void replicateUnitCell(const std::string & inFileName,
                             const std::string & outFileName,
                             IntVec<D> const & replicas);
 
 
      void writeTimers(std::ostream& out);
 
      void clearTimers();
 
 
      WFieldContainer<D> const & w() const;
 
      CFieldContainer<D> const & c() const;
 
      WFieldContainer<D>& h();
 
      WFieldContainer<D> const & h() const;
 
      Mask<D>& mask();
 
      Mask<D> const & mask() const;
 
 
      Mixture<D>& mixture();
 
      Mixture<D> const & mixture() const;
 
      Interaction & interaction();
   
      Interaction const & interaction() const;
   
      Domain<D> & domain();
 
      Domain<D> const & domain() const;
 
      Iterator<D>& iterator();
 
      Iterator<D> const & iterator() const;
 
      Simulator<D>& simulator();
 
      UnitCell<D> const & unitCell() const;
 
      Mesh<D> const & mesh() const;
 
      Basis<D> const & basis() const;
 
      FieldIo<D> const & fieldIo() const;
 
      FFT<D> const & fft() const ;
 
      Homogeneous::Mixture& homogeneous();
 
      FileMaster& fileMaster();
 
 
      bool hasIterator() const;
 
      bool hasSweep() const;
 
      bool hasSimulator() const;
 
      bool hasExternalFields() const;
 
      bool hasMask() const;
 
      bool hasCFields() const;
 
      bool hasFreeEnergy() const;
 
 
   private:
 
      Mixture<D> mixture_;
 
      Domain<D> domain_;
 
      FileMaster fileMaster_;
 
      Homogeneous::Mixture homogeneous_;
 
      Interaction* interactionPtr_;
   
      Iterator<D>* iteratorPtr_;
 
      IteratorFactory<D>* iteratorFactoryPtr_;
 
      Sweep<D>* sweepPtr_;
 
      SweepFactory<D>* sweepFactoryPtr_;
 
      Simulator<D>* simulatorPtr_;
 
      SimulatorFactory<D>* simulatorFactoryPtr_;
 
      WFieldContainer<D> w_;
 
      CFieldContainer<D> c_;
 
      WFieldContainer<D> h_;
 
      Mask<D> mask_;
 
      mutable DArray< DArray<double> > tmpFieldsBasis_;
 
      mutable DArray< RField<D> > tmpFieldsRGrid_;
 
      mutable DArray<RFieldDft<D> > tmpFieldsKGrid_;
 
      double fHelmholtz_;
 
      double fIdeal_;
 
      double fInter_;
 
      double fExt_;
 
      double pressure_;
 
      bool isAllocatedGrid_;
 
      bool isAllocatedBasis_;
 
      bool hasMixture_;
 
      bool hasCFields_;
 
      bool hasFreeEnergy_;
 
      IntVec<D> kMeshDimensions_;
 
      RField<D> workArray_;
 
      // Private member functions
 
      void allocateFieldsGrid();
 
      void allocateFieldsBasis();
 
      void readFieldHeader(std::string filename);
 
      void readEcho(std::istream& in, std::string& string) const;
 
      void readEcho(std::istream& in, double& value) const;
 
      void initHomogeneous();
 
   };
 
   // Inline member functions
 
   // Get the associated Mixture<D> object.
   template <int D>
   inline Mixture<D>& System<D>::mixture()
   { return mixture_; }
 
   // Get the associated Mixture<D> object by const reference.
   template <int D>
   inline Mixture<D> const & System<D>::mixture() const
   { return mixture_; }
 
   // Get the Interaction (excess free energy) by non-const reference.
   template <int D>
   inline Interaction & System<D>::interaction()
   {
      UTIL_ASSERT(interactionPtr_);
      return *interactionPtr_;
   }
 
   // Get the Interaction (excess free energy ) by const reference.
   template <int D>
   inline Interaction const & System<D>::interaction() const
   {
      UTIL_ASSERT(interactionPtr_);
      return *interactionPtr_;
   }
 
   // Get the Domain by non-const reference.
   template <int D>
   inline Domain<D> & System<D>::domain()
   {  return domain_; }
 
   // Get the Domain by const reference.
   template <int D>
   inline Domain<D> const & System<D>::domain() const
   {  return domain_; }
 
   template <int D>
   inline UnitCell<D> const & System<D>::unitCell() const
   {  return domain_.unitCell(); }
 
   // Get the Mesh<D> by const reference.
   template <int D>
   inline Mesh<D> const & System<D>::mesh() const
   {  return domain_.mesh(); }
 
   // Get the Basis<D>  by const reference.
   template <int D>
   inline Basis<D> const & System<D>::basis() const
   {  return domain_.basis(); }
 
   // Get the FFT<D> object by const reference.
   template <int D>
   inline FFT<D> const & System<D>::fft() const
   {  return domain_.fft(); }
 
   // Get the const FieldIo<D> object by const reference.
   template <int D>
   inline FieldIo<D> const & System<D>::fieldIo() const
   {  return domain_.fieldIo(); }
 
   // Get the Iterator by non-const reference.
   template <int D>
   inline Iterator<D>& System<D>::iterator()
   {
      UTIL_ASSERT(iteratorPtr_);
      return *iteratorPtr_;
   }
 
   // Get the Iterator by const reference.
   template <int D>
   inline Iterator<D> const & System<D>::iterator() const
   {
      UTIL_ASSERT(iteratorPtr_);
      return *iteratorPtr_;
   }
 
   // Get the Simulator by non-const reference.
   template <int D>
   inline Simulator<D>& System<D>::simulator()
   {
      UTIL_ASSERT(simulatorPtr_);
      return *simulatorPtr_;
   }
 
   // Get the FileMaster.
   template <int D>
   inline FileMaster& System<D>::fileMaster()
   {  return fileMaster_; }
 
   // Get the Homogeneous::Mixture object.
   template <int D>
   inline
   Homogeneous::Mixture& System<D>::homogeneous()
   {  return homogeneous_; }
 
   // Get container of w fields (const reference)
   template <int D>
   inline
   WFieldContainer<D> const & System<D>::w() const
   {  return w_; }
 
   // Get container of c fields (const reference)
   template <int D>
   inline
   CFieldContainer<D> const & System<D>::c() const
   {  return c_; }
 
   // Get container of external potential fields (reference)
   template <int D>
   inline WFieldContainer<D>& System<D>::h()
   {  return h_; }
 
   // Get container of external potential fields (const reference)
   template <int D>
   inline WFieldContainer<D> const & System<D>::h() const
   {  return h_; }
 
   // Get mask field (reference)
   template <int D>
   inline Mask<D>& System<D>::mask()
   {  return mask_; }
 
   // Get mask field (const reference)
   template <int D>
   inline Mask<D> const & System<D>::mask() const
   {  return mask_; }
 
   // Get precomputed Helmoltz free energy per monomer / kT.
   template <int D>
   inline double System<D>::fHelmholtz() const
   {
      UTIL_CHECK(hasFreeEnergy_);
      return fHelmholtz_;
   }
 
   // Get precomputed pressure (units of kT / monomer volume).
   template <int D>
   inline double System<D>::pressure() const
   {
      UTIL_CHECK(hasFreeEnergy_);
      return pressure_;
   }
 
   // Has the free energy been computed for the current w fields?
   template <int D>
   inline bool System<D>::hasFreeEnergy() const
   {  return hasFreeEnergy_; }
 
   // Have the c fields been computed for the current w fields?
   template <int D>
   inline bool System<D>::hasCFields() const
   {  return hasCFields_; }
 
   // Does the system have an Iterator object?
   template <int D>
   inline bool System<D>::hasIterator() const
   {  return (iteratorPtr_); }
 
   // Does this system have an associated Sweep object?
   template <int D>
   inline bool System<D>::hasSweep() const
   {  return (sweepPtr_); }
 
   // Does the system have an associated Simulator ?
   template <int D>
   inline bool System<D>::hasSimulator() const
   {  return (simulatorPtr_); }
 
   // Does this system have external potential fields?
   template <int D>
   inline bool System<D>::hasExternalFields() const
   {  return h_.hasData(); }
 
   // Does this system have a mask?
   template <int D>
   inline bool System<D>::hasMask() const
   {  return mask_.hasData(); }
 
   #ifndef RPG_SYSTEM_TPP
   // Suppress implicit instantiation
   extern template class System<1>;
   extern template class System<2>;
   extern template class System<3>;
   #endif
 
} // namespace Rpg
} // namespace Pscf
#endif