rpg/solvers/Block.h

#ifndef RPG_BLOCK_H
#define RPG_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 <prdc/cuda/RField.h>            // member
#include <prdc/cuda/RFieldDft.h>         // member
#include <prdc/cuda/FFT.h>               // member
#include <prdc/cuda/FFTBatched.h>        // member
#include <rpg/solvers/WaveList.h>
#include <prdc/crystal/UnitCell.h>
#include <pscf/solvers/BlockTmpl.h>      // base class template
#include <util/containers/FSArray.h>
 
namespace Pscf {
 
   template <int D> class Mesh;
 
namespace Rpg {
 
   using namespace Util;
   using namespace Pscf::Prdc;
   using namespace Pscf::Prdc::Cuda;
 
   // Wrappers for the CUDA kernels used internally
   
   void realMulVConjVV(DeviceArray<cudaReal>& a, 
                       DeviceArray<cudaComplex> const & b,
                       DeviceArray<cudaComplex> const & c,
                       DeviceArray<cudaReal> const & d);
   
   void richardsonEx(DeviceArray<cudaReal>& qNew, 
                     DeviceArray<cudaReal> const & qr,
                     DeviceArray<cudaReal> const & qr2, 
                     DeviceArray<cudaReal> const & expW2);
   
   void addEqMulVVc(DeviceArray<cudaReal>& a, DeviceArray<cudaReal> const & b,
                    DeviceArray<cudaReal> const & c, cudaReal const d);
 
   template <int D>
   class Block : public BlockTmpl< Propagator<D> >
   {
 
   public:
 
      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 step(RField<D> const & q, RField<D>& qNew);
 
      void computeConcentration(double prefactor);
 
      void computeStress(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> >::length;
      using BlockTmpl< Pscf::Rpg::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:
 
      FFTBatched<D> fftBatchedAll_;
 
      FFTBatched<D> fftBatchedPair_;
 
      FSArray<double, 6> stress_;
 
      RField<D> expKsq_;
 
      RField<D> expKsq2_;
 
      RField<D> expW_;
 
      RField<D> expW2_;
 
      DeviceArray<cudaReal> qrPair_;
 
      DeviceArray<cudaComplex> qkPair_;
 
      DeviceArray<cudaComplex> qkBatched_;
 
      DeviceArray<cudaComplex> qk2Batched_;
 
      // Propagators on r-grid
      RField<D> q1_;
      RField<D> q2_;
 
      HostDArray<cudaReal> expKsq_h_;
 
      HostDArray<cudaReal> expKsq2_h_;
 
      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