c++c++11data-structuresstdarrayvalarray

Using std::valarray in numerical simulation


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) {

Solution

  • 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