dbpp_MUSCLReconstr.hpp

Go to the documentation of this file.

#pragma once
 
// C++ includes
#include <map>
#include <list>
#include <numeric>
#include <functional>
// C++20 include
#include <ranges>
// SfxBase19 inlcudes
#include "include/Sfx_StateVector.h"
#include "include/Sfx_StateVariables.h"
#include "include/Sfx_cellFaceVariables.h"
// DamBreak++ includes
#include "../SfxBase/Sfx_MathFunctions.hpp" // computeDU
// VS15 include
//#include "../May_VS2015/Test/GlobalDiscretization.h"
// 
//#include "dbpp_Omega.h"
//#include "dbpp_CellFace.h"
#include "dbpp_Enumerations.h"
//#include "dbpp_StateVectorField.h"
#include "dbpp_SimulationUtilities.hpp"
#include "dbpp_FiniteVolumeDiscretization.h"
 
//namespace emcil { class GlobalDiscretization; }
 
namespace dbpp 
{
  // Design Note
  //  could be implemented as template, not sure yet
     // Create a list of global cell face (deprecated used for prototyping) according to grid
     // Reconstruction of state variables at global cell - face j+1/2
     //const auto w_gblFaceVar = reconstr(vectorField, GlobalDiscretization);

    class MUSCLReconstr
    {
        using CellId = short;
        //using LRvariables = std::pair<CellFaceVariables/*A*/, CellFaceVariables/*Q*/>;
        //using mapCellFaceULR = std::map<CellId, LRvariables>;
        using slopeLimiterFunc = std::function<float64(float64, float64)>;
        //using pairdblvec = std::pair<std::vector<float64>, std::vector<float64>>;
        using mapcellfaceVar = std::map<short/*face idx*/, std::pair<Sfx::StateVector/*UL*/, Sfx::StateVector/*UR*/>>;
 
    public:
        //default 
        MUSCLReconstr() = default;
 
        // ctor (under construction) ???
        MUSCLReconstr(slopeLimiterFunc aFunction, unsigned aOrder)
      : m_dUOrder{1},
            m_faceVarOrder{static_cast<short>(aOrder)} // need to check 
        {}
 
        void setFaceVariableOrder(short aOrder) { m_faceVarOrder = aOrder; }
 
        // Reconstruction of state variables at cell face by using slope limiter and TVD method
        void setSlopeLimiter(slopeLimiterFunc aSlopeLimiterFunc/*HydroUtils::minmod*/) {} // slope limiter gradient 
        void setDUOrder(short aDUOrder) { m_dUOrder = aDUOrder; } // gradient over each cell at 1st order
        short getDUOrder() const { return m_dUOrder; }
//     void computeDU1stOrder( const emcil::vectorField& aVec, 
//       const emcil::GlobalDiscretization* aGblDiscr); // gradient 
        void setReconstrVariableOrder(short aVarOrder) { m_faceVarOrder = aVarOrder; }
        short getReconstrVariableOrder() const { return m_faceVarOrder; }
        void setPhysicalBndCnd() {} // to be implemented!!!
             // enum class eFluxType {incomplete, complete};    
    public:
 
