rpc/fts/simulator/Simulator.h

#ifndef RPC_SIMULATOR_H
#define RPC_SIMULATOR_H
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <util/param/ParamComposite.h>     // base class
 
#include <rpc/fts/simulator/SimState.h>    // member
#include <prdc/cpu/RField.h>               // member (template arg)
#include <util/random/Random.h>            // member
#include <util/containers/DArray.h>        // member (template)
#include <util/containers/DMatrix.h>       // member (template)
 
namespace Pscf {
namespace Rpc {
 
   // Forward declarations 
   template <int D> class System;
   template <int D> class Compressor;
   template <int D> class CompressorFactory;
   template <int D> class Perturbation;
   template <int D> class PerturbationFactory;
   template <int D> class Ramp;
   template <int D> class RampFactory;
 
   using namespace Util;
   using namespace Prdc;
   using namespace Prdc::Cpu;

   template <int D>
   class Simulator : public ParamComposite
   {
 
   public:
 
      Simulator() = delete;

      Simulator(System<D>& system);

      ~Simulator();

      void allocate();

      virtual void readParameters(std::istream &in);


      virtual void simulate(int nStep);

      virtual void analyze(int min, int max,
                           std::string classname,
                           std::string filename);

      void clearData();


      virtual void outputTimers(std::ostream& out) const;

      virtual void outputMdeCounter(std::ostream& out) const;

      virtual void clearTimers();

      long iStep();
      
      long iTotalStep();


      void analyzeChi();

      DArray<double> const & chiEvals() const;

      double chiEval(int a) const;

      DMatrix<double> const & chiEvecs() const;

      double chiEvecs(int a, int i) const;

      DArray<double> const & sc() const;

      double sc(int a) const;


      void computeHamiltonian();

      double hamiltonian() const;

      double idealHamiltonian() const;

      double fieldHamiltonian() const;

      double perturbationHamiltonian() const;

      bool hasHamiltonian() const;


      void computeWc();

      DArray< RField<D> > const & wc() const;

      RField<D> const & wc(int a) const;

      bool hasWc() const;


      void computeCc();

      DArray< RField<D> > const & cc() const;

      RField<D> const & cc(int a) const;

      bool hasCc() const;


      void computeDc();

      DArray< RField<D> > const & dc() const;

      RField<D> const & dc(int i) const;

      bool hasDc() const;
      
      
      void saveState();
      
      void restoreState();
      
      void clearState();
      

      System<D>& system();

      Random& random();
      
      bool hasCompressor() const;

      Compressor<D>& compressor();

      Compressor<D> const & compressor() const;

      bool hasPerturbation() const;

      Perturbation<D>& perturbation();

      Perturbation<D> const & perturbation() const;

      bool hasRamp() const;

      Ramp<D> const & ramp() const;
      
      Ramp<D>& ramp();

 
   protected:
 
      // Protected member functions
 
      using Util::ParamComposite::setClassName;

      void readRandomSeed(std::istream& in);

      CompressorFactory<D>& compressorFactory();

      void readCompressor(std::istream& in, bool& isEnd);

      PerturbationFactory<D>& perturbationFactory();

      void readPerturbation(std::istream& in, bool& isEnd);

      void setPerturbation(Perturbation<D>* ptr);

      RampFactory<D>& rampFactory();

      void readRamp(std::istream& in, bool& isEnd);

      void setRamp(Ramp<D>* ptr);
 
      // Protected data members

      Random random_;

      DArray< RField<D> > wc_;

      DArray< RField<D> > cc_;

      DArray< RField<D> > dc_;
      
      mutable SimState<D> state_;

      double hamiltonian_;

      double idealHamiltonian_;

      double fieldHamiltonian_;

      double perturbationHamiltonian_;

      long iStep_;
      
      long iTotalStep_;

      long seed_;
 
      // Boolean status indicators

      bool hasHamiltonian_;

