dbpp_FluxAlgorithmFunctions.hpp

Go to the documentation of this file.

#pragma once
 
#include <valarray>
 
#include "include/Sfx_UniversalConstants.h"
 
namespace dbpp 
{
  template<typename T, typename ArrayType = std::valarray<T>>
  void hllFluxAlgorithm( ArrayType& aFF1, ArrayType& aFF2,  // return 
    const ArrayType& aU1, const ArrayType& aU2,      // in 
    ArrayType&& dU1, ArrayType&& dU2) noexcept // (forward reference)
  {
    std::cout << "Debugging purpose\n";
 
    static_assert(std::is_same_v<decltype(ArrayType), const std::valarray<double, 3>&>);
    static_assert(std::size(aU1) == std::size(aU2), "Numerical array must be same size");
 
    // return global domain faces
    //auto w_glbFaces = createGlobalFaces();
 
    for( auto i = 0; i < std::size(aU1) - 1; ++i)
    {
      // U1 variable
      const auto w_U1L = aU1[i] + dU1[i];
      const auto w_U1R = aU1[i + 1] - dU1[i + 1];
 
      // U2 variable
      const auto w_U2L = aU2[i] + dU2[i];
      const auto w_U2R = aU2[i + 1] + dU2[i + 1];
      /*     StateVector UL{ U1[i] + dU1[i], U2[i] + dU2[i] };
           StateVector UR{ U1[i + 1] - dU1[i + 1], U2[i + 1] - dU2[i + 1] };*/
           //  std::cout << "StateVector not supported yet\n";
             // compute physical flux
             // aFF1[i]=1.22;
             // aFF2[i]=0/23;
    }
  }

  template< typename ArrayType,
    auto NbSections = Sfx::EMCNEILNbSections::value>
    void hllFluxAlgorithm(ArrayType& aFF1, ArrayType& aFF2, // face flux
      const  ArrayType& aU1, const ArrayType& aU2) noexcept // ?? not sure
  {
    // C++20 range concept and compile-time if  (if with initialization and follow by condition C++17)
    if constexpr( range_t<ArrayType> numrng; !std::ranges::range<decltype(numrng)>) // condition if it is a range
    {
      static_assert(std::size(ArrayType) == NbSections);
      std::cout << "numrng is not a range with exact value\n";
    }
 
    // std::ranges::distance(); just check
    using sect_sizet = decltype(NbSections);
    using arr_valuetype = typename ArrayType::value_type;
 
    // distance
    auto rngFF1Dist = std::ranges::distance(aFF1);
 
    // aU over global domain (include ghost node)
  //  static_assert( std::size(aFF1) == NbSections - 1);         //<< "Face Flux and global U array incompatible";
  //  static_assert( std::size(aFF2) == NbSections - 1);         //<< "Face Flux and global U array incompatible";
    static_assert(std::is_same_v<double, ArrayType::value_type>); //<< "";
 
#if 0
            // slope limiter stencil: 
    // Do some reconstruction of state variables at cell face (MUSCLreconstr)
    ArrayType dU1{ aU1.size() }; // gradient over cell
    std::adjacent_difference(std::begin(aU1), std::end(aU1), std::begin(dU1));
    dU1[0] = 0.; // force to zero
    ArrayType dU2{ aU2.size() }; // gradient over cell 
    std::adjacent_difference(std::begin(aU2), std::end(aU2), std::begin(dU2));
    dU2[0] = 0.; // force to zero (adjacent_difference leave first element as is)
 
    // Design Note
    //  we need to use range view for this slope limiter because
    //  compare adjacent value with stencil [i,i+1]
    //  [0,N[ and ]0,N]
    // slope limiter
    MinMod w_minModSLopeLimiter{};
    std::adjacent_difference(dU1.cbegin(), dU1.cend(), dU2.cbegin(), dU1.begin(), w_minModSLopeLimiter);
 
    // ... not sure about this one
    std::vector<arr_valuetype> w_dU1(std::begin(dU1), std::end(dU1));
    w_dU1.push_back(0.);
 
    // limiter function gradient sharp profile
    // remove first element ( boost range ]begin(),end()]) 
 
    std::vector<typename ArrayType::value_type> w_dU2(std::begin(dU2), std::end(dU2));
    w_dU2.push_back(0.);
    // limiter function gradient sharp profile
    // remove first element ( boost range ]begin(),end()]) 
 
#endif // 0
 
    // loop on cell faces
 /*   auto UL = 0.;
    auto UR = 0.;
 
    auto w_glbFaces = createGlobalFaces();
 
    for( auto i = 0; i < w_glbFaces.length(); ++i)
    {
      StateVector UL{ U1[i] + dU1[i], U2[i] + dU2[i] };
      StateVector UR{ U1[i+1] - dU1[i+1], U2[i+1] - dU2[i+1] };*/
 
      // compute physical flux
  }
} // End of namespace