forked from NSCCN/hpc-app
52 lines
1.3 KiB
C++
52 lines
1.3 KiB
C++
#include "stencil.hpp"
|
|
|
|
#include <cmath>
|
|
#include <numeric>
|
|
#include <stdexcept>
|
|
|
|
namespace hpc_demo {
|
|
|
|
std::vector<double> make_initial_field(std::size_t size, std::size_t global_offset) {
|
|
std::vector<double> values(size);
|
|
constexpr double pi = 3.14159265358979323846;
|
|
|
|
for (std::size_t i = 0; i < size; ++i) {
|
|
const double x = static_cast<double>(global_offset + i);
|
|
values[i] = std::sin(x * pi / 180.0) + 0.1 * std::cos(x * pi / 17.0);
|
|
}
|
|
|
|
return values;
|
|
}
|
|
|
|
void stencil_step(const std::vector<double>& current,
|
|
std::vector<double>& next,
|
|
double left_ghost,
|
|
double right_ghost,
|
|
double alpha) {
|
|
if (current.empty()) {
|
|
next.clear();
|
|
return;
|
|
}
|
|
|
|
if (next.size() != current.size()) {
|
|
next.assign(current.size(), 0.0);
|
|
}
|
|
|
|
if (alpha < 0.0 || alpha > 0.5) {
|
|
throw std::invalid_argument("alpha must be in [0, 0.5] for a stable diffusion step");
|
|
}
|
|
|
|
for (std::size_t i = 0; i < current.size(); ++i) {
|
|
const double left = (i == 0) ? left_ghost : current[i - 1];
|
|
const double right = (i + 1 == current.size()) ? right_ghost : current[i + 1];
|
|
next[i] = current[i] + alpha * (left - 2.0 * current[i] + right);
|
|
}
|
|
}
|
|
|
|
double checksum(const std::vector<double>& values) {
|
|
return std::accumulate(values.begin(), values.end(), 0.0);
|
|
}
|
|
|
|
} // namespace hpc_demo
|
|
|