Sfx_MathFunctions.hpp

Go to the documentation of this file.

#pragma once
 
// C++ includes
#include <cassert>
#include <type_traits> 
#include <ranges> // C++20 
// numeric includes
#include <numeric>
#include <valarray>
// SfxBase19 includes
#include "include/Sfx_Utility.h"
#include "include/BaseMacros.h"
#include "include/Sfx_UniversalConstants.h"
// App include
#include "../Nov_VS2019/dbpp_SimulationUtilities.hpp"
 
namespace Sfx 
{
    std::valarray<float64> cell_face_average(std::ranges::view auto aRng)
    {
        // C++20 function declaration
        //auto viewBeg = std::ranges::cbegin(aRng);
        //auto viewSiz = std::ranges::size(aRng);
        assert( Sfx::EMCNEILNbSections::value == std::ranges::size(aRng));
 
        std::valarray<double> w_valArray(Sfx::DIM::value); // initialize to 0
        // Based on Nujic paper (1995)    
        for( auto i = 0; i < Sfx::DIM::value; ++i)
        {
            if (i == 0) // boundaary point
            {
                w_valArray[i] = (0.5 *(aRng[i + 1] + aRng[i]));
            }
            else
            {
                w_valArray[i] = (0.5 * (aRng[i + 1] + aRng[i - 1]));
            }
        }
        return w_valArray;
    }
 
        // Alias template
        //  not valid
        //template<typename Range1, typename Range2>
        //using EnableIfIsameRngType = std::is_same_t<Range1, Range2>;
 
        /*template<typename Range1, typename Range2>
        using EnableIfIsameRngSiz = std::is_same_v<std::size(Range1), std::size(Range2)>;*/
 
        // DESIGN NOTE
//   can be done with the STL Adjacent_difference algorithm.
// I have an implementation in the VS15 project, should it. 
// Delta function for the positive part computed by the  
// ENO stencil which is a three point stencil [j-1,j,j+1]

        template<typename Range1, typename Range2>
        std::common_type_t<Range1, Range2> computeLimitedGradient( const Range1& aRng1, const Range2& aRng2)
        {
            // some basic check before proceeding (sanity check)
            static_assert( dbpp::isHomogeneous(aRng1,aRng2));
            static_assert( std::is_same_v<decltype(aRng1), const std::valarray<double>&>, "Not using numeric valarray");
            static_assert( std::is_same_v<decltype(aRng2), const std::valarray<double>&>, "Not using numeric valarray");
 
            if constexpr (std::ranges::size(aRng1) != std::ranges::size(aRng2))
            {
                return { Range1 {}, Range2 {} };
            }
            else if constexpr (std::ranges::empty(aRng1) && std::ranges::empty(aRng2))
            {
                return { Range1 {}, Range2 {} };
            }
            else
            {
                // gradient according to limiter function
                Range1 w_df1p{ aRng1.size() };
                Range2 w_df2p{ aRng2.size() };
 
                // boundary node i=0
                w_df1p[0] = Sfx::minmod(aRng1[1] - aRng1[0], 0.);
                w_df2p[0] = Sfx::minmod(aRng2[1] - aRng2[0], 0.);
 
                // interior node (exclude i=0 and i=100 global node index)
                for (uint32 i2 = 1; i2 < aRng1.size() - 1; ++i2)
                {
                    w_df1p[i2] = Sfx::minmod(aRng1[i2 + 1] - aRng1[i2], aRng1[i2] - aRng1[i2 - 1]);
                    w_df2p[i2] = Sfx::minmod(aRng2[i2 + 1] - aRng2[i2], aRng2[i2] - aRng2[i2 - 1]);
                }
 
                // boundary cnd i=100 ghost node (avoid index out of range)
                w_df1p[aRng1.size() - 1] = Sfx::minmod(0., aRng1[aRng1.size() - 1] - aRng1[aRng1.size() - 2]);
                w_df2p[aRng2.size() - 1] = Sfx::minmod(0., aRng2[aRng2.size() - 1] - aRng2[aRng2.size() - 2]);
 
                return { w_df1p, w_df2p };
            }
        }
        