        template<std::ranges::contiguous_range Range> // must be DIM+1 size
        std::list<Sfx::cellFaceVariables> reconstr( const dbpp::Omega& aOmega,
            Range&& aA, Range&& aQ)  // nodal values
        {
            // 'DIM' metaprogramming constant
            assert(std::ranges::size(aA) == Sfx::DIM::value + 1);
            assert(std::ranges::size(aA)==std::ranges::size(aQ));
 
            // check value category
            
            // retrieve geometric node (global domain)
            const auto& w_gnode = aOmega.omega_gnode();
            using enum dbpp::eNodalField;
 
            // check value type
            static_assert(std::is_same_v<Range, std::vector<float64>>);
 
            // support to different version of 'Range'. Vector has a push_back
            // ... to be completed
            if constexpr (!std::is_same_v<Range, std::vector<float64>>)
            {
        // Apply physical boundary cnd (half-open domain)
        if (auto&& node = w_gnode.front(); node.isTiedNode())
        {
          //state variable set
          aA[node.nodeNo()] = std::get<dbpp::toEnumType(A)>(m_upstreamBC.Values());
          aQ[node.nodeNo()] = std::get<dbpp::toEnumType(Q)>(m_upstreamBC.Values());
        }
        if (auto&& node = w_gnode.back(); node.isGhostNode())
        {
          aA[node.nodeNo()] = std::get<dbpp::toEnumType(A)>(m_downstreamBC.Values());
          aQ[node.nodeNo()] = std::get<dbpp::toEnumType(Q)>(m_downstreamBC.Values());
        }
            }
 
            // algo steps
      // --- compute dU
      // --- compute left/right state variables for each faces
      // --- return map (i,UL,UR)
      // 
      // gradient dU (O(2))
            //std::function<float64(float64, float64)> w_minmodFunc = MinMod19;
            auto w_dA = Sfx::computeDU(aA, MinMod19<float64>{}); // first state variable conserved quantity (A)
            auto w_dQ = Sfx::computeDU(aQ, MinMod19<float64>{}); // second state variable (Q)
 
        //  mapcellfaceVar w_mapULUR;
            std::list<Sfx::cellFaceVariables> w_listReconstrVar;
            
            // Create a list of global cell face (deprecated used for prototyping)
            auto w_omegaGlbFaces = aOmega.getGblFaces();
            auto begListFaces = w_omegaGlbFaces.begin();
            //auto i = 0;
            for( const auto& w_gnode : aOmega.omega_gnode())
            {
            //  ++i;
                //tied node and ghost (no part of computational domain)
                if( w_gnode.isTiedNode() && !w_gnode.isGhostNode()) // first node is tied
                {
                    //++begListFaces;  next in the list
                    continue;
                }
 
                // 'cellFace' type is a representation of the cell interface (j+1/2)
                const cellFace& w_cFace = *begListFaces++;
                auto w_leftIdxNode = w_cFace.getLeftNodeI();    // stencil [i-1,i]
                auto w_rightIdxNode = w_cFace.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(w_cFace.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
                //);
                // rvalue refernce version of push_back
                w_listReconstrVar.push_back( Sfx::cellFaceVariables(w_cFace.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 state variable
                    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 state variable
                ));
            }//while-loop
 
            return w_listReconstrVar;
        }
        
        // boundary condition  
        void setPhysicalBC( const dbpp::PhysicalBoundaryCnd& aPhysBC) 
        { 
            m_upstreamBC = aPhysBC.getLeftEnd();
            m_downstreamBC = aPhysBC.getRightEnd();
        }
 
        // reconstruct a scalar field to the given order
        // MUSCL reconstruction algorithm (Monolitic Upstream System Conservation Law)
        virtual std::vector<Sfx::cellFaceVariables> reconstr( const Sfx::StateVectorField& aU,
            const std::shared_ptr<FiniteVolumeDiscretization>& aGblDiscr)
        {
            //auto U1vec = aU.getAField().asStdVector();
            //auto U2vec = aU.Q().values().to_stdVector();
            std::vector<float64> U1vec(std::ranges::begin(aU.Avalues()), std::ranges::end(aU.Avalues()));
            std::vector<float64> U2vec(std::ranges::begin(aU.Qvalues()), std::ranges::end(aU.Qvalues()));
 
        //  using enum eNodalVarComp; //C++20 using enum declaration
            using enum eNodalField; // scoped enum
 
            // testing what exactly?
            const auto& w_uh = aGblDiscr->U_h();
            const auto& arrayNodal = w_uh.vecNodalVarArray().front();
            //auto checkNo = arrayNodal.node_no(); // debug purpose
 
            // apply B.C. upstream state variables
            // implicit conversion from 'int' to 'size_t' (std::get<> take size_t as arg)  
            U1vec[0] = std::get<toEnumType(A)>(m_upstreamBC.m_tpl);// A
            U2vec[0] = std::get<toEnumType(Q)>(m_upstreamBC.m_tpl); // Q
 
            // ReconstrUtility.h version
            const auto& dU1vec = Sfx::computeDU(U1vec, MinMod19<float64>{});
            const auto& dU2vec = Sfx::computeDU(U2vec, MinMod19<float64>{});
 
            // ... to be completed
            std::vector<Sfx::cellFaceVariables> w_retVarReconstr;
            // main steps of the algorithm
            w_retVarReconstr.reserve(aGblDiscr->Omega_e().getGblFaces().size());
 
            // Compute cell face flux j+1/2 (global faces) 
            // solve for face of the cell and then compute dF = F_i+1 - F_i (bias upwind)  
            for (const auto& w_cellFace : aGblDiscr->Omega_e().getGblFaces())   //emcil::createListGlobalFaces(aGblDiscr->omega().omega_eh_array().size()))
            {
                // left face 'L' (x_i+1/2 global index) U_i + 0.5*dU_i
                auto U1L = U1vec[w_cellFace.getLeftNodeI()] + 0.5 * dU1vec[w_cellFace.getLeftNodeI()];
                auto U2L = U2vec[w_cellFace.getLeftNodeI()] + 0.5 * dU2vec[w_cellFace.getLeftNodeI()];
 
                // right face 'R' (x_i+1/2 global index) U_i+1 - 0.5*dU_i+1
                auto U1R = U1vec[w_cellFace.getRightNodeI()] - 0.5 * dU1vec[w_cellFace.getRightNodeI()];
                auto U2R = U2vec[w_cellFace.getRightNodeI()] - 0.5 * dU2vec[w_cellFace.getRightNodeI()];
 
                // psuh_back move version (pass a 'rvalue' reference)
                w_retVarReconstr.push_back(Sfx::cellFaceVariables{ static_cast<unsigned>(w_cellFace.getGblFaceI()),
                    Sfx::StateVector{ U1L,U2L }, Sfx::StateVector{ U1R,U2R } });
            }//for-loop  
 
            return w_retVarReconstr; //return by copy
        }
    private:
        short m_dUOrder{1};                
        short m_faceVarOrder{2};           
        PhyBCNdlConstraint m_upstreamBC;   
        PhyBCNdlConstraint m_downstreamBC; 
    };
} // End of namespace