pspg/solvers/Block.h

#ifndef PSPG_BLOCK_H
#define PSPG_BLOCK_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.
*/
 
#include "Propagator.h"                   // base class argument
#include <pscf/solvers/BlockTmpl.h>       // base class template
#include <pspg/field/RDField.h>           // member
#include <pspg/field/RDFieldDft.h>        // member
#include <pspg/field/FFT.h>               // member
#include <pspg/field/FFTBatched.h>        // member
#include <util/containers/FArray.h>
#include <pscf/crystal/UnitCell.h>
#include <pspg/solvers/WaveList.h>
 
 
namespace Pscf { 
 
   template <int D> class Mesh; 
 
namespace Pspg { 
 
   using namespace Util;
 
   template <int D>
   class Block : public BlockTmpl< Propagator<D> >
   {
 
   public:
 
      Block();
 
      ~Block();
 
      void setDiscretization(double ds, const Mesh<D>& mesh);
 
      void setupUnitCell(UnitCell<D> const & unitCell, 
                         WaveList<D> const & waveList);
 
      void setLength(double length);
 
      void setKuhn(double kuhn);
 
      void setupSolver(RDField<D> const & w);
 
      void setupFFT();
 
      void step(cudaReal const * q, cudaReal* qNew);
 
      void computeConcentration(double prefactor);
 
      void computeStress(WaveList<D> const & waveList, double prefactor);
 
      double stress(int n);
 
      Mesh<D> const & mesh() const;
 
      double ds() const;
 
      int ns() const;
 
      // Functions with non-dependent names from BlockTmpl<Propagator<D>>
      using BlockTmpl< Pscf::Pspg::Propagator<D> >::setKuhn;
      using BlockTmpl< Pscf::Pspg::Propagator<D> >::propagator;
      using BlockTmpl< Pscf::Pspg::Propagator<D> >::cField;
      using BlockTmpl< Pscf::Pspg::Propagator<D> >::length;
      using BlockTmpl< Pscf::Pspg::Propagator<D> >::kuhn;
 
      // Functions with non-dependent names from BlockDescriptor
      using BlockDescriptor::setId;
      using BlockDescriptor::setVertexIds;
      using BlockDescriptor::setMonomerId;
      using BlockDescriptor::setLength;
      using BlockDescriptor::id;
      using BlockDescriptor::monomerId;
      using BlockDescriptor::vertexIds;
      using BlockDescriptor::vertexId;
      using BlockDescriptor::length;
 
   private:
 
      int nBlocks_;
 
      int nThreads_;
 
      FFT<D> fft_;
 
      FFTBatched<D> fftBatched_;
      
      FArray<double, 6> stress_;
      
      // Array of elements containing exp(-K^2 b^2 ds/6) on k-grid
      RDField<D> expKsq_;
      
      // Array of elements containing exp(-K^2 b^2 ds/12) on k-grid
      RDField<D> expKsq2_;
 
      // Array of elements containing exp(-W[i] ds/2) on r-grid
      RDField<D> expW_;
 
      // Array of elements containing exp(-W[i] ds/4) on r-grid
      RDField<D> expW2_; 
 
      // Work arrays for r-grid fields
      RDField<D> qr_;
      RDField<D> qr2_;
 
      // Work arrays for wavevector space (k-grid) field
      RDFieldDft<D> qk_;
      RDFieldDft<D> qk2_;
 
      // Batched FFTs of q
      cudaComplex* qkBatched_;
      cudaComplex* qk2Batched_;
 
      // Propagators on r-grid
      RDField<D> q1_;
      RDField<D> q2_;
 
      cudaReal* d_temp_;
      cudaReal* temp_;
 
      cudaReal* expKsq_host;
      cudaReal* expKsq2_host;
 
      Mesh<D> const * meshPtr_;
 
      UnitCell<D> const* unitCellPtr_;
 
      WaveList<D> const * waveListPtr_;
      
      IntVec<D> kMeshDimensions_;
 
      int kSize_;
 
      double ds_;
 
      int ns_;
 
      bool isAllocated_;
 
      bool hasExpKsq_;
 
      UnitCell<D> const & unitCell() const 
      {  return *unitCellPtr_; }
 
      WaveList<D> const & wavelist() const 
      {  return *waveListPtr_; }
 
      int nParams_;
 
      void computeExpKsq();
 
   };
 
   // Inline member functions
 
   template <int D>
   inline int Block<D>::ns() const
   {  return ns_; }
 
   template <int D>
   inline double Block<D>::ds() const
   {  return ds_; }
 
   template <int D>
   inline double Block<D>::stress(int n)
   {  return stress_[n]; }
 
   template <int D>
   inline Mesh<D> const & Block<D>::mesh() const
   {   
      UTIL_ASSERT(meshPtr_);
      return *meshPtr_;
   }
 
   #ifndef PSPG_BLOCK_TPP
   // Suppress implicit instantiation
   extern template class Block<1>;
   extern template class Block<2>;
   extern template class Block<3>;
   #endif
 
}
}
//#include "Block.tpp"
#endif