        template<typename Range, typename FuncLimiter>        // let compiler deduce type (moved to "Dsn_MathUtils.hpp")
        requires std::same_as<Range, std::vector<float64>>
        auto computeDU(Range aRng, FuncLimiter&& aFuncLimiter) // deduced return type according to values categories C++17
        {
             // NOTE range must be constant to call distance
            if constexpr (std::ranges::range<decltype(aRng)>)
                assert(std::ranges::distance(aRng) == Sfx::EMCNEILNbSections::value);
 
            // gradient over each cell
            std::vector<double> w_deltaU(std::ranges::size(aRng));
 
            // compute the gradient of dU (U_i+1-U_i) over each cell (1st order)
            std::adjacent_difference(std::ranges::begin(aRng), std::ranges::end(aRng), w_deltaU.begin());
 
            // Overwriting first value not part of computation according to 
            // numerical algo "adjacent_difference" and set to 0. for slope limiter
            if (w_deltaU[0] != 0.) // boundary node (right)
            {
                w_deltaU[0] = 0.;   // force to zero
            }
 
            // last node i=100 (ghost node) compare with 0, need to add 
            w_deltaU.push_back(0.); // fwd stencil [i+1-i]
 
            // Apply slope limiter function gradient 1st order (Total Variation Diminishing)
            std::adjacent_difference(w_deltaU.begin(), w_deltaU.end(), // range
                w_deltaU.begin(), // result
                aFuncLimiter);   //function 
 
            return std::move(w_deltaU) | std::views::drop(1);
        }

     // let compiler deduce type (auto always decay)
        template<typename Range, typename FuncLimiter>
        decltype(auto) computeDU( Range& aRng, FuncLimiter&& aFuncLimiter)
        {
#if _DEBUG
        //   static_assert( std::size(aRng) == 101,"Range not correct size");
            assert(Sfx::EMCNEILNbSections::value == std::distance(std::begin(aRng), std::end(aRng)));
#endif
            std::vector<float64> w_deltaU(std::ranges::cbegin(aRng), std::ranges::cend(aRng));
 
            // compute the gradient of U (U_i+1-U_i) over each cell
            std::adjacent_difference(std::ranges::cbegin(aRng), std::ranges::cend(aRng), w_deltaU.begin());
 
            if (w_deltaU[0] != 0.) // boundary node (right)
            {
                w_deltaU[0] = 0.;   // force to zero
            }
 
            // last node i=100 (ghost node) compare with 0, need to add 
            w_deltaU.push_back(0.); //[i+1-i]
 
            std::adjacent_difference(w_deltaU.begin(), w_deltaU.end(), w_deltaU.begin(),
                std::forward<FuncLimiter>(aFuncLimiter));
 
            //auto w_rngDiff = make_iterator_range( w_deltaU.begin(), w_deltaU.end()).advance_begin(1);
 
            // subrange of computational range (without ) 
         // return iterator_range(w_deltaU.begin(), w_deltaU.end()).advance_begin(1);
         // return std::move(w_deltaU) | std::views::drop(w_deltaU,1);
            return std::vector<double>(std::next(w_deltaU.begin()), w_deltaU.end());
        }

        auto computeDU_v( std::ranges::view auto aView)
        {
            // sanity check (concept)
            if constexpr (std::ranges::view<decltype(aView)>)
                static_assert( std::ranges::distance(aView) == Sfx::EMCNEILNbSections::value);
 
            auto view_type = std::ranges::range_value_t<decltype(aView)>;
            std::vector<view_type> w_dU;
            w_dU.reserve(std::ranges::size(aView));
 
            // Compute the gradient of dU (U_i+1-U_i) over each cell (1st order)
            std::adjacent_difference(std::ranges::begin(aView), std::ranges::end(aView), 
                w_dU.begin());
            
            w_dU.push_back(0.); // comparison purpose 
 
            // Apply slope limiter function gradient 1st order (Total Variation Diminishing)
            std::adjacent_difference(w_dU.begin(), w_dU.end(), // range
                w_dU.begin(),   // result
                dbpp::MinMod19{}); // limiter function 
 
            return std::valarray<view_type>( std::next(w_dU.begin()), w_dU.size()-1); // temp fix for debugging
        }
 
        // compute the average over the cell
        template<typename Range>
        auto cell_face_average( const Range& aRng) // consider to pass it b value
        {                                          // std::views such as span is cheap to copy
            // Design Note
            //  call in SweRhsAlgorithm::calculate()  "auto w_Avrg = Sfx::cell_face_average(w_A);"
            //  but w_A not constant.
#if _DEBUG
            // Since we are in the prototyping phase using 
            // E. McNeil dambreak data, particularly nbSections=101
            //static_assert( std::size(aRng) == EMCNEILNbSections::value, "Range not valid size"); error msg "expr did not evaluate to constant"???
            assert( std::ranges::size(aRng) == EMCNEILNbSections::value);
#endif
 
            // concept of a range C++20
            // NOTE
            //  not sure about this one. If we pass valarray I think it would not work.
            //  valrray doesn't provide a begin and end operations
            if constexpr( std::ranges::range<decltype(aRng)>)
            {
                //auto w_rngSiz = std::size(aRng);
              // create array of given size
              std::valarray<float64> w_valArray(std::ranges::size(aRng) - 1);
        // Based on Nujic paper (1995)    
        for (auto i = 0; i < std::ranges::size(aRng) - 1; ++i)
        {
          if (i == 0)
          {
            w_valArray[i] = (0.5 * Sfx::cGravity<float64> *(aRng[i + 1] + aRng[i]));
          }
          else
          {
            w_valArray[i] = (0.5 * Sfx::cGravity<float64> *(aRng[i + 1] + aRng[i - 1]));
          }
        }
              return w_valArray;
            }//if
            else
            {
                throw "Not a range"; // not sure yet, but at least it must do something
                return std::valarray<double>{};
            }
        }
 
    // ===============================================================
    //
    //                      Math Utilities
    //
    // ===============================================================

    inline float64 wrapPi(float64 aTheta)
    {
      aTheta += cPi<float64>;
      aTheta -= static_cast<float>(floor(aTheta * c1Over2Pi<float64>) * c2Pi<float64>);
      aTheta -= cPi<float64>;
 
      return aTheta;
    }

        inline float64 safeAcos(float64 x)
    {
      // Check limit conditions
      if (x <= -1.) {
        return cPi<float64>;
      }
      if (x >= 1.) {
        return 0.;
      }
 
      // Value is in the domain - use standard C function
      return static_cast<float>(acos(x));
    }

        constexpr float64 degToRad( float64 aDeg) { return aDeg * cPiOver180<float64>; }
        constexpr float64 radToDeg( float64 aRad) { return aRad * c180OverPi<float64>; }

        inline void sinCos( float* aReturnSin, float* aReturnCos, float aTheta)
        {
            // For simplicity, we'll just use the normal trig functions.
            // Note that on some platforms we may be able to do better
            *aReturnSin = static_cast<float>(sin(aTheta));
            *aReturnCos = static_cast<float>(cos(aTheta));
        }
 
        // shall return std::optional<double>
        template<typename Range>
        auto StVenant1D_Incomplete_Flux(Range& aU1, Range& aU2)
        {
            static_assert(0 != std::size(aU1));
            static_assert(std::is_same_v<std::size(aU1), std::size(aU2)>);
 
            if constexpr (!std::is_same_v<std::size(aU1), std::size(aU2)>)
            {
                return -1.;
            }
            else
            {
                return (aU2 * aU2) / aU1;
            }
 
            //static_assert(std::extent_v<decltype(aU1[0],0)> == 10.);
            //static_assert(std::extent_v<decltype(aU1[99],0)> == 1.);
            //static_assert( std::size(aU1) == std::size(decltype(aU2)),   "Incomplete Flux range empty");
            //static_assert( !std::size(aU1) && !std::size(aU2)  , "Incomplete Flux range empty");
 
            //return (aU2 * aU2) / aU1;
 
            //  Fonction de calcul de F2 incomplet (excluant le terme de pression hydrostatique)
 
    //      float64 F2;
 
    //      F2 = U2 * U2 / U1;
 
    //      return (F2);
        }
} // End of namespace