dbpp_Reconstruction.hpp

Go to the documentation of this file.

#pragma once
 
//C+ includes
#include <type_traits>
#include <ranges>
#include <vector>
#include <numeric>
#include <functional>
#include <valarray>
// App include
#include "dbpp_SimulationUtilities.hpp" // MinMod19
// SfxBase19 includes
#include "include/Sfx_StateVector.h"
#include "include/Sfx_StateVariables.h"
#include "include/Sfx_DefineTypes.h"
 
namespace dbpp
{
  template<typename FuncLimiter>
  std::pair<std::valarray<float64>/*dU1*/, std::valarray<float64>/*dU2*/>
    computeDU12( const std::valarray<float64>& aU1,
      const std::valarray<float64>& aU2, FuncLimiter&& aFuncLimtr)
  {
    // Sanity checks
    assert(std::size(aU1) == Sfx::EMCNEILNbSections::value);
    assert(std::size(aU2) == Sfx::EMCNEILNbSections::value);
 
    //
    // Reconstruction process at the cell interface x_j+1/2
    // 
 
    // Need documentation abpout each step of the algorithm 
    std::valarray<float64> dU1(size(aU1));
    std::valarray<float64> dU2(size(aU2));
 
    // compute gradient over each cell
    std::adjacent_difference(std::ranges::begin(aU1), std::ranges::end(aU1),
      std::ranges::begin(dU1));
 
    std::adjacent_difference(std::ranges::begin(aU2), std::ranges::end(aU2),
      std::ranges::begin(dU2));
 
    dU1[0] = 0.; // force first element to 0 (adjacent_difference keep first element unchanged)
    dU2[0] = 0.; // force first element to 0 (adjacent_difference keep first element unchanged)
 
    // apply slope limiter function for each gradient (minimum slope)
    std::adjacent_difference(std::ranges::begin(dU1), std::ranges::end(dU1),  // range to apply limiter
      std::ranges::begin(dU1),                // result 
      aFuncLimtr); // slope limiter function 
 
    std::adjacent_difference(std::ranges::begin(dU2), std::ranges::end(dU2),  // range to apply limiter 
      std::ranges::begin(dU2),                // result 
      aFuncLimtr); // slope limiter function 
 
    auto w_dU1Shifted = dU1.shift(1);  // circular shift (shifted element at the end with element T{})  
    auto w_dU2Shifted = dU2.shift(1);  // circular shift (shifted element at the end with element T{})
 
    // last global discretization node (ghost node) overwrite
    w_dU1Shifted[aU1.size() - 1] = aFuncLimtr(0., aU1[aU1.size() - 1] - aU1[aU1.size() - 2]);  // i+1/2
    w_dU2Shifted[aU2.size() - 1] = aFuncLimtr(0., aU2[aU2.size() - 1] - aU2[aU2.size() - 2]);  // i+1/2
 
    // don't need anymore of shifted array, might as well to move it 
    dU1 = std::move(w_dU1Shifted); // move semantic supported by valarray
    dU2 = std::move(w_dU2Shifted); // ditto
 
    return { dU1,dU2 };
  }

  template<typename F, typename NumArrayType, // state variables components
    typename = std::enable_if<std::is_same_v<NumArrayType, std::valarray<double>>>>
    auto MUSCLReconstruction(const NumArrayType& aU1, const NumArrayType& aU2,
      F&& aSlopeLimiter) // slope limiter function (forward/Universal reference)
  {
    // compute gradient over each cell by applying slope limiter function 
    auto [dU1, dU2] = computeDU12(aU1, aU2, std::forward<F>(aSlopeLimiter));
 
    const NumArrayType w_dU1(std::ranges::begin(dU1), std::ranges::size(w_dU1)); // gradient of the first state variable
    const NumArrayType w_dU2(std::ranges::begin(dU2), std::ranges::size(w_dU2)); // gradient of the second state variable
 
    //  MUSCL (Monotone Upwind Scheme Conservation System Law) 
    //  Second-order linear extrapolation (polynomial reconstruction at the interface).
    //  Total Variation Diminishing (TVD) property guarantees convergence of such schemes
    //  Uses the Minmod limiter for the slope reconstruction.
 
    //
    // Reconstruction process at the cell interface x_j+1/2
    // Second-order achieve by: U + 0.5*dU where dU derivative at first order
    // expand state variable in a serie (dU has been computed by using a function
    // to limit slope called "Slope limiter"), physically similar to adding dissipation or diffusion  
 
    // Compute the left  and right states: UR, UL(state variables) j+1/2
    const NumArrayType UL1{ std::begin(aU1 + 0.5 * w_dU1), Sfx::DIM::value };  // j
    const NumArrayType UL2{ std::begin(aU2 + 0.5 * w_dU2), Sfx::DIM::value };  // j
    const NumArrayType UR1{ std::next(std::begin(aU1 - 0.5 * dU1)), Sfx::DIM::value }; // j+1
    const NumArrayType UR2{ std::next(std::begin(aU2 - 0.5 * dU2)), Sfx::DIM::value }; // j+1
      
    // template argument deduction
    return std::make_tuple(UL1, UL2, UR1, UR2);
   }

