Scientific Modern C++ Features for Scientific Programmer

Scientific Modern C++

Scientific Modern C++ is the application of Modern C++ new programming techniques in developing scientific software applications. Scientific Modern C++ added new features (language upgrade) and many of these newly introduced features can be applied in the development of scientific applications. Scientific application demands efficiency and flexibility, and the new features of Scientific Modern C++ help to achieve these challenges.

I have been programming using Scientific Modern C++ for the last 5 years and I would like to share some thoughts about what I learned and how I see scientific programming in the next years or so. Below I discuss some of the features of Scientific Modern C++ that changed my way of programming or developing scientific applications.

AAA: Almost Always Auto

‘auto’ keyword, ‘decltype(auto)’ is at the heart of modern C++. To program in modern C++ needs a good understanding of how types are deduced. Deep understanding will help you to program in Modern C++. In the book of S. Meyers “Effective Modern C++” item ‘Prefer ‘auto’ to explicit type declaration’, remember of this advice.

First, understand how ‘auto’ deduce type, for that you must have a good grasp how types are deduced template function, this is a must!!! Good reference See book of Scott Meyer (“Effective Modern C++”) and Nicols Josuttis “C++ Template: The Complete Guide” second Edition. Same applies to decltype(auto) introduced in C++14.

There are so many contexts where you must use these keywords (almost everywhere). You have to be aware of all these situations when you need to use one or the other. For example, return of a function with ‘auto’ (let compiler deduce the return type) never return a reference. The reason, ‘auto’ decay, implicit conversion when argument is passed by value, top-level qualifiers such as const and reference are ignored, it means that never return a reference.

Using ‘decltype(auto)’ as a returning type (let compiler deduce type) if you really want returns a reference, a placeholder that deduce return type of an expression according to value category. By using as return type, references are returned by reference and temporaries by value.

rvalue (reference)

‘rvalue’ is one of the new features of Modern C++ and it comes in 2 flavors: standard ‘rvalue reference’ and ‘forward reference’, the terminology is the same, but the behavior is different and their applications as well. Passing arguments as ‘rvalue reference’ allow you to modify these arguments without affecting the correctness of a program. In scientific programming, sometimes we need to modify the passing arguments to compute derivative of a scalar field (adding boundary condition values).

Optimization with ‘rvalue’ or ‘move semantic’, ‘move’ keyword tell the compiler I don’t need this anymore, might as well swap resources instead of copying, which could be more expensive and not necessary. This can have a huge difference when efficiency is important.

Function template with ‘rvalue reference’ for parameters is different from standard ‘rvalue’ reference. Template function rvalue reference are often called forward/universal reference and are called forward/universal reference. It uses ‘automatic template argument type deduction’ to deduces the type. It can bind to any value category and keep valueness and constness of the type. This allows a single function template to accept both lvalue and rvalue parameters (object can be moved into internal storage rather than copied if the parameter is an rvalue.

Main application of these forward/universal reference is in the use of ‘Perfect forwarding’ forward the call to function based on the signature type (forward refence deduction type). It is the overload counterpart of runtime overload (forward/universal reference).

// support C++20 range concept (used with our physics algorithm)
template<typename Range, typename ...Types> 
void hllSchemeFlux( Range&& aRange, Types&&... args) 
{
  namespace rng = std::ranges;
  ...
  // pass argument using forward which in turn is a move since 'aRange' rvalue.
  if constexpr ( std::is_rvalue_reference_v<decltype(aRange)>) // check the type of 'aRange'
  {                                                            // it's an 'rvalue' reference                                                             
     // check the value category of 'aRange' used as expression
     // inside this function it is an lvalue (Range&& aRange) 
     static_assert(std::is_lvalue_reference_v<decltype((aRange))>);
     auto fmt = std::format("Range type is a lvalue: {} ", std::is_lvalue_reference_v<decltype(aRange)>);
     std::cout << std::boolalpha << fmt << '\n';
     ...
     // C++20 Range algorithm 
     // use std::forward to restore the value category (Range&& aRange) 
     // since inside this function it is an 'lvalue'.
     std::ranges::copy( std::forward<Range>(aRange), std::back_inserter(w_out));
      ...
     // contiguous storage (span view is created) w_out is an 'lvalue'
     auto w_spn = std::views::counted(rng::begin(w_out), rng::size(w_out));
     // Compute first derivative (stencil bias upwind backward) 
     auto w_dU = DerivativeBiasUpwind(w_spn); // pass by-value (cheap to copy/move)
     ...
}

Generic programming (variadic template)

Variadic template, is another feature of modern C++ allows you to pass an arbitrary number of template arguments (parameters). You can overload a function with many different signatures (algorithms share many parameters depending on the behavior). Overload can be implemented through what is called “perfect forwarding’, dispatching the call to the right function (signature). This has changed how to program physics algorithms in the mathematical library. From my point view, many of these algorithms are best implemented as function rather than class. It’s more natural and built-in (functional programming).

   /** C++20 template function with universal/forward reference
     *  perfect forwarding preserve lvalue/rvalue-ness (forward to right function)
    */ 
    void fwdSrcAlgo(auto&&... aFwdAlgo)
    {
        // done at compile time!!!
        if constexpr (sizeof...(aFwdAlgo) > 0)
        {
            // calling corresponding function according 
            // to value categorie of each argument
           TreatmentSrcTerms(std::forward<decltype(aFwdAlgo)>(aFwdAlgo)...);
        }
    }

Structured Binding

Structured binding exactly what we need!!Many of our functions (algorithms, solver) return structures such as tensor, pair of values. Structured binding is a perfect tool for functions that return structures or array. Behind the hood structured binding introduced a new name, anonymous entity which is a copy of the initialization type, and elements of this new entity are binding to original object they were initialized from (these are alias name of what they bind to).

It’s important to note that structured binding is extensible, you can add support to structured binding to any type. You can use a tuple-like API to bind names to whatever the API defines as elements. C++20 ranges library type ‘subrange’ is tuple-like, we can bundle together two iterators’ pair using structured bindings.

auto [begIter,endIter] = std::ranges::subrange(....);

Guaranteed Copy Elision

In C++, exchanging data with a function through arguments and return values typically requires copies. In C++17 and above, using return values should almost always be preferred over modifying reference arguments. This holds even for multiple return values, thanks to structured bindings.

//  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

// C++17 structured binding multiples return values
// compute gradient over each cell by applying slope limiter function
auto [dU1,dU2] = computeDU(aU1, aU2, std::forward<F>(aSlopeLimiter));

C++20 Range library (Concept)

Generic programming and scientific programming are tightly bound.Modern C++ provides features such as ‘Trait’, ‘Concept’ to customize the algorithm implementation according to value type, value category and so on. I have already mentioned that our library supports legacy code (C function), since C++ 17 a new trait called ‘invocable’ can be used to query (call) that function doesn’t throw exception (calling pointer-to-function) calling this function no throw exception.

const auto w_ptr2funcName = Sfx::SimulationMgr::getSingleton().getAlgorithmName();
//check equality
if( auto strOrder = w_ptr2funcName <=> std::string{ "HLL_Scheme" }; strOrder == 0) 
{
  std::cout << "Calling the right pointer flux function (HLL_Scheme) validate\n";
}
else {
  std::cerr << "Problem retrieving the pointer function for flux algorithm\n";
}

// Numerical Flux basic algorithm (support legacy code pointer-to-function)
CalculFF w_calculFF = Sfx::SimulationMgr::getSingleton().getAlgorithm(std::string{ "HLL_Scheme" });
assert( nullptr != w_calculFF);

// sanity check (legacy code pointer-2-function are not declare noexcept)
static_assert(!std::is_nothrow_invocable_v<CalculFF, decltype(w_FF1), decltype(w_FF2), decltype(w_U1), decltype(w_U2)>);
    
// std::invoke wrapper functions C++17 !!!! contraint on calling pointer-to-function legacy code
if( std::is_invocable<CalculFF,decltype(w_FF1), decltype(w_FF2), decltype(w_U1), decltype(w_U2)>::value)
{
   // Compute F_i+1/2 face flux (legacy code call pointer-to-function)
   std::invoke(w_calculFF, w_FF1, w_FF2, w_U1, w_U2);
}

Scientific Modern C++ Features for Scientific Programmer