dbpp_NujicIntegrator.h

Go to the documentation of this file.

 
#pragma once
 
// C++ include
#include <ranges>
// App includes
#include "../SfxBase/Sfx_Logger.h"
#include "../include/Sfx_Simulation.h"
// SfxBase19 include
#include "include/Sfx_UniversalConstants.h"
// VS2019 includes
#include "../Nov_VS2019/dbpp_SweRhsData.hpp"
#include "../Nov_VS2019/SimulationConfig.h"
 
// forward declaration
namespace dbpp { class SweRhsAlgorithm; }
 
namespace dbpp 
{
  class NujicIntegrator  
  {
    using slnpairvarray = std::pair<std::valarray<float64>, std::valarray<float64>>;
    using slnpairvector = std::pair<std::vector<float64>, std::vector<float64>>;
  public:
    enum class eIntegratorStep {
      predictor=0, 
      corrector=1  
    };
  public:
    NujicIntegrator();
    NujicIntegrator( std::size_t aNbGridPts, float64 aSpatialStep /*DamBreakData*/);
    void step( dbpp::SweRhsData<float64, std::valarray<float64>>&& aRhs, float64 aDt);
    void TwoInOneStep(dbpp::SweRhsAlgorithm* aRhsAlgo, const float64 aDt);
    template<typename Range>
    void initialSolution( const Range& aRng1, const Range& aRng2)
    {
      assert(std::ranges::size(aRng1) == std::ranges::size(aRng2));
 
      if( std::ranges::size(aRng1) == Sfx::DIM::value+1)
      {
        if constexpr( std::is_same_v<Range,std::valarray<float64>>)
        {
          auto&& w_vw1 = std::views::counted(std::ranges::begin(aRng1), std::ranges::size(aRng1) - 1);
          auto&& w_vw2 = std::views::counted(std::ranges::begin(aRng2), std::ranges::size(aRng2) - 1);
          // resize first
          m_currU1Values.resize(w_vw1.size());
          m_currU2Values.resize(w_vw2.size());
          std::ranges::copy(w_vw1, std::ranges::begin(m_currU1Values));
          std::ranges::copy(w_vw2, std::ranges::begin(m_currU2Values));
        }
        else
        {
          auto&& w_vw1 = std::views::counted(aRng1.cbegin(), std::ranges::size(aRng1) - 1);
          auto&& w_vw2 = std::views::counted(aRng2.cbegin(), std::ranges::size(aRng2) - 1);
          // resize first
          m_currU1Values.resize(w_vw1.size());
          m_currU2Values.resize(w_vw2.size());
          std::ranges::copy(w_vw1, std::ranges::begin(m_currU1Values));
          std::ranges::copy(w_vw2, std::ranges::begin(m_currU2Values));
        }
      }
      else if ( std::ranges::size(aRng1) == Sfx::DIM::value)
      {
        // resize first
        m_currU1Values.resize(std::ranges::size(aRng1));
        m_currU2Values.resize(std::ranges::size(aRng2));
        std::ranges::copy(aRng1, std::ranges::begin(m_currU1Values));
        std::ranges::copy(aRng2, std::ranges::begin(m_currU2Values));
      }
      else{
        Sfx::Logger::instance()->OutputError(std::string{ "Not supported container (NujicIntegrator)" }.data());
      }
      allocField();
    }

    void setInitialState( StateVector aField);  // pass by-value (move semantic??)
    StateVector getCurrentState() const
    {
      return m_currState;
    }

    StateVector getMidState() const
    {
       return m_prevState;
    }

    slnpairvarray getSolutionAsVarray() const 
    {
      if( m_integratorStep==eIntegratorStep::predictor)
      {
        return { m_prevU1Values, m_prevU2Values };
      }
      return {m_currU1Values, m_currU2Values}; 
    }

