Using std::valarray in numerical simulation

Question

I posted a simple n-body class that I have written in C++ here in Code Review.

There I was told to use std::valarray instead of plain std::array with the goal that I can rewrite some code that looks currently like this

void Particle::update_position() {
    for (unsigned int d = 0; d < DIM; ++d) {
        position[d] += dt*(velocity[d] + a*force_new[d]);
        force_old[d] = force_new[d];
    }
}

to this

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

Since I am a beginner in C++ I wrote a short program that uses the std::valarray such that I can learn how to use this data structure. The compiler, however, throws a lot of errors and I do not know why. I hope you can help me with this. Here is the small program I wrote:

#include <iostream>
#include <vector>
#include <array>
#include <valarray>

constexpr unsigned int DIM = 2;

struct Particle{
    std::array<std::valarray<double>, DIM> position; 
    std::array<std::valarray<double>, DIM> velocity;
    std::array<std::valarray<double>, DIM> force_new;
    std::array<std::valarray<double>, DIM> force_old;
    void update_position();
};

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

void compute_position(std::vector<Particle>& particles) {
    for (auto& particle: particles) { 
        particle.update_position();
    }
}

void print_data(const std::vector<Particle>& particles) {
    for (const auto& particle: particles) {
        for (const auto& x: particle.position) std::cout << x << " ";
        for (const auto& v: particle.position) std::cout << v << " ";
        for (const auto& F: particle.position) std::cout << F << " ";
        std::cout << std::endl;
    }
}

void init_data(std::vector<Particle>& particles) {
    for (auto& particle: particles) {
        for (const auto& p: particle) {
            p.position = 1.0
            p.velocity = 2.0
            p.force_new = 3.0
            p.force_old = 4.0
        }
    }
}

int main() { 
    const unsigned int n = 10;
    std::vector<Particle> particles(n);
    init_data(particles);
    compute_position(particles);
    print_data(particles); 
    return 0;
}

When I try to compile this code, I get the following errors:

so.cpp: In member function ‘void Particle::update_position()’:
so.cpp:17:31: error: no match for ‘operator+’ (operand types are ‘std::array<std::valarray<double>, 2>’ and ‘std::array<std::valarray<double>, 2>’)
     position += 0.1*(velocity + force_new);

so.cpp: In function ‘void print_data(const std::vector<Particle>&)’:
so.cpp:29:58: error: no match for ‘operator<<’ (operand types are ‘std::ostream’ {aka ‘std::basic_ostream<char>’} and ‘const std::valarray<double>’)
         for (const auto& x: particle.position) std::cout << x << " ";