      bool hasWc_;

      bool hasCc_;

      bool hasDc_;
 
   private:

      DMatrix<double> chiP_;

      DMatrix<double> chiEvecs_;

      DArray<double>  chiEvals_;

      DArray<double>  sc_;
      
      System<D>* systemPtr_;

      CompressorFactory<D>* compressorFactoryPtr_;

      Compressor<D>* compressorPtr_;

      PerturbationFactory<D>* perturbationFactoryPtr_;

      Perturbation<D>* perturbationPtr_;

      RampFactory<D>* rampFactoryPtr_;

      Ramp<D>* rampPtr_;

      bool isAllocated_;
 
   };
 
   // Inline functions
 
   // Management of owned and associated objects
 
   // Get the parent System by reference.
   template <int D>
   inline System<D>& Simulator<D>::system()
   {
      UTIL_ASSERT(systemPtr_);  
      return *systemPtr_; 
   }
 
   // Get the random number generator by reference.
   template <int D>
   inline Random& Simulator<D>::random()
   {  return random_; }
 
   // Get the compressor factory by reference.
   template <int D>
   inline CompressorFactory<D>& Simulator<D>::compressorFactory()
   {
      UTIL_CHECK(compressorFactoryPtr_);  
      return *compressorFactoryPtr_; 
   }
 
   // Does this Simulator have an associated Compressor?
   template <int D>
   inline bool Simulator<D>::hasCompressor() const
   {  return (bool)compressorPtr_; }
   
   // Get the Compressor by non-const reference.
   template <int D>
   inline Compressor<D>& Simulator<D>::compressor()
   {
      UTIL_CHECK(hasCompressor());
      return *compressorPtr_;
   }
 
   // Get the Compressor by const reference.
   template <int D>
   inline Compressor<D> const& Simulator<D>::compressor() const
   {
      UTIL_CHECK(hasCompressor());
      return *compressorPtr_;
   }
 
   // Does this Simulator have an associated Perturbation?
   template <int D>
   inline bool Simulator<D>::hasPerturbation() const
   {  return (bool)perturbationPtr_; }
   
   // Get the perturbation (if any) by const reference.
   template <int D>
   inline Perturbation<D> const & Simulator<D>::perturbation() const
   {
      UTIL_CHECK(perturbationPtr_);  
      return *perturbationPtr_; 
   }
 
   // Get the perturbation (if any) by non-const reference.
   template <int D>
   inline Perturbation<D>& Simulator<D>::perturbation()
   {
      UTIL_CHECK(perturbationPtr_);  
      return *perturbationPtr_; 
   }
 
   // Get the perturbation factory by reference.
   template <int D>
   inline PerturbationFactory<D>& Simulator<D>::perturbationFactory()
   {
      UTIL_CHECK(perturbationFactoryPtr_);  
      return *perturbationFactoryPtr_; 
   }
 
   // Does this Simulator have an associated Ramp?
   template <int D>
   inline bool Simulator<D>::hasRamp() const
   {  return (bool)rampPtr_; }
   
   // Get the ramp (if any) by const reference.
   template <int D>
   inline Ramp<D> const & Simulator<D>::ramp() const
   {
      UTIL_CHECK(rampPtr_);  
      return *rampPtr_; 
   }
 
   // Get the ramp (if any) by non-const reference.
   template <int D>
   inline Ramp<D>& Simulator<D>::ramp()
   {
      UTIL_CHECK(rampPtr_);  
      return *rampPtr_; 
   }
 
   // Get the ramp factory.
   template <int D>
   inline RampFactory<D>& Simulator<D>::rampFactory()
   {
      UTIL_CHECK(rampFactoryPtr_);  
      return *rampFactoryPtr_; 
   }
 
   // Projected Chi Matrix
 
   // Return an array of eigenvalues of the projected chi matrix.
   template <int D>
   inline DArray<double> const & Simulator<D>::chiEvals() const
   {  return chiEvals_; }
 
