456 * #include <algorithm>
459 * #include <iostream>
462 * #include <opencv2/core.hpp>
463 * #include <opencv2/videoio.hpp>
465 * #include <ginkgo/ginkgo.hpp>
468 *
void set_val(
unsigned char* data,
double value)
470 *
double col_r[] = {255, 221, 129, 201, 249, 255};
471 *
double col_g[] = {255, 220, 130, 161, 158, 204};
472 *
double col_b[] = {255, 220, 133, 93, 24, 8};
473 * value = std::max(0.0, value);
474 *
auto i = std::max(0, std::min(4,
int(value)));
475 *
auto d = std::max(0.0, std::min(1.0, value - i));
476 * data[2] =
static_cast<unsigned char>(col_r[i + 1] * d + col_r[i] * (1 - d));
477 * data[1] =
static_cast<unsigned char>(col_g[i + 1] * d + col_g[i] * (1 - d));
478 * data[0] =
static_cast<unsigned char>(col_b[i + 1] * d + col_b[i] * (1 - d));
482 * std::pair<cv::VideoWriter, cv::Mat> build_output(
int n,
double fps)
484 * cv::Size videosize{n, n};
486 * std::make_pair(cv::VideoWriter{}, cv::Mat{videosize, CV_8UC3});
487 *
auto fourcc = cv::VideoWriter::fourcc(
'a',
'v',
'c',
'1');
488 * output.first.open(
"nls.mp4", fourcc, fps, videosize);
493 *
void output_timestep(std::pair<cv::VideoWriter, cv::Mat>& output,
int n,
494 *
const std::complex<double>* data)
496 *
for (
int i = 0; i < n; i++) {
497 *
auto row = output.second.ptr(i);
498 *
for (
int j = 0; j < n; j++) {
499 * set_val(&row[3 * j],
abs(data[i * n + j]));
502 * output.first.write(output.second);
506 *
int main(
int argc,
char* argv[])
512 *
const auto t0 = 15.0;
513 *
const auto nonlinear_scale = 1.0;
514 *
const auto potential_scale = 3.0;
515 *
const auto time_scale = 0.25;
516 *
const auto n = 256;
517 *
const auto steps_per_sec = 1000;
518 *
const auto fps = 25;
519 *
const auto n2 = n * n;
520 *
const auto h = 2.0 * gko::pi<double>() / n;
521 *
const auto h2 = h * h;
522 *
const auto tau = 1.0 / steps_per_sec;
523 *
const auto idx = [&](
int i,
int j) {
return i * n + j; };
526 * std::ifstream initial_stream(
"data/gko_logo_2d.mtx");
527 * std::ifstream potential_stream(
"data/gko_text_2d.mtx");
528 *
auto amplitude = gko::read<vec>(initial_stream, exec);
529 *
auto potential = gko::read<real_vec>(potential_stream, exec);
530 *
auto frequency =
vec::create(exec, amplitude->get_size());
531 *
auto fft = fft2::create(exec, n, n);
532 *
auto ifft = fft->conj_transpose();
533 *
auto output = build_output(n, fps);
534 *
double last_t = -t0;
535 *
for (
double t = 0; t < t0; t += tau) {
536 *
if (t - last_t > 1.0 / fps) {
538 * std::cout << t << std::endl;
539 * output_timestep(output, n, amplitude->get_const_values());
541 * fft->apply(amplitude, frequency);
542 *
for (
int i = 0; i < n; i++) {
543 *
for (
int j = 0; j < n; j++) {
544 * frequency->at(idx(i, j)) *=
545 * std::polar(1.0, -h2 * (i * i + j * j) * tau * time_scale);
546 * frequency->at(idx(i, j)) *= 1.0 / n2;
549 * ifft->apply(frequency, amplitude);
550 *
for (
int i = 0; i < n; i++) {
551 *
for (
int j = 0; j < n; j++) {
552 * amplitude->at(idx(i, j)) *= std::polar(
553 * 1.0, -(nonlinear_scale *
555 * potential_scale * potential->at(idx(i, j))) *
static std::shared_ptr< OmpExecutor > create(std::shared_ptr< CpuAllocatorBase > alloc=std::make_shared< CpuAllocator >())
Definition executor.hpp:1396
static std::unique_ptr< Dense > create(std::shared_ptr< const Executor > exec, const dim< 2 > &size={}, size_type stride=0)
constexpr std::enable_if_t<!is_complex_s< T >::value, T > abs(const T &x)
Definition math.hpp:931
constexpr auto squared_norm(const T &x) -> decltype(real(conj(x) *x))
Definition math.hpp:913
1.9.8