heat-equation.hpp

 
 *    
 *    #include <chrono>
 *    #include <fstream>
 *    #include <iostream>
 *    
 *    #include <opencv2/core.hpp>
 *    #include <opencv2/videoio.hpp>
 *    
 *    #include <ginkgo/ginkgo.hpp>
 *    
 *    
 *    void set_val(unsigned char* data, double value)
 *    {
 *        double col_r[] = {255, 221, 129, 201, 249, 255};
 *        double col_g[] = {255, 220, 130, 161, 158, 204};
 *        double col_b[] = {255, 220, 133, 93, 24, 8};
 *        value = std::max(0.0, value);
 *        auto i = std::max(0, std::min(4, int(value)));
 *        auto d = std::max(0.0, std::min(1.0, value - i));
 *        data[2] = static_cast<unsigned char>(col_r[i + 1] * d + col_r[i] * (1 - d));
 *        data[1] = static_cast<unsigned char>(col_g[i + 1] * d + col_g[i] * (1 - d));
 *        data[0] = static_cast<unsigned char>(col_b[i + 1] * d + col_b[i] * (1 - d));
 *    }
 *    
 *    
 *    std::pair<cv::VideoWriter, cv::Mat> build_output(int n, double fps)
 *    {
 *        cv::Size videosize{n, n};
 *        auto output =
 *            std::make_pair(cv::VideoWriter{}, cv::Mat{videosize, CV_8UC3});
 *        auto fourcc = cv::VideoWriter::fourcc('a', 'v', 'c', '1');
 *        output.first.open("heat.mp4", fourcc, fps, videosize);
 *        return output;
 *    }
 *    
 *    
 *    void output_timestep(std::pair<cv::VideoWriter, cv::Mat>& output, int n,
 *                         const double* data)
 *    {
 *        for (int i = 0; i < n; i++) {
 *            auto row = output.second.ptr(i);
 *            for (int j = 0; j < n; j++) {
 *                set_val(&row[3 * j], data[i * n + j]);
 *            }
 *        }
 *        output.first.write(output.second);
 *    }
 *    
 *    
 *    int main(int argc, char* argv[])
 *    {
 *        using mtx = gko::matrix::Csr<>;
 *        using vec = gko::matrix::Dense<>;
 *    
 *        auto t0 = 5.0;
 *        auto diffusion = 0.0005;
 *        auto source_scale = 2.5;
 *        auto n = 256;
 *        auto steps_per_sec = 500;
 *        auto fps = 25;
 *        auto n2 = n * n;
 *        auto h = 1.0 / (n + 1);
 *        auto h2 = h * h;
 *        auto tau = 1.0 / steps_per_sec;
 *    
 *        auto exec = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
 *        std::ifstream initial_stream("data/gko_logo_2d.mtx");
 *        std::ifstream source_stream("data/gko_text_2d.mtx");
 *        auto source = gko::read<vec>(source_stream, exec);
 *        auto in_vector = gko::read<vec>(initial_stream, exec);
 *        auto out_vector = in_vector->clone();
 *        auto tau_source_scalar = gko::initialize<vec>({source_scale * tau}, exec);
 *        auto stencil_matrix = gko::share(mtx::create(exec));
 *        gko::matrix_data<> mtx_data{gko::dim<2>(n2, n2)};
 *        for (int i = 0; i < n; i++) {
 *            for (int j = 0; j < n; j++) {
 *                auto c = i * n + j;
 *                auto c_val = diffusion * tau * 4.0 / h2 + 1.0;
 *                auto off_val = -diffusion * tau / h2;
 *                if (i > 0) {
 *                    mtx_data.nonzeros.emplace_back(c, c - n, off_val);
 *                }
 *                if (j > 0) {
 *                    mtx_data.nonzeros.emplace_back(c, c - 1, off_val);
 *                }
 *                mtx_data.nonzeros.emplace_back(c, c, c_val);
 *                if (j < n - 1) {
 *                    mtx_data.nonzeros.emplace_back(c, c + 1, off_val);
 *                }
 *                if (i < n - 1) {
 *                    mtx_data.nonzeros.emplace_back(c, c + n, off_val);
 *                }
 *            }
 *        }
 *        stencil_matrix->read(mtx_data);
 *        auto output = build_output(n, fps);
 *        auto solver =
 *            gko::solver::Cg<>::build()
 *                .with_preconditioner(gko::preconditioner::Ic<>::build())
 *                .with_criteria(gko::stop::ResidualNorm<>::build()
 *                                   .with_baseline(gko::stop::mode::rhs_norm)
 *                                   .with_reduction_factor(1e-10))
 *                .on(exec)
 *                ->generate(stencil_matrix);
 *        double last_t = -t0;
 *        for (double t = 0; t < t0; t += tau) {
 *            if (t - last_t > 1.0 / fps) {
 *                last_t = t;
 *                std::cout << t << std::endl;
 *                output_timestep(
 *                    output, n,
 *                    gko::make_temporary_clone(exec->get_master(), in_vector)
 *                        ->get_const_values());
 *            }
 *            in_vector->add_scaled(tau_source_scalar, source);
 *            solver->apply(in_vector, out_vector);
 *            std::swap(in_vector, out_vector);
 *        }
 *    }
 * @endcode
*/