c++ - 在数值模拟中使用 std::valarray

Question

我发布了一个用 C++ 编写的简单的 n-body 类 here在代码审查中。

我被告知要使用 std::valarray 而不是普通的 std::array ，目的是我可以重写一些目前看起来像这样的代码

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

对此

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

因为我是 C++ 的初学者，所以我写了一个使用 std::valarray 的小程序，这样我就可以学习如何使用这个数据结构。但是，编译器会抛出很多错误，我不知道为什么。我希望你能帮我解决这个问题。这是我写的小程序:

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

当我尝试编译这段代码时，出现以下错误:

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

首先，当您编写或更改代码时，始终从第一个工作版本开始，并确保代码在每个步骤之间进行编译。这将使隔离错误编译代码变得更加容易。

为什么我要告诉你这个？这是因为您的部分代码从未正确编译过。无论是在引入 valarray 之前还是之后。

例如，这个:

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

单个粒子不是可迭代类型，行尾没有分号。

一步一步来，确保代码在每个步骤之间编译。

---

其次，我认为 valarray 不是您要找的东西。除非您希望每个粒子的每个属性都具有动态的维数，否则这将是非常令人惊讶的。

我建议您引入一个 vec 类型，它包含您进行简化所需的组件。这样的 vec 类型可以在库中找到，例如 glm提供一个 vec2 类，其中包含 +-/* 等运算符。

即使没有库，也可以创建简单的 vector 类型。

这是一个使用 vector (在数学意义上)而不是 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