Scientific Programming Using Modern C++ Features

Scientific Programming Modern C++

Scientific programming using Modern C++ requires a new way of thinking about programming, at first sight it may often be more abstract and more difficult to understand than conventional codes. Adopting these new features means not only learning a whole programming style but learning new unfamiliar features as well.

I’m currently working on a migration project of a scientific library to C++20, in this kind of project you must prototype (refactor) and want to use latest features. But you also want to support legacy code, these codes have been tested and used for years and gives good result.

Modern C++ is a lot about generic programming, there are many enhancements about templates, and one of the features I like is the variadic template which allows a generic number of template parameters. Template function can be used, and it is more powerful by combining forward (universal) reference and variadic template functions. Recently I started to use templated functions for prototyping (migration project), support different signature (perfect forwarding), in some context may need different version of the same algorithm and by using template functions provides a simple way of implementing this requirement (support both legacy code and these algorithms are implemented as functions).

Briefly, perfect forwarding allows to forward a call to a given function depending on the value category of the arguments. You can have many implementations of the same functions with different signatures. What makes it possible to do such a thing?  The mechanism that takes place behind the hood is the Automatic Template Argument Type Deduction (property of the forward reference). Template parameter type is deduced based on value categories. 

In my last blog (perfect forwarding), I showed an example of perfect forwarding in selecting different physics algorithms based on signature. 

• Algorithm may depend (use) to implement behavior; by checking type and value categories you can switch between implementations. 
• Algorithm with different signature, some taking lvalue reference, others rvalue reference or both 
• templated function taking forward reference (overload resolution) 

This is a small subset where these concepts take place and need a good understanding of what is going on behind the hood.

Variadic Template

Experimenting an implementation for nodal variable concept used in our application by making use of variadic template.
Nodal values concept, a node can be represented by a set 2,3,4,… ndf (number degree of freedom) the number of attributes is not constant.I think variadic template is best suited to model this.

	template<typename... Ts>
    class NodalTpl 
    {
    public:
        // ctor
        NodalTpl(bool isTied=false) : m_data<Ts...>{}, m_isTied{isTied} {}

        // want to use it with structured binding
        auto getData() const { return m_data; }

        // tuple support this operator
        auto operator<=>( const Nodaltpl &aOther) = default;
        size_t size() const { return std::tuple_size_v<decltype(m_data)>; }
        int nodeNo() { return 1; } // default value
        bool isTiedNode() { return false; } // default value

    protected:
    private:
        // Hold data at node of the global grid
        // For example, can be  variables such state
        // variables etc...
        std::tuple<Ts...> m_data;
        bool m_isTied;
    };

// Example of usage
// create a nodal variable
auto w_nvalue = std::make_from_tuple<NodalTpl>(std::make_tuple(1,2.,3.2,.4));

Templated Functions

Support legacy code using templated function.

// C++20 concept
  template<typename T>
  concept IsPointer = std::is_pointer_v<T>;

  template<IsPointer T, auto N=101>
    void hllSchemeFlux( T aU1, T aU2, T aFF1, T aFF2)
    {
      static_assert( N == vsc19::EMCNEILNbSections::value);
    
      // do some numerical stuff (function-ptr from legacy code)
      // with pointer signature (double*,double*,double*,double*)
      EMcNeilCalculFF(aU1,aU2,aFF1,aFF2); 
    }

Metaprogramming (AAA: Almost Always Auto)

Modern C++ the use of template auto is also a very nice feature that I am using specifically when declaring global functionality of the simulation. Can be used also as return value of function. Because auto decay (implicit conversion when passing argument by-value, it strips const and reference), can’t never return a reference. It may be what you want in some applications. In other cases, you should use decltype(auto), a placeholder for return value of a function based on value categories (can return a reference).

template<auto v>
	struct AppConstant {
		static constexpr auto value = v; // inline variable and a definition since C++17
	};

	/**
	*   @brief Constant for the number of sections (scenario)
	*/
	using EMCNEILNbSections = AppConstant<101>;
	using TESTSIZE5 = AppConstant<5>;
	using DIM = AppConstant<100>;

Support legacy code using variadic template

Prototyping this kind of implementation (migration project). Below some the algorithms signature form legacy code, these algo neeed to be supported by new version of our library

void TreatmentTermeS0( std::vector<double> aVec, double dx, int NbSections) 
  { 
    std::cout << "TreatmentTermeS0 with 3 arguments\n"; 
  }
  
  void TreatmentTermeS0( std::vector<double> aVec, double dx, double grav, int NbSections) 
  { 
    std::cout << "TreatmentTermeS0 with 4 arguments\n"; 
  }

In our programming environment we could use this kind of implenentationto support legacy code (many implementations of this function)
—- NumArray: is a numerical array for fast-floating point operation
—- range_t: is a numerical interface for our numerical container
support C++20 range concept (used by our physics algorithm)

template<typename NumArray, typename... Types> // using universal reference (lvalue/rvalue)
  void numSchemeFlux( NumArray&& aNumArray, Types... args) 
  {
     // C++20 range concept and compile-time if  
    // (initialization and follow by condition C++17)
    if constexpr( range_t<NumArray> numrng; std::ranges::range<decltype(numrng)>)
    { 
      // C++20 utility for ranges
      if (!std::ranges::empty(aNumArray))
      {
        assert(std::ranges::size(aNumArray) == vsc19::DIM::value);
        std::cout << "range with exact value\n";
      }
      // call function from legacy code 
      TreatmentTermeS0(aNumArray, std::forward<Types>(args)...);
    }
    std::cout << "Test varidic algo functions completed\n";
  }

An example taken from our programming environment

// create a grid lattice form a string representation
  auto grid1Dptr = std::make_shared<vsc19::gridLattice1D>(std::string{"d=1 [0,1] [0:99]"});
  // one-dimensional scalar field
  vsc19::scalarField1D w_scal1D{grid1Dptr, std::string{"scalar field 1D"}};
  std::fill_n(std::begin(w_scal1D.values()),w_scal1D.values().size()/2, 2.5);
  // 50 shall be the first index with value 1.0
  std::fill_n(std::next(std::begin(w_scal1D.values()), w_scal1D.values().size()/2), w_scal1D.values().size()/2,1.);

  // calling the right function according to the number of args 
  // cGravity is a variable template to declare universal constant
  vsc19::numSchemeFlux( w_scal1D.asStdVector(), w_scal1D.grid().Delta(), //prvalue: pure reading value
                        vsc19::cGravity<double>, vsc19::DIM::value); // variable template

For some who are interested in scientific programming modern C++, below a link

https://github.com/chavid/ModernScientificCpp/blob/main/1-ClassRoom/README.md