rpg/solvers/Block.h

#ifndef RPG_BLOCK_H
#define RPG_BLOCK_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 <pscf/solvers/BlockTmpl.h>      // base class template
#include "Propagator.h"                  // base class argument
 
#include <prdc/cuda/RField.h>            // member
#include <prdc/cuda/RFieldDft.h>         // member
#include <prdc/cuda/FFTBatched.h>        // member
#include <pscf/cuda/DeviceArray.h>       // member
#include <util/containers/FSArray.h>     // member
 
#include <prdc/cuda/resources.h>        
 
// Forward declarations
namespace Pscf {
   template <int D> class Mesh;
   namespace Prdc {
      template <int D> class UnitCell;
      namespace Cuda {
         template <int D> class WaveList;
         template <int D> class FFT;
      }
   }
   namespace Rpg {
      template <int D> class FieldIo;
   }
}
 
namespace Pscf {
namespace Rpg {
 
 
   using namespace Util;
   using namespace Pscf::Prdc;
   using namespace Pscf::Prdc::Cuda;

   template <int D>
   class Block : public BlockTmpl< Propagator<D> >
   {
 
   public:
 
      // Public type name aliases

      using Base = BlockTmpl< Propagator<D> > ;

      using PropagatorT = Propagator<D>;

      using FFTT = FFT<D>;

      using WaveListT = WaveList<D>;

      using FieldIoT = FieldIo<D>;
 
      // Public member functions

      Block();

      ~Block();

      void associate(Mesh<D> const & mesh, 
                     FFT<D> const & fft, 
                     UnitCell<D> const & cell, 
                     WaveList<D>& wavelist);

      void allocate(double ds, bool useBatchedFFT = true);

      void clearUnitCellData();

      void setLength(double newLength);

      void setKuhn(double kuhn);

      void setupSolver(RField<D> const & w);

      void stepThread(RField<D> const & qin, RField<D>& qout);

      void stepBead(RField<D> const & qin, RField<D>& qout);

      void stepFieldBead(RField<D> & q);

      void stepBondBead(RField<D> const & qin, RField<D>& qout);
 
      void stepHalfBondBead(RField<D> const & qin, RField<D>& qout);

      void computeConcentrationThread(double prefactor);

      void computeConcentrationBead(double prefactor);

      double averageProduct(RField<D> const& q0, RField<D> const& q1);

      double averageProductBead(RField<D> const& q0, RField<D> const& q1);

      void computeStressThread(double prefactor);

      void computeStressBead(double prefactor);

      double stress(int n);

      Mesh<D> const & mesh() const;

      FFT<D> const & fft() const;

      double ds() const;

      int ns() const;
 
      // Functions with non-dependent names from BlockTmpl<Propagator<D>>
      using BlockTmpl< Pscf::Rpg::Propagator<D> >::setKuhn;
      using BlockTmpl< Pscf::Rpg::Propagator<D> >::propagator;
      using BlockTmpl< Pscf::Rpg::Propagator<D> >::cField;
      using BlockTmpl< Pscf::Rpg::Propagator<D> >::kuhn;
 
      // Functions with non-dependent names from Edge
      using Edge::setId;
      using Edge::setVertexIds;
      using Edge::setMonomerId;
      using Edge::setLength;
      using Edge::id;
      using Edge::monomerId;
      using Edge::vertexIds;
      using Edge::vertexId;
      using Edge::length;
      using Edge::nBead;
 
   private:

      FFTBatched<D> fftBatchedAll_;

      FFTBatched<D> fftBatchedPair_;

      FSArray<double, 6> stress_;

      RField<D> expKsq_;

      RField<D> expW_;

      RField<D> expKsq2_;

      RField<D> expW2_;

      RField<D> expWInv_;

      DeviceArray<cudaReal> qrPair_;

      DeviceArray<cudaComplex> qkPair_;
 
      // R-grid work space (used in productAverage)
      RField<D> qr_; 
 
      // K-grid work space (used for FFT of q in stepBondBead)
      RFieldDft<D> qk_; 

      DeviceArray<cudaComplex> q0kBatched_;

      DeviceArray<cudaComplex> q1kBatched_;
 
      // Slices of forward and reverse propagator on a k-grid (for stress)
      RFieldDft<D> q0k_; 
      RFieldDft<D> q1k_;

      Mesh<D> const * meshPtr_;

      FFT<D> const * fftPtr_;

      UnitCell<D> const * unitCellPtr_;

      WaveList<D> * waveListPtr_;

      IntVec<D> kMeshDimensions_;

      int kSize_;

      double ds_;

      double dsTarget_;

      int ns_;

      bool isAllocated_;

      bool hasExpKsq_;

      bool useBatchedFFT_;

      UnitCell<D> const & unitCell() const
      {  return *unitCellPtr_; }

      WaveList<D> const & wavelist() const
      {  return *waveListPtr_; }

      int nParams_;

      void computeExpKsq();
 
   };
 
   // Inline member functions
 
   // Get number of contour steps.
   template <int D>
   inline int Block<D>::ns() const
   {  return ns_; }
 
   // Get contour length step size.
   template <int D>
   inline double Block<D>::ds() const
   {  return ds_; }
 
   // Get derivative of free energy w/ respect to a unit cell parameter.
   template <int D>
   inline double Block<D>::stress(int n)
   {  return stress_[n]; }
 
   // Get Mesh by reference.
   template <int D>
   inline Mesh<D> const & Block<D>::mesh() const
   {
      UTIL_ASSERT(meshPtr_);
      return *meshPtr_;
   }
 
   // Get FFT by reference.
   template <int D>
   inline FFT<D> const & Block<D>::fft() const
   {
      UTIL_ASSERT(fftPtr_);
      return *fftPtr_;
   }
 
   #ifndef RPG_BLOCK_TPP
   // Suppress implicit instantiation
   extern template class Block<1>;
   extern template class Block<2>;
   extern template class Block<3>;
   #endif
 
}
}
#endif