Использование std :: valarray в численном моделировании - PullRequest
Новая
       86

Использование std :: valarray в численном моделировании

0 голосов
/ 28 октября 2019

Я опубликовал простой класс n-body, который я написал на C ++ здесь в Code Review.

Там мне сказали использовать 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) {

Ответы [ 3 ]

1 голос
/ 28 октября 2019

Прежде всего, когда вы пишете или изменяете код, всегда начинайте с первой рабочей версии и убедитесь, что код компилируется между каждым шагом. Это значительно облегчит изоляцию нескомпилируемого кода.

Почему я говорю вам это? Это потому, что есть части вашего кода, которые никогда не компилировались правильно. Неважно, до введения 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, который имеет операторы, такие как +-/* и более.

Даже без библиотеки, простой векторТип может быть создан.

Вот пример вашего кода, использующего векторы (в математическом смысле) вместо 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;
}

Живой пример с пользовательским (неполным)Тип VEC2

1 голос
/ 28 октября 2019

Хорошо, давайте объясним некоторые понятия, чтобы вы могли понять причину ошибок.

std :: valarray против простого массива c ++ против std :: array

Когда мы говорим о простом c ++массив, мы ссылаемся на что-то вроде следующего: 

double myArray[DIM];

- это просто нативная форма массива в c ++. Тогда у нас есть std :: array 

std::array<double, DIM> myStdArry;

, это просто класс-оболочка вокруг простого массива c ++ с той разницей, что у вас есть доступ к некоторым служебным функциям для более простого манипулирования данными, например 

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 (size_t i = 0; i < DIM /*or myStdArry.size()*/; ++i;) {
  myArray[i] = ...;
} 

//Or with range-based for loop

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

Из-за этого вы не можете использовать арифметикуоператоры непосредственно на контейнерах (массиве), так как они не перегружены для них. Это когда valarray входит в картину, так как он был специально разработан для такого рода операций. Valarray - это просто еще один тип массива, который перегружает арифметические операторы, применяя их к каждому элементу массива.

т. Е. Давайте предположим, что мы хотели бы возвести в квадрат каждый элемент массива. С массивом plain и std мы могли бы добиться этого, выполнив: 

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

Но с помощью valarray мы можем просто сделать: 

myValArray *= myValArray;

и мы получим тот же результат.

Другое существенное отличие состоит в том, что хотя обычный массив и массив std имеют фиксированный размер (вы должны установить его размер во время компиляции), массив val может увеличиваться в размерах динамически, и поэтому вам не нужно указывать его размер при объявлении, нодо строительства или позже 

myValArray{DIM} //Or
myValArray.resize(DIM)
1 голос
/ 28 октября 2019

Ваши Particle участники std::array из std::valarray из double , а не просто std::valarray. Это означает, что стандартная библиотека предоставляет только operator+, operator*, [...] для std::valarray, остальное вы должны предоставить самостоятельно.

Например, следующеерешит первую ошибку компилятора из строки position += 0.1*(velocity + force_new); ( Смотрите в прямом эфире ) 

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

Или оберните std::array<std::valarray<double>, DIM> в класси предоставить необходимых операторов для работы. ( Смотреть онлайн в прямом эфире ) 

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