so.cpp: In function ‘void init_data(std::vector<Particle>&)’:
so.cpp:38:29: error: ‘begin’ was not declared in this scope
         for (const auto& p: particle) {
                             ^~~~~~~~
so.cpp:38:29: note: suggested alternative:
In file included from so.cpp:4:
/usr/include/c++/8/valarray:1211:5: note:   ‘std::begin’
     begin(const valarray<_Tp>& __va)
     ^~~~~
so.cpp:38:29: error: ‘end’ was not declared in this scope
         for (const auto& p: particle) {

Answer 1

Your Particle members are std::array of std::valarray of double 's, not just a simple std::valarray. Meaning, the standard library only provides the operator+, operator*, [...] for the std::valarray, the rest you need to provide yourself.

For instance, the following will solves the first compiler error from the line position += 0.1*(velocity + force_new); (See live online)

template<typename ValArray, std::size_t N>
auto operator+(std::array<ValArray, N> lhs, const std::array<ValArray, N>& rhs)
{
   for (std::size_t idx{}; idx < N; ++idx)
      lhs[idx] += rhs[idx];
   return lhs;
}


template<typename ValArray, std::size_t N, typename T>
auto operator*(T constant, std::array<ValArray, N> rhs)
{
   for (std::size_t idx{}; idx < N; ++idx)
      rhs[idx] *= constant;
   return rhs;
}

---

Or wrap the std::array<std::valarray<double>, DIM> to a class and provide necessary operators, to work on. (See live online)

#include <array>
#include <valarray>    

template<typename Type, std::size_t N>
class ValArray2D final
{
   std::valarray<std::valarray<Type>> _arra2D;
public:
   explicit constexpr ValArray2D(const std::valarray<Type>& valArray) noexcept
      : _arra2D{ valArray , N }
   {}

   ValArray2D& operator+=(ValArray2D rhs) noexcept
   {
      _arra2D += rhs._arra2D;
      return *this;
   }

   ValArray2D& operator*(Type constant) noexcept
   {
      for(auto& valArr: this->_arra2D)
         valArr = valArr * constant;
      return *this;
   }

   void operator=(Type constant) noexcept
   {
      for (auto& valArr : this->_arra2D)
         valArr = constant;
   }

   friend ValArray2D operator+(ValArray2D lhs, const ValArray2D& rhs) noexcept
   {
      lhs += rhs;
      return lhs;
   }

   friend std::ostream& operator<<(std::ostream& out, const ValArray2D& obj) noexcept
   {
      for (const std::valarray<Type>& valArray : obj._arra2D)
         for (const Type element : valArray)
            out << element << " ";
      return out;
   }
};

template<typename Type, std::size_t N>
ValArray2D<Type, N> operator*(Type constant, ValArray2D<Type, N> lhs)
{
   return lhs = lhs * constant;
}

Answer 2

First of all, when you write or change code, always start with a first working version, and ensure the code is compiling between each steps. That would make isolating miscompiling code much easier.

Why Am I telling you this? It's because there are parts of your code that has never compiled correctly. No matter before the introduction of valarray or after.

For example, this:

for (auto& particle : particles) {
    for (const auto& p: particle) {
        p.position = 1.0
        p.velocity = 2.0
        p.force_new = 3.0
        p.force_old = 4.0
    }
}

A single particle is not an iterable type, and there is no semicolon at the end of the lines.

Just take it step by step, and ensure the code is compiling between each steps.

---

Second, I don't think valarray is what you're looking for. Unless you want each particle to have a dynamic number of dimension for each property, which would be very surprising.

I'd suggest you to introduce a vec type that would have the component you need to do your simplation. Such vec type can be found in libraries such as glm provide a vec2 class that has operators such as +-/* and more.

Even without a library, a simple vector type can be created.

Here's an example of your code using vectors (in the matematical sense) instead of std::valarray:

struct Particle{
    glm::vec2 position; 
    glm::vec2 velocity;
    glm::vec2 force_new;
    glm::vec2 force_old;
    void update_position();
};

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

void compute_position(std::vector<Particle>& particles) {
    for (auto& particle: particles) { 
        particle.update_position();
    }
}

void print_data(const std::vector<Particle>& particles) {
    for (const auto& particle : particles) {
        std::cout << particle.position.x << ", " << particle.position.y << " ";
        std::cout << particle.velocity.x << ", " << particle.velocity.y << " ";
        std::cout << particle.force_new.x << ", " << particle.force_new.y << " ";
        std::cout << std::endl;
    }
}

void init_data(std::vector<Particle>& particles) {
    for (auto& particle : particles) {
        particle.position = {1, 2};
        particle.velocity = {2, 24};
        particle.force_old = {1, 5};
        particle.force_new = {-4, 2};
    }
}

int main() { 
    const unsigned int n = 10;
    std::vector<Particle> particles(n);
    init_data(particles);
    compute_position(particles);
    print_data(particles); 
    return 0;
}

Live example with custom (imcomplete) vec2 type

Answer 3

Okay, let's explain some concepts so you can understand the reason for the errors.

std::valarray vs plain c++ array vs std::array

When we talk about plain c++ array we refer to something like the follow

double myArray[DIM];

this is just the native form of an array in c++. Then we have std::array

std::array<double, DIM> myStdArry;

this is just a wrapper class around a plain c++ array with the difference that you have access to some utility functions for easier manipulation of the data i.e

myStdArry.size() //get the number of elements
myStdArry.begin() & myStdArry.end() //iterators to the begin and end of the array
myStdArry.fill(1.4) //replace all the values of the array to 1.4

For both plain and std array you will have to use subscript operator ([]) in order to access each element of the array i.e

for (size_t i = 0; i < DIM /*or myStdArry.size()*/; ++i;) {
  myArray[i] = ...;
} 

//Or with range-based for loop

for (auto& element : myArray) {
  element = ...;
}

Due to this, you can't use arithmetic operators directly on the containers (array) since they are not overloaded for them. This is when valarray comes into the picture since it is was specifically designed for this kind of operation. Valarray is just another type of array which overloads arithmetic operators by applying them to each element of the array.

i.e let's suppose that we would like to square each element of an array. With the plain and std array we could achieve it by doing:

for (auto& element : myStdArray) {
  element *= element;
}

But with valarray we can simply do:

myValArray *= myValArray;

and we get the same result.

Other major difference is that while plain and std array are both fixed size (you have to set its size at compile time) val array can grow in size dynamically, and that's why you don't have to specify its size at declaration but until construction or later

myValArray{DIM} //Or
myValArray.resize(DIM)