   template<typename T, typename CONT, // valid for rvalue only
   typename std::enable_if_t<!std::is_lvalue_reference<T>::value>>
     decltype(auto) MUSCLReconstrv( T&& aA, T&& aQ, CONT aListCellFaces)
   {
     // algo steps
     // --- compute dU
     // --- compute left/right state variables for each faces
     // --- return map (i,UL,UR)
     // 
   //  using namespace emcil;
 
     using mapcellfaceVar = 
       std::map<short/*face idx*/, std::pair<Sfx::StateVector/*UL*/, Sfx::StateVector/*UR*/>>;
     // gradient dU (O(2))
    // std::function<float64(float64, float64)> w_minmodFunc = emcil::HydroUtils::minmod;
     auto w_dA = computeDU(aA, dbpp::MinMod19<float64>{}); // first state variable conserved quantity (A)
     auto w_dQ = computeDU(aQ, dbpp::MinMod19<float64>{}); // second state variable (Q)
 
     mapcellfaceVar w_mapULUR;
     // don't know signature interface of the container, use auto&&
     // no additional copies are made of the values we are iterating through
     for( auto&& cellFace : aListCellFaces)
     { 
       auto w_leftIdxNode  = cellFace.getLeftNodeI();  // stencil [i-1,i]
       auto w_rightIdxNode = cellFace.getRightNodeI(); // stencil [i,i+1]
 
       // C++17 structured binding (bind alias) and try_emplace() if already in, do nothing!!
       auto [pos, succeed] = w_mapULUR.try_emplace(cellFace.getGblFaceI(), // global face index
         Sfx::StateVector{ static_cast<float64>(Sfx::StateVariables{ aA[w_leftIdxNode] + 0.5 * w_dA[w_leftIdxNode]}),
        static_cast<float64>(Sfx::StateVariables{ aQ[w_leftIdxNode] + 0.5 * w_dQ[w_leftIdxNode] }) }, //A
         Sfx::StateVector{ static_cast<float64>(Sfx::StateVariables{ aA[w_rightIdxNode] - 0.5 * w_dA[w_rightIdxNode]}),
         static_cast<float64>(Sfx::StateVariables{ aQ[w_rightIdxNode] - 0.5 * w_dQ[w_rightIdxNode]}) } //Q
       );
   //    if (!succeed) {
    //     pos->
     //  }
         //std::cerr << "Couldn't add the cell face with id: " << 
     }//for-loop
 
     return w_mapULUR;
   }
 
#if 0 //not in use for now
    // Some reconstruction procedure (perfect forwarding)
    // these are templated functions and range can be vector,
    // valarray, numerical array from Elligno math library
    // they all support opertor [] which is a requirement
    //
    //  Reference
    //   Check reconstr files in VS2015 project, lot of functions
    //   with different signatures ... and create a seperate file
    //   for the implementation of these algorithms
    template<typename... Args>
    auto MUSCLReconstrFwd(Args&&... args) // Args&& ...args??
    {
        //std::is_empty<Range>(aRng);
 
        return MUSCLReconstr(std::forward<Args>(args)...);
    }
 
    template<typename Range>
    auto MUSCLReconstr(Range&& aRng1, Range&& aRng2)
    {}
 
    template<typename Range>
    auto MUSCLReconstr(Range& aRng1, Range& aRng2)
    {}
 
    template<typename Range>
    auto MUSCLReconstr(const Range& aRng1, const Range& aRng2)
    {}
 
    template<typename Range>
    auto MUSCLReconstr(const Range& aRng1, const Range& aRng2, Range& aU1LR, Range& aU2LR)
    {}
#endif // 0 //not in use for now
} // End of namespace