   // Return a single eigenvalue of the projected chi matrix.
   template <int D>
   inline double Simulator<D>::chiEval(int a) const
   {  return chiEvals_[a]; }
 
   // Return a matrix of eigenvectors of the projected chi matrix.
   template <int D>
   inline DMatrix<double> const & Simulator<D>::chiEvecs() const
   {  return chiEvecs_; }
 
   // Return an element of an eigenvector of the projected chi matrix.
   template <int D>
   inline double Simulator<D>::chiEvecs(int a, int i) const
   {  return chiEvecs_(a, i); }
 
   // Return the vector S in eigenvector basis.
   template <int D>
   inline DArray<double> const & Simulator<D>::sc() const
   {  return sc_; }
 
   // Return one component of vector S in an eigenvector basis.
   template <int D>
   inline double Simulator<D>::sc(int a) const
   {  return sc_[a]; }
 
   // Hamiltonian and its derivatives
 
   // Get the precomputed Hamiltonian.
   template <int D>
   inline double Simulator<D>::hamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return hamiltonian_;
   }
 
   // Get the ideal gas component of the precomputed Hamiltonian.
   template <int D>
   inline double Simulator<D>::idealHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return idealHamiltonian_;
   }
 
   // Get the W field component of the precomputed Hamiltonian.
   template <int D>
   inline double Simulator<D>::fieldHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return fieldHamiltonian_;
   }
 
   // Get the perturbation component of the precomputed Hamiltonian.
   template <int D>
   inline double Simulator<D>::perturbationHamiltonian() const
   {
      UTIL_CHECK(hasHamiltonian_);
      return perturbationHamiltonian_;
   }
 
   // Has the Hamiltonian been computed for the current w fields ?
   template <int D>
   inline bool Simulator<D>::hasHamiltonian() const
   {  return hasHamiltonian_; }
 
   // Fields
 
   // Return all eigencomponents of the w fields.
   template <int D>
   inline DArray< RField<D> > const & Simulator<D>::wc() const
   {  return wc_; }
 
   // Return a single eigenvector component of the w fields.
   template <int D>
   inline RField<D> const & Simulator<D>::wc(int a) const
   {  return wc_[a]; }
 
   // Have eigenvector components of current w fields been computed?
   template <int D>
   inline bool Simulator<D>::hasWc() const
   {  return hasWc_; }
 
   // Return all eigenvector components of the current c fields.
   template <int D>
   inline DArray< RField<D> > const & Simulator<D>::cc() const
   {  return cc_; }
 
   // Return a single eigenvector component of the current c fields.
   template <int D>
   inline RField<D> const & Simulator<D>::cc(int a) const
   {  return cc_[a]; }
 
   // Have eigenvector components of current c fields been computed?
   template <int D>
   inline bool Simulator<D>::hasCc() const
   {  return hasCc_; }
 
   // Return all eigenvector components of the current d fields.
   template <int D>
   inline DArray< RField<D> > const & Simulator<D>::dc() const
   {  return dc_; }
 
   // Return a single eigenvector component of the current d fields.
   template <int D>
   inline RField<D> const & Simulator<D>::dc(int a) const
   {  return dc_[a]; }
 
   // Have eigenvector components of current d fields been computed?
   template <int D>
   inline bool Simulator<D>::hasDc() const
   {  return hasDc_; }
 
   // Return the current converged simulation step index.
   template <int D>
   inline long Simulator<D>::iStep()
   {  return iStep_; }
   
   // Return the current simulation step index.
   template <int D>
   inline long Simulator<D>::iTotalStep()
   {  return iTotalStep_; }
 
   #ifndef RPC_SIMULATOR_TPP
   // Suppress implicit instantiation
   extern template class Simulator<1>;
   extern template class Simulator<2>;
   extern template class Simulator<3>;
   #endif
 
}
}
#endif