Programmation Scientifique C++ Moderne: Data Structure

Programmation Scientifique C++ Moderne

La programmation scientifique est l’une des premières applications de l’ordinateur et demeure une des plus importantes. Aujourd’hui nous sommes en mesure de simuler un très grand éventail de phénomènes physique. C++ est le langage dominant en programmation scientifique. Depuis quelques années le C++ a connu une évolution très rapide et de nouvelles fonctionnalités rendent le langage encore plus performant.

Implementation du Champ Scalaire

Notre environnement de programmation scientifique prend en charge de nombreux types de données scientifiques, tels que les champs scalaires, les champs vectoriels, les tenseurs, etc., types mathématiques couramment utilisés en programmation scientifique.

C++20 contient plusieurs nouvelles fonctionnalités qui trouvent une application en programmation scientifique. Dans ce blog je présente un exemple d’application dans le cas particulier d’un type mathématique appelé champ scalaire.

Parmi toutes ces nouvelles fonctionnalités il y a le concept de ‘Range ». Défini comme étant une séquence de nombre sur lesquels je peux itérer. Un « range » à un début et une fin ou appelé sentinelle dans la nouvelle mouture C++20. Notre librairie mathématique supporte de nombreux types (supporte STL-like API). Pourquoi? parce que je peux utiliser ceux-ci avec les algorithmes de la STL. Mais depuis C++20 ceux-ci satisfont au concept de Range. Je donne un exemple ci-dessous. 

STL API

Le champ scalaire est une API de type STL fournissant un type d’itérateur pour parcourir les valeurs du champ. Qu’est-ce qui motive ce choix ? Je souhaite que notre structure de données soit utilisée avec les algorithmes de la STL et la nouvelle bibliothèque de « Range »C++20.

using pointer = value_type *;
using iterator = pointer;
using const_pointer = const value_type *;
using const_iterator = const_pointer;
uing 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();

scalar field ctor

Nous souhaitons pouvoir initialiser un champ scalaire avec des valeurs initiales. J’ai implémenté un constructeur basé sur un modèle utilisant le langage C++20 « Concept » (contigu). Je peux initialiser le champ scalaire avec un std::span (C++20), mais aussi avec un conteneur standard tel que std::vector ou avec l’un de nos types numériques. Toutes nos structures de données scientifiques, telles que les modèles d’expression et les tableaux à virgule flottante rapide, sont compatibles avec l’interface STL. Je détaillerai ci-dessous les avantages de cette approche.

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

Ci-dessous un exemple d’initialization d’un champ scalaire avec differents containeurs contigus.

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

D’autres types d’initialization

Le champ scalaire étant une plage et pouvant être initialisé à partir d’une plage contiguë, nous pouvons initialiser à partir de n’importe quelle plage satisfaisant à cette exigence. Ci-dessous, nous utilisons un vecteur de la bibliothèque standard et un tableau numérique issu de notre environnement de programmation.

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

Subrange est probablement l’utilitaire de plage le plus pratique. Il s’agit d’un couple itérateur/sentinelle qui modélise le concept de view C++20. Vous pouvez l’utiliser pour regrouper deux itérateurs, ou un itérateur et une sentinelle.

L’une des fonctionnalités intéressantes de subrange est tuple-like API : nous pouvons décomposer le couple itérateur/sentinelle à l’aide de structured binding C++17.

namespace rng = std::ranges;

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

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