dbpp_HLLParametersAlgo.h

Go to the documentation of this file.

#pragma once
 
// SfxBase19 includes
#include "include/Sfx_UniversalConstants.h"
#include "include/Sfx_DefineTypes.h"
#include "include/Sfx_StateVector.h"
// DamBreak include
//#include "../../Nov_VS2019/dbpp_HLLParameterStruct.h"
 
namespace dbpp 
{
    struct HLLParameterStruct
    {
        float64 uL;    
        float64 uR;    
        float64 CL;  
        float64 CR;  
    };
    HLLParameterStruct computeHLLParameters( const Sfx::StateVector& aUL,
        const Sfx::StateVector& aUR, float64 aSectionWidth = 1.)
    {
        // celerity wave on the left and right (unit width (=1))
        // In the case of a flat bed A=H (water level) but when Z!=0 
        // A = ?? water depth (not sure but need to think about that)
        auto CR = std::sqrt(Sfx::cGravity<float64>*aUR.A())/aSectionWidth;
        auto CL = std::sqrt(Sfx::cGravity<float64>*aUL.A())/aSectionWidth;
    
        // uL and uR are velocities on the left and right states
        auto uL = aUL.Q() / aUL.A();
        auto uR = aUR.Q() / aUR.A();
 
        return HLLParameterStruct{uL,uR,CL,CR};
    }

    std::pair<float64, float64> computeHLLShockSpeed( HLLParameterStruct aPrmStruct)
    {
        // IMPORTANT sign error in these formula below (need to be fix!!!)
        //  see article (J. Belanger Technical Report 2008) 
        // "Validating Shock Capturing Schemes On The DamBreak Problem") 
        // Also other references:
        //  Alcrudo F. "Advanced Mathematical Modeling Techniques 
        //  For FloodPropagation in Natural Topographies", pp. 15
        //  Xinya Y., Sam S. Y. Wang
        //  "Improved HLL Scheme for 1D Dam-Break Flows Over Complex Topography" pp. 34
        // E. McNeil Note (original code)
    //  Une attention particulière devra être ici portée sur
      //    les signes (+/-) relativement au produit scalaire avec Si+1/2
        auto CS = 0.5*( aPrmStruct.CL + aPrmStruct.CR) - 0.25*(aPrmStruct.uL - aPrmStruct.uR); // sign error to be fix (... '+' 0.25*())
        auto uS = 0.5*( aPrmStruct.uL - aPrmStruct.uR) + aPrmStruct.CL - aPrmStruct.CR;        // sign error to be fix 0.5*(uL '+' uR) + ...
        // shock speed at left/right face
        auto SL = std::min( aPrmStruct.uL - aPrmStruct.CL, uS - CS);
        auto SR = std::max( aPrmStruct.uR + aPrmStruct.CR, uS + CS);
 
        return{ SL, SR };
    }
    void computeHLLPhysicalFlux(const Sfx::StateVector& aUL,
        const Sfx::StateVector& aUR, float64 aSectionWidth = 1.)
    {
        // retrieve conservative system   
        // Factory to create math equation  
    //  const auto w_stVen1D = aRprob.getMthEquations();
//      StVenant1D w_stVen1D;
 
        // compute state variables physical flux at cell face
        //const auto FL1 = aRprob.getUL().Q();
        //const auto FR1 = aRprob.getUR().Q();
 
        // physical flux (math model in use)
        //const auto FL2 = w_stVen1D->flux(Sfx::StateVector{ aRprob.getUL().A(), aRprob.getUL().Q() }); //FL2
        //const auto FR2 = w_stVen1D->flux(Sfx::StateVector{ aRprob.getUR().A(), aRprob.getUR().Q() }); //FR2
    }
} //End of namespace