r1d/sweep/Sweep.cpp

/*
* 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 "Sweep.h"
//#include <pscf/sweep/SweepTmpl.h>
#include <r1d/System.h>
#include <r1d/domain/Domain.h>
#include <r1d/solvers/Mixture.h>
#include <r1d/iterator/Iterator.h>
#include <util/misc/ioUtil.h>
#include <util/format/Int.h>
#include <util/format/Dbl.h>
 
namespace Pscf {
namespace R1d
{
 
   using namespace Util;
 
   // Define maximum number of stored states = order of continuation + 1
   #define R1D_HISTORY_CAPACITY 3
 
   /*
   * Constructor.
   */
   Sweep::Sweep(System& sys)
    : Base(R1D_HISTORY_CAPACITY),
      SystemAccess(sys),
      homogeneousMode_(-1),
      comparison_(sys),
      fieldIo_()
   {
      setClassName("Sweep"); 
      fieldIo_.associate(sys.domain(), sys.fileMaster());
 
      addParameterTypes(system().iterator().getParameterTypes());
   }
 
   /*
   * Destructor.
   */
   Sweep::~Sweep()
   {}
 
   /*
   * Read parameters.
   */
   void Sweep::readParameters(std::istream& in)
   {
      Base::readParameters(in);
      homogeneousMode_ = -1; // default value
      readOptional<int>(in, "homogeneousMode", homogeneousMode_);
   }
 
   /*
   * Check allocation of a state, allocate if necessary.
   */
   void Sweep::checkAllocation(Sweep::State& state) 
   {
      int nm = mixture().nMonomer();
      int nx = domain().nx();
      UTIL_CHECK(nm > 0);
      UTIL_CHECK(nx > 0);
 
      if (state.isAllocated()) {
         UTIL_CHECK(state.capacity() == nm);
      } else {
         state.allocate(nm);
      }
      for (int j = 0; j < nm; ++j) {
         if (state[j].isAllocated()) {
            UTIL_CHECK(state[j].capacity() == nx);
         } else {
            state[j].allocate(nx);
         }
      }
 
   };
 
   /*
   * Setup operation at beginning sweep.
   *
   * Must call initializeHistory.
   */
   void Sweep::setup() 
   {
      // Initialize history
      initialize();
 
      // Open log summary file
      std::string fileName = baseFileName_;
      fileName += "log";
      fileMaster().openOutputFile(fileName, logFile_);
 
   };
 
   /*
   * Set non-adjustable system parameters to new values.
   *
   * \param s path length coordinate, in range [0,1]
   */
   void Sweep::setParameters(double s) 
   {  
      UTIL_THROW("Called un-implemented R1d::Sweep::setParameters");
   };
 
   /*
   * Create guess for adjustable variables by polynomial extrapolation.
   */
   void Sweep::extrapolate(double sNew) 
   {
      int nm = mixture().nMonomer();
      int nx = domain().nx();
      UTIL_CHECK(nm > 0);
      UTIL_CHECK(nx > 0);
 
      if (historySize() > 1) {
 
         // Compute an array of coefficients for polynomial extrapolation
         setCoefficients(sNew);
 
         // Set System::wField to a linear combination of previous values
         // Note: wField(i) is accessed through SystemAccess base class.
         double coeff;
         int i, j, k;
         for (i = 0; i < nm; ++i) {
            coeff = c(0);
            for (j = 0; j < nx; ++j) {
               wField(i)[j] = coeff*state(0)[i][j];
            }
            for (k = 1; k < historySize(); ++k) {
               coeff = c(k);
               for (j = 0; j < nx; ++j) {
                  wField(i)[j] += coeff*state(k)[i][j];
               }
            }
         }
      }
   };

   int Sweep::solve(bool isContinuation) 
   {  return system().iterate(isContinuation); };

   void Sweep::reset() 
   {  assignFields(wFields(), state(0)); };

   void Sweep::getSolution() 
   {
      // Assign current wFields to state(0)
      assignFields(state(0), wFields());
 
      if (homogeneousMode_ >= 0) {
         comparison_.compute(homogeneousMode_);
      }
 
      // Output solution
      int i = nAccept() - 1;
      std::string fileName = baseFileName_;
      fileName += toString(i);
      outputSolution(fileName);
 
      // Output brief summary to log file
      outputSummary(logFile_);
   };
 
   void Sweep::outputSolution(std::string const & fileName)
   {
      std::ofstream out;
      std::string outFileName;
 
      // Write parameter file, with thermodynamic properties at end
      outFileName = fileName;
      outFileName += ".stt";
      fileMaster().openOutputFile(outFileName, out);
      system().writeParam(out);
      out << std::endl;
      system().writeThermo(out);
 
      if (homogeneousMode_ >= 0) {
         comparison_.output(homogeneousMode_, out);
      }
      out.close();
 
      // Write concentration fields
      outFileName = fileName;
      outFileName += ".c";
      fieldIo_.writeFields(cFields(), outFileName);
 
      // Write chemical potential fields
      outFileName = fileName;
      outFileName += ".w";
      fieldIo_.writeFields(wFields(), outFileName);
   }
 
   void Sweep::outputSummary(std::ostream& out)
   {
      int i = nAccept() - 1;
      double sNew = s(0);
      out << Int(i,5) << Dbl(sNew)
          << Dbl(system().fHelmholtz(),16)
          << Dbl(system().pressure(),16);
      #if 0
      if (homogeneousMode_ != -1) {
         if (homogeneousMode_ == 0) {
            double dF = system().fHelmholtz()
                      - system().homogeneous().fHelmholtz();
            out << Dbl(dF, 16);
         } else {
            double dP = system().pressure()
                      - system().homogeneous().pressure();
            double dOmega = -1.0*dP*domain().volume();
            out << Dbl(dOmega, 16);
         }
      }
      #endif
      out << std::endl;
   }
 
   void Sweep::cleanup()
   {  logFile_.close(); }
 
   void Sweep::assignFields(DArray<System::Field>& lhs,
                            DArray<System::Field> const & rhs) const
   {
 
      int nm = mixture().nMonomer();
      int nx = domain().nx();
 
      UTIL_CHECK(lhs.capacity() == nm);
      UTIL_CHECK(rhs.capacity() == nm);
      int i, j;
      for (i = 0; i < nm; ++i) {
         UTIL_CHECK(rhs[i].isAllocated());
         UTIL_CHECK(rhs[i].capacity() == nx);
         UTIL_CHECK(lhs[i].isAllocated());
         UTIL_CHECK(lhs[i].capacity() == nx);
         for (j = 0; j < nx; ++j) {
            lhs[i][j] = rhs[i][j];
         }
      }
   }
 
} // namespace R1d
} // namespace Pscf