r1d/system/System.h

#ifndef R1D_SYSTEM_H
#define R1D_SYSTEM_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.
*/
 
// Header file includes
#include <util/param/ParamComposite.h>    // base class
#include <r1d/field/FieldIo.h>             // member
#include <r1d/solvers/Mixture.h>          // member
#include <r1d/field/Domain.h>            // member
#include <util/misc/FileMaster.h>         // member
#include <util/containers/DArray.h>       // member template
#include <util/containers/Array.h>        // function parameter
 
// Forward declarations
namespace Pscf {
   class Interaction;
   namespace R1d {
      class Iterator;
      class IteratorFactory;
      class Sweep;
      class SweepFactory;
      using namespace Util;
   }
}
 
namespace Pscf {
namespace R1d {

   class System : public ParamComposite
   {
 
   public:

      typedef DArray<double> Field;

      typedef DArray<double> WField;

      typedef DArray<double> CField;


      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 compute();
   
      int iterate(bool isContinuation = false);
   
      void sweep();



      void computeFreeEnergy();

      double fHelmholtz() const;

      double pressure() const;



      void writeParamNoSweep(std::ostream& out) const;

      void writeThermo(std::ostream& out) const;


      
      void writeW(std::string const & filename) const;
   
      void writeC(std::string const & filename) const;
   
      void writeBlockC(std::string  const & 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) const;



      DArray<WField>& wFields();

      DArray<WField> const & wFields() const;

      WField& wField(int monomerId);

      WField const & wField(int monomerId) const;

      DArray<Field>& cFields();

      DArray<Field> const & cFields() const;

      Field& cField(int monomerId);

      Field const & cField(int monomerId) const;


      Mixture& mixture();

      Mixture const & mixture() const;

      Interaction & interaction();

      Interaction const & interaction() const;

      Domain& domain();

      Iterator& iterator();

      FileMaster& fileMaster();

 
   private:

      Mixture mixture_;

      Domain domain_;

      FileMaster fileMaster_;

      FieldIo fieldIo_;

      Interaction* interactionPtr_;

      Iterator* iteratorPtr_;

      IteratorFactory* iteratorFactoryPtr_;

      Sweep* sweepPtr_;

      SweepFactory* sweepFactoryPtr_;

      DArray<WField> wFields_;

      DArray<Field> cFields_;

      mutable DArray<double> f_;

      mutable DArray<double> c_;

      double fHelmholtz_;

      double fIdeal_;

      double fInter_;

      double pressure_;

      bool hasMixture_;

      bool hasDomain_;

      bool hasFields_;

      bool hasSweep_;
 
      // Private member functions

      void allocateFields();

      void readEcho(std::istream& in, std::string& string) const;
 
   };
 
   // Inline member functions
 
   /*
   * Get the associated Mixture object by reference.
   */
   inline Mixture& System::mixture()
   { return mixture_; }
 
   /*
   * Get the associated Mixture object by const reference.
   */
   inline Mixture const & System::mixture() const
   { return mixture_; }
 
   /*
   * Get the Interaction (excess free energy model).
   */
   inline Interaction & System::interaction()
   {
      UTIL_ASSERT(interactionPtr_);
      return *interactionPtr_;
   }
 
   /*
   * Get the Interaction (excess free energy) by const reference.
   */
   inline Interaction const & System::interaction() const
   {
      UTIL_ASSERT(interactionPtr_);
      return *interactionPtr_;
   }
 
   /*
   * Get the spatial Domain.
   */
   inline Domain& System::domain()
   { return domain_; }
 
   /*
   * Get the Iterator.
   */
   inline Iterator& System::iterator()
   {
      UTIL_ASSERT(iteratorPtr_);
      return *iteratorPtr_;
   }
 
   /*
   * Get the FileMaster.
   */
   inline FileMaster& System::fileMaster()
   {  return fileMaster_; }
 
   /*
   * Get an array of all monomer excess chemical potential fields.
   */
   inline 
   DArray< System::WField >& System::wFields()
   {  return wFields_; }
 
   /*
   * Get a const array of all monomer excess chemical potential fields.
   */
   inline 
   DArray< System::WField > const & System::wFields() const
   {  return wFields_; }
 
   /*
   * Get a single monomer excess chemical potential field.
   */
   inline 
   System::WField& System::wField(int id)
   {  return wFields_[id]; }
 
   /*
   * Get a single const monomer excess chemical potential field.
   */
   inline 
   System::WField const & System::wField(int id) const
   {  return wFields_[id]; }
 
   /*
   * Get array of all monomer concentration fields.
   */
   inline
   DArray< System::Field >& System::cFields()
   {  return cFields_; }
 
   /*
   * Get const array of all monomer concentration fields.
   */
   inline
   DArray< System::Field > const & System::cFields() const
   {  return cFields_; }
 
   /*
   * Get a single monomer concentration field.
   */
   inline System::Field& System::cField(int id)
   {  return cFields_[id]; }
 
   /*
   * Get a single const monomer concentration field.
   */
   inline System::Field const & System::cField(int id) const
   {  return cFields_[id]; }
 
   /*
   * Get precomputed Helmoltz free energy per monomer / kT.
   */
   inline double System::fHelmholtz() const
   {  return fHelmholtz_; }
 
   /*
   * Get precomputed pressure (units of kT / monomer volume).
   */
   inline double System::pressure() const
   {  return pressure_; }
 
} // namespace R1d
} // namespace Pscf
#endif