finite_difference_baseline.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <sys/time.h>
#include <vector>
#include <cmath>
#include <chrono>

#include <silo.h>

template<typename T>
class fd_grid
{
    std::vector<T>  m_data;
    size_t m_rows, m_cols;

public:
    fd_grid()
    {}

    fd_grid(size_t rows, size_t cols) :
        m_data(rows*cols), m_rows(rows), m_cols(cols)
    {}

    fd_grid(fd_grid&& other)
    {
        std::swap(m_data, other.m_data);
    }

    fd_grid& operator=(fd_grid&& other)
    {
        std::swap(m_data, other.m_data);
        return *this;
    }

    T operator()(size_t i, size_t j) const
    {
        return m_data[i*m_cols + j];
    }

    T& operator()(size_t i, size_t j)
    {
        return m_data[i*m_cols + j];
    }

    size_t rows() const { return m_rows; }
    size_t cols() const { return m_cols; }

    T* data() { return m_data.data(); }
    const T* data() const { return m_data.data(); }
};

int
visit_dump(const fd_grid<double>& g, const std::string& fn)
{
    DBfile *db = nullptr;
    db = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL, "Kokkos test", DB_HDF5);
    if (!db)
    {
        std::cout << "Cannot write simulation output" << std::endl;
        return -1;
    }

    auto size_x = g.cols();
    auto size_y = g.rows();
    
    std::vector<double> x(size_x);
    std::vector<double> y(size_y);
    
    for (size_t i = 0; i < size_x; i++)
        x.at(i) = double(i)/(size_x-1);
        
    for (size_t i = 0; i < size_y; i++)
        y.at(i) = double(i)/(size_y-1);
        
    int dims[] = { int(size_y), int(size_y) };
    int ndims = 2;
    double *coords[] = {x.data(), y.data()};
    
    DBPutQuadmesh(db, "mesh", NULL, coords, dims, ndims, DB_DOUBLE, DB_COLLINEAR, NULL);
        
    DBPutQuadvar1(db, "solution", "mesh", g.data(), dims, ndims, NULL, 0, DB_DOUBLE, DB_NODECENT, NULL);
    
    DBClose(db);
    return 0;
}

#define HIGH_ORDER

#ifdef HIGH_ORDER
    #define HALO_SIZE 4
#else
    #define HALO_SIZE 1
#endif

int main( int argc, char* argv[] )
{
    using namespace std::chrono;

    int     Nx          = 500;
    int     Ny          = 500;
    int     halo        = HALO_SIZE;

    int     rows        = Ny + 2*halo;
    int     cols        = Nx + 2*halo;
    int     timesteps   = 5000;
    double  dx          = 1./(Nx-1);
    double  dy          = 1./(Ny-1);
    double  dt          = 0.001;

    double  velocity    = 1;
    double  damping     = 0.1;

    auto c2 = velocity*velocity;
    auto dx2 = dx*dx;
    auto dy2 = dy*dy;
    auto dt2 = dt*dt;

    auto D1 = 1./(1. + 0.5*damping*dt);
    auto D2 = 0.5*damping*dt - 1;

    typedef fd_grid<double>  matrix_view_type;

    matrix_view_type u_prev(rows, cols);
    matrix_view_type u_curr(rows, cols);
    matrix_view_type u_next(rows, cols);

    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < cols; j++)
        {
            if ( i < halo or i > rows-halo or j < halo or j > cols-halo)
            {
                u_prev(i,j) = 0.0;
                u_curr(i,j) = 0.0;
            }
            else
            {
                double x = dx*i - 0.5;
                double y = dy*j - 0.1;
                u_prev(i,j) = -std::exp(-2400*(x*x + y*y));
                u_curr(i,j) = 2*dt*u_prev(i,j);
            }
        }
    }

    double iter_time = 0.0;

    double w0 = -205.0/72.0;
    double w1 = 8.0/5.0;
    double w2 = -1.0/5.0;
    double w3 = 8.0/315.0;
    double w4 = -1.0/560.0;
    double w[9] = { w4, w3, w2, w1, w0, w1, w2, w3, w4 };

    for (int ts = 0; ts < timesteps; ts++)
    {
        auto t_begin = high_resolution_clock::now();

        for (size_t i = halo; i < rows - halo; i++)
        {
            for (size_t j = halo; j < cols - halo; j++)
            {
                /* 2nd order in space */
#ifndef HIGH_ORDER
                u_next(i,j) = D1*(2.0*u_curr(i,j) + D2*u_prev(i,j) +
                        + (c2*dt2/dy2)*(u_curr(i+1,j) - 2.0*u_curr(i,j) + u_curr(i-1,j))
                        + (c2*dt2/dx2)*(u_curr(i,j+1) - 2.0*u_curr(i,j) + u_curr(i,j-1)));

                if (j == 1)
                    u_next(i,j-1) = 0;//u_next(i,j);

                if (j == cols-2)
                    u_next(i,j+1) = 0;//u_next(i,j);

                if (i == 1)
                    u_next(i-1,j) = 0;//u_next(i,j);

                if (i == rows-2)
                    u_next(i+1,j) = 0;//u_next(i,j);
#else

                /* 8th order in space */
                double deriv_x = 0.0;
                double deriv_y = 0.0;
                for (int k = -4; k <= 4; k++)
                {
                    deriv_x += w[k+4]*u_curr(i+k,j);
                    deriv_y += w[k+4]*u_curr(i,j+k);
                }
                u_next(i,j) = D1*(2.0*u_curr(i,j) + D2*u_prev(i,j) +
                        + (c2*dt2/dy2)*(deriv_y)
                        + (c2*dt2/dx2)*(deriv_x));

                if (j == halo)
                    for (size_t k = 0; k < halo; k++)
                        u_next(i,k) = 0;

                if (j == cols-(halo+1))
                    for (size_t k = 0; k < halo; k++)
                        u_next(i,j+k+1) = 0;

                if (i == 1)
                    for (size_t k = 0; k < halo; k++)
                        u_next(k,j) = 0;

                if (i == rows-(halo+1))
                    for (size_t k = 0; k < halo; k++)
                        u_next(i+k+1,j) = 0;
#endif
            }
        }

        std::swap(u_prev, u_curr);
        std::swap(u_curr, u_next);

        auto t_end = high_resolution_clock::now();

        duration<double> time_span = duration_cast<duration<double>>(t_end - t_begin);
        iter_time += time_span.count();

        /*
        if ( (ts % 100) == 0 )
        {
            std::stringstream ss;
            ss << "solution_" << ts << ".silo";
            visit_dump(u_curr, ss.str());
        }
        */
    }

    std::cout << "Average iteration time: " << iter_time/timesteps << std::endl;

    return 0;
}