    slnpairvector getSolutionAsVector() const
    {
      if (m_integratorStep == eIntegratorStep::predictor)
      {
        /*auto vec1 = std::vector<float64>(std::begin(m_prevU1Values), std::end(m_prevU1Values));
        auto vec2 = std::vector<float64>(std::begin(m_prevU2Values), std::end(m_prevU2Values));*/
 
        // RVO (Return Value Optimization) copy elison C++17
        return { std::vector<float64>(std::begin(m_prevU1Values), std::end(m_prevU1Values)),
                 std::vector<float64>(std::begin(m_prevU2Values), std::end(m_prevU2Values)) };
      };
 
      // RVO (Return Value Optimization) copy elison C++17
      return { std::vector<float64>(std::begin(m_currU1Values), std::end(m_currU1Values)),
      std::vector<float64>(std::begin(m_currU2Values), std::end(m_currU2Values)) };
    }

    StateVector getSolutionAsField() const noexcept
    {
      if (m_integratorStep == eIntegratorStep::predictor)
      {
        return m_prevState;
      }
      return m_currState;
    }

    void setStep( float64 aStep) noexcept { m_step=aStep;}
    float64 getStep() const noexcept { return m_step;}
    void setIntegratorStep( const eIntegratorStep aIntegratorStep)
    {
      if( m_integratorStep!=aIntegratorStep)
      {
        m_integratorStep = aIntegratorStep;
      }
    }

    eIntegratorStep getIntegratorStep() const noexcept { return m_integratorStep; }
    void setPhysicalBoundaryCnd() {/*under construction*/ }
    void setLeftNodeBCValues( const std::pair<float64, float64>& aPairLeftBC)
    {
      m_leftBCValues = aPairLeftBC;
    }
  protected:
    void predictor( dbpp::SweRhsData<float64,std::valarray<float64>>&& aRhs, float64 aDt);
    void corrector( dbpp::SweRhsData<float64,std::valarray<float64>>&& aRhs, float64 aDt);
 
  private:
    std::size_t m_numOfGridPts;                 
    float64 m_step;                             
    float64 m_spatialStep;                      
    std::pair<float64, float64> m_leftBCValues; 
    eIntegratorStep m_integratorStep;           
    std::valarray<float64> m_prevU1Values;      
    std::valarray<float64> m_prevU2Values;      
    StateVector m_prevState;                    
    std::valarray<float64> m_currU1Values;      
    std::valarray<float64> m_currU2Values;      
    StateVector m_currState;                    
    void allocField();
  };
#if 0
  //reserve must be done before assignment
  template<typename Range>
  void NujicIntegrator::initialSolution( const Range& aRng1, const Range& aRng2)
  {
    // for now we support valarray (under construction)??
    assert(std::ranges::size(aRng1) == std::ranges::size(aRng2));
 
    // concept of viewable range
    if constexpr (std::ranges::view<Range>)
    {
  //    if( std::ranges::size(aRng1) == Sfx::DIM::value + 1)
   //   {
         auto&& w_vw1 = std::views::counted(aRng1.cbegin(), std::ranges::size(aRng1) - 1);
         auto&& w_vw2 = std::views::counted(aRng2.cbegin(), std::ranges::size(aRng2) - 1);
         // resize first
         m_currU1Values.resize(w_vw1.size());
         m_currU2Values.resize(w_vw2.size());
         std::ranges::copy(w_vw1, std::ranges::begin(m_currU1Values));
         std::ranges::copy(w_vw2, std::ranges::begin(m_currU2Values));
  //    }
 //     else // DIM computational domain
 //     {
        // resize first
  //      m_currU1Values.resize(std::ranges::size(aRng1));
  //      m_currU2Values.resize(std::ranges::size(aRng2));
  //     std::ranges::copy(aRng1, std::ranges::begin(m_currU1Values));
  //      std::ranges::copy(aRng2, std::ranges::begin(m_currU2Values));
   //   }
    }
    else if constexpr (std::is_same_v<Range, std::valarray<float64>> )
    {
      m_currU1Values = aRng1; //move?? No because 
      m_currU2Values = aRng2;
    }
    else
    {
      std::copy(aRng1.cbegin(), aRng1.cend(), std::begin(m_currU1Values));
      std::copy(aRng2.cbegin(), aRng2.cend(), std::begin(m_currU2Values));
    }
#endif
 
} // end of namespace