“Modern C++ Scientific Data Type Use case “

Scientific Data Type(scalar field)

Our scientific programming environment supports many scientific data types such as scalar field, vector field, tensor and so on, usual mathematical types used in scientific programming. Recently I started to refactor scalar field using Modern C++20 features (particularly ‘Range’ library). I thought it would be a good example to show how these new C++20 features can be used in our scientific library (programming environment).

Briefly, scalar field type represents a function over a grid in a finite difference grid. It is made up of a grid (holding information about discretization). For example, grid spacing, grid, number of grid points, … and so on. Since in our simulation we use scalar fields in many places, grid must be shared by all other instances, we don’t want the create a grid each time we make a copy or assign to another field. 

Some physics algorithms are templated functions and take range or numerical container as arguments. C++20 ‘Concept’ makes these algorithms more robust by using ‘Concept’. Range ‘Concept’ is a constraint on template argument, for example may require this algorithm only works for type of container (fast-floating point). I am still prototyping with these new features and the more I dig into it, the more I think these can be applied in the development of scientific libraries (more robust and easier to understand).

The best way to learn something is by doing it! and the best way to understand something is to explain it!! 

When I stated the implementation of the scalar field data type, I had in mind Scott Meyer’s advice  

“Easy to use correctly, hard to use incorrectly”, Scott Meyers (More Effective C++) 

 It says it all! That’s it, that’s all!!! End of the story!!! 

Below I present some code snippets showing some implementation with a short description.

STL-like API

Scalar field is STL-like API providing iterator type to traverse field values.  What motivates this choice? I want our data structure to be used with STL algorithm and the new C++20 range library (support the concept of range). C++20 ‘Concept’ is a constraint on the template argument.

using pointer = value_type *;
using iterator = pointer;
using const_pointer = const value_type *;
using const_iterator = const_pointer;
using reference = value_type &;
using const_reference = const value_type &;
using range_diff_type = ptrdiff_t;
using reverse_iterator = std::reverse_iterator<iterator>;
using const_reverse_iterator = std::reverse_iterator<const_iterator>;

/**
 * @brief Iterator access
 * @return iterator
 */
const_iterator cbegin() const;
iterator begin();
const_iterator cend() const;
iterator end();
const_reverse_iterator rbegin() const;

scalar field ctor

We want to be able to initialize a scalar field with initial values. I have implemented a templated constructor that uses C++20 ‘Concept’ (contiguous). I can initialize the scalar field with a std::span (C++20 view which is very cheap to copy and move) but also with standard container such std::vector or with one our numerical types (programming environment) for fast-floating point operations.  All our scientific data structure such as expression template, fast-floating point arrays are STL compatible like interface. I will show below the benefits of this approach. 

// scalar field initialization
template<std::ranges::contiguous_range ContRng> // contiguous_range concept C++20
	scalarField1D( const std::shared_ptr<gridLattice1D>& aGrid, 
	               std::string aName, ContRng aInitValues)
	: m_gridLattice{aGrid},
	  m_fieldName{std::move(aName)}
	  {
			   ....
	  }

Scalar field intialization

Can intitilaize the scalar field using different contiguous storage.

// create shared ptr grid (quick utility to create a one-diemnsional grid)
auto grid1Dptr = std::make_shared<gridLattice1D>(std::string{"d=1 [0,1] [0:99]"});
// one-dimensional scalar field 
scalarField1D w_scal1D{ grid1Dptr, std::string{"scalar field 1D"}};
// sanity check
static_assert( std::is_same_v<decltype(w_scal1D),scalarField1D> == true ); // value type
assert(std::is_lvalue_reference_v<decltype((w_scal1D))>);        // value category
assert(rng::range<decltype(w_scal1D)>);                          // check range
       
// since our field is a range, we can use the range library algorithm
rng::fill( w_scal1D, 2.3); // initialize with values

More initialization (C++20 contiguous range concept)

Since the scalar field is a range and can be initialized from contiguous range, we can initialize from any range that satisfiy this requirement. Below we use vector from standard library and a numerical array taken from our programming environment.

// std vector initialization
std::vector<double> w_contRng(w_scal1D.grid().getNoPoints()); // zero initialized
rng::fill(w_contRng, 1.);
// initial value (contiguous range concept)
scalarField1D w_contRngScal1D{grid1Dptr, std::string{"scalar field 1D"}, w_contRng};
	   
// in-house fast-floating point array (numerical types library)
RealNumArray<double> w_numArray(DIM::value, ptrArr); // 'DIM' metaprogramming constant
if( rng::sized_range<decltype(w_numArray)> && // 'RealNumArray' satisfied C++20 concepts
    rng::contiguous_range<RealNumArray<double>>) // for fast-floating point operations
{
  // initial values for scalar field
  scalarField1D w_numRngScal1D{grid1Dptr, std::string{"scalar field 1D"}, w_numArray};
          ....
}

C++20 Range Utilities

C++20 span view provide an interface to read/write elements stored in contiguous memory.
Contiguous storage (span view is created) to initialize the scalar field (view is cheap to copy/move).

auto w_dynExtent = std::views::counted(rng::begin(w_numArray), rng::size(w_numArray));
if( !rng::empty(w_dynExtent))
{
    scalarField1D w_spnRngScal1D{grid1Dptr, std::string{"scalar field 1D"}, w_dynExtent};         
}

Our numerical types are STL-like API (provides begin/end) and can be manipulated as view or range.

namespace rng = std::ranges; // C++20 range library
// range library C++20 support
rng::fill_n( rng::begin(w_scal1D), rng::size(w_scal1D)/2, 2.5);
rng::fill_n(rng::next( rng::begin(w_scal1D.values()), w_scal1D.values().size()/2),
rng::size(w_scal1D)/2, 2.);

// create a view from adaptor 'all', which convert a range into a ref_view
auto w_vws =  std::views::all(w_mvField); //lvalue

C++20 subrange

Subrange is probably the most handy of the range utilities. It is an iterator/sentinel pair
that models the C++20 View concept. You can use it to bundle together two iterators, or an iterator
and a sentinel, for when you want to return a range or call an API that expects a range.

C++17/20 great feature of subrange is tuple-like, we can unpack the iterator/sentinel pair using
structured bindings (use deduction guide).

// define a subrange from scalar field values (scalar field is a range)
auto [rngFieldBeg, rngFieldEnd] =  // iterator pair
rng::subrange{rng::next(rng::begin(w_contRngScal1D)),  // begin range 
                        rng::prev(rng::end(w_contRngScal1D)),    // end range
                        rng::size(w_contRngScal1D)-2};           // number of elements
			
// can pass it to our physics algorithm (HLL solver part of our generic library) 
hllSchemeFluxAlgo(rngFieldBeg, rngFieldEnd);

Math Operations

Scalar field must provide basic mathematical operations such add, subtract and so on. This is the bread and butter of scientific programming. Also, copy and move semantic constructor must be supported as well. 

// create shared ptr grid
auto grid1Dptr = // unit move VS19 project
std::make_shared<vsc19::gridLattice1D>(std::string{"d=1 [0,1] [0:99]"});
// one-dimensional scalar field from grid
scalarField1D w_scal1D{ grid1Dptr, std::string{"scalar field 1D"}};
		
// copy construct
scalarField1D w_copyField = w_scal1D;
      
// move copy ctor semantic
scalarField1D w_mvField = std::move(w_copyField);

// add 2-scalar field
auto w_sumField = w_scal1D + w_mvField;