Add device support to implementation

RUN apt-get update

RUN apt-get install -y --no-install-recommends g++ gcc make cmake zlib1g-dev libjpeg-dev libpng-dev libtinyxml2-dev

RUN apt-get install -y --no-install-recommends git openssh-client subversion procps

RUN apt-get install -y --no-install-recommends vim valgrind gpg man-db

RUN apt-get install -y --no-install-recommends vim valgrind gpg man-db gnuplot

RUN apt-get install -y --no-install-recommends zsh ca-certificates curl netbase wget

RUN git config --global alias.ll 'log --oneline --graph --all --decorate'

#include "main.h"

TNL::Config::ConfigDescription HeatmapSolver::Parameters::makeInputConfig()

{

  TNL::Config::ConfigDescription config;

  config.addEntry<TNL::String>("device", "Device the computation will run on.", "host");

  config.addEntryEnum<TNL::String>("host");

#ifdef HAVE_CUDA

  config.addEntryEnum<TNL::String>("cuda");

#endif

  config.addEntry<int>("grid-x-size", "Grid size along x-axis.", 100);

  config.addEntry<int>("grid-y-size", "Grid size along y-axis.", 100);

  config.addEntry<double>("domain-x-size", "Domain size along x-axis.", 2.0);

  config.addEntry<double>("domain-y-size", "Domain size along y-axis.", 2.0);

  config.addEntry<double>("sigma", "Sigma in exponential initial condition.", 2.0);

  config.addEntry<double>("time-step", "Time step. By default it is proportional to one over space step square.", 0.0);

  config.addEntry<double>("final-time", "Final time of the simulation.", 1.0);

  config.addEntry<bool>("verbose", "Verbose mode.", true);

  return config;

}

HeatmapSolver::Parameters::Parameters(const TNL::Config::ParameterContainer& parameters):

  xSize(parameters.getParameter<int>("grid-x-size")),

  ySize(parameters.getParameter<int>("grid-y-size")),

  xDomainSize(parameters.getParameter<double>("domain-x-size")),

  yDomainSize(parameters.getParameter<double>("domain-y-size")),

  sigma(parameters.getParameter<double>("sigma")),

  timeStep(parameters.getParameter<double>("time-step")),

  finalTime(parameters.getParameter<double>("final-time")),

  verbose(parameters.getParameter<bool>("verbose")) {}

bool HeatmapSolver::Solve(const HeatmapSolver::Parameters& params, TNL::Timer& timer) const {

  // This is always an external storage for grid.

  TNL::Container::Array ux(params.xSize * params.ySize),  // data at step u

                        aux(params.xSize * params.ySize); // data at step u + 1

  // Invalidate ux/aux

  ux = 0;

  aux = 0;

  const double hx = params.xDomainSize / (double)params.xSize;

  const double hy = params.yDomainSize / (double)params.ySize;

  auto uxView = ux.getView(), auxView = aux.getView();

  timer.reset();

  // TODO: - Initial Condition

  auto horizontalBoundaryCondition = [=] __cuda_callable__ (int i) {

    aux[i] = 0;

    aux[(params.ySize - 1) * params.xSize + i] = 0;

  };

  auto verticalBoundaryCondition = [=] __cuda_callable__(int i) {

    aux[j * params.ySize] = 0;

    aux[j * params.ySize + params.xSize - 1] = 0;

  };

  auto next = [=] __cuda_callable__(int i, int j) {

    auto index = j * params.ySize + i;

    aux[index] = (u[c - 1] - 2 * u[c] + u[c + 1]) / hx +

                 (u[c - params.xSize] - 2 * u[c] + u[c + params.xSize]) / hy;

  };

  double time = 0;

  while (time < params.finalTime) {

    // TODO: - Do we really need this

    TNL::Algorithm::ParallelFor(0, params.xSize, horizontalBoundaryCondition);

    TNL::Algorithm::ParallelFor(0, params.ySize, verticalBoundaryCondition);

    TNL::Algorithm::ParallelFor2D(1, 1, params.xSize - 1, params.ySize - 1, next);

    time += params.timeStep;

  }

  return false;

}

#include "main.h"

#include "main.h"

#include<iostream>

#include<fstream>

#include<TNL/Config/parseCommandLine.h>

#include<TNL/Containers/Array.h>

#include <TNL/Timer.h>

#include<TNL/Algorithms/ParallelFor.h>

#pragma once

template<typename Real>

class HeatmapSolver {

  public:

    class Parameters {

      public:

        const int xSize, ySize;

        const double xDomainSize, yDomainSize;

        const double sigma;

        const double timeStep, finalTime;

        const Real xDomainSize, yDomainSize;

        const Real sigma;

        const Real timeStep, finalTime;

        const bool verbose;

        Parameters(const TNL::Config::ParameterContainer &parameters);

        static TNL::Config::ConfigDescription makeInputConfig();

      private:

    };

    bool solve(const Parameters &parameters, TNL::Timer &timer) const;

    template <typename Device>

    bool solve(const Parameters &parameters) const;

  private:

    template <typename Device>

    bool writeGNUPlot(const std::string &filename,

                      const Parameters& parameters,

                      const TNL::Containers::Array<Real, Device>& map) const;

};

template <typename Real>

TNL::Config::ConfigDescription HeatmapSolver<Real>::Parameters::makeInputConfig()

{

  TNL::Config::ConfigDescription config;

  config.addEntry<TNL::String>("device", "Device the computation will run on.", "host");

  config.addEntryEnum<TNL::String>("host");

#ifdef HAVE_CUDA

  config.addEntryEnum<TNL::String>("cuda");

#endif

  config.addEntry<int>("grid-x-size", "Grid size along x-axis.", 100);

  config.addEntry<int>("grid-y-size", "Grid size along y-axis.", 100);

  config.addEntry<Real>("domain-x-size", "Domain size along x-axis.", 2.0);

  config.addEntry<Real>("domain-y-size", "Domain size along y-axis.", 2.0);

  config.addEntry<Real>("sigma", "Sigma in exponential initial condition.", 2.0);

  config.addEntry<Real>("time-step", "Time step. By default it is proportional to one over space step square.", 0.0);

  config.addEntry<Real>("final-time", "Final time of the simulation.", 1.0);

  config.addEntry<bool>("verbose", "Verbose mode.", true);

  return config;

}

template<typename Real>

HeatmapSolver<Real>::Parameters::Parameters(const TNL::Config::ParameterContainer &parameters) : xSize(parameters.getParameter<int>("grid-x-size")),

                                                                                                 ySize(parameters.getParameter<int>("grid-y-size")),

                                                                                                 xDomainSize(parameters.getParameter<Real>("domain-x-size")),

                                                                                                 yDomainSize(parameters.getParameter<Real>("domain-y-size")),

                                                                                                 sigma(parameters.getParameter<Real>("sigma")),

                                                                                                 timeStep(parameters.getParameter<Real>("time-step")),

                                                                                                 finalTime(parameters.getParameter<Real>("final-time")),

                                                                                                 verbose(parameters.getParameter<bool>("verbose")) {}

/***

 *

 * ySize|j                                          (ySize - 1) * xSize + xSize - 1

 *  |------------------------------------------------------

 *  |

 *  |

 *  |

 *  |

 *  |

 *  |

 *  |

 *  |

 *  |

 *  |------------------------------------------------------>

 *

 *  0                                                xSize|i

 *

 *  j * xSize + i

 ***/

template<typename Real>

template<typename Device>

bool HeatmapSolver<Real>::solve(const HeatmapSolver<Real>::Parameters &params) const

{

  // This is always an external storage for grid.

  TNL::Containers::Array<Real, Device> ux(params.xSize * params.ySize),  // data at step u

                                       aux(params.xSize * params.ySize); // data at step u + 1

  // Invalidate ux/aux

  ux = 0;

  aux = 0;

  const Real hx = params.xDomainSize / (Real)params.xSize;

  const Real hy = params.yDomainSize / (Real)params.ySize;

  const Real hx_inv = 1 / (hx * hx);

  const Real hy_inv = 1 / (hy * hy);

  auto timestep = params.timeStep ? params.timeStep : std::min(hx * hx, hy * hy);

  auto uxView = ux.getView(),

       auxView = aux.getView();

  auto init = [=] __cuda_callable__(int i, int j) mutable {

    auto index = j * params.xSize + i;

    uxView[index] = exp(-params.sigma * (i * i + j * j));

  };

  TNL::Algorithms::ParallelFor2D<Device>::exec(1, 1, params.xSize - 1, params.ySize - 1, init);

  if (!writeGNUPlot("data.txt", params, ux))

    return false;

  // auto horizontalBoundaryCondition = [=] __cuda_callable__ (int i) mutable {

  //   auxView[i] = 0;

  //   auxView[(params.ySize - 1) * params.xSize + i] = 0;

  // };

  // auto verticalBoundaryCondition = [=] __cuda_callable__(int j) mutable {

  //   auxView[j * params.xSize] = 0;

  //   auxView[j * params.xSize + params.xSize - 1] = 0;

  // };

  // TNL::Algorithms::ParallelFor<Device>::exec(0, params.xSize, horizontalBoundaryCondition);

  // TNL::Algorithms::ParallelFor<Device>::exec(0, params.ySize, verticalBoundaryCondition);

  auto next = [=] __cuda_callable__(int i, int j) mutable {

    auto index = j * params.ySize + i;

    auxView[index] = (uxView[index - 1] - 2 * uxView[index] + uxView[index + 1]) * hx_inv +

                     (uxView[index - params.xSize] - 2 * uxView[index] + uxView[index + params.xSize]) * hy_inv;

  };

  auto update = [=] __cuda_callable__(int i) mutable {

    uxView[i] += auxView[i] * timestep;

  };

  Real start = 0;

  while (start < params.finalTime)

  {

    TNL::Algorithms::ParallelFor2D<Device>::exec(1, 1, params.xSize - 1, params.ySize - 1, next);

    TNL::Algorithms::ParallelFor<Device>::exec(0, params.xSize * params.ySize, update);

    std::cout << "Iteration: " << start << std::endl;

    start += timestep;

  }

  return writeGNUPlot("data_final.txt", params, ux);

}

template <typename Real>

template <typename Device>

bool HeatmapSolver<Real>::writeGNUPlot(const std::string &filename,

                                       const HeatmapSolver<Real>::Parameters &params,

                                       const TNL::Containers::Array<Real, Device> &map) const {

  std::ofstream out(filename, std::ios::out);

  if (!out.is_open())

    return false;

  for (int j = 0; j < params.ySize; j++)

    for (int i = 0; i < params.xSize; i++)

      out << i << " " << j << " " << map[j * params.xSize + i] << std::endl;

  return out.good();

}

int main(int argc, char * argv[]) {

  TNL::Timer timer;

  auto config = HeatmapSolver::Parameters::makeInputConfig();

  using Real = double;

  auto config = HeatmapSolver<Real>::Parameters::makeInputConfig();

  TNL::Config::ParameterContainer parameters;

  if (!parseCommandLine(argc, argv, config, parameters))

    return EXIT_FAILURE;

  auto params = HeatmapSolver::Parameters(parameters);

  auto device = parameters.getParameter<TNL::String>("device");

  auto params = HeatmapSolver<Real>::Parameters(parameters);

  HeatmapSolver<Real> solver;

  if (device == "host" && !solver.solve<TNL::Devices::Host>(params))

    return EXIT_FAILURE;

#ifdef HAVE_CUDA

  if (device == "cuda" && !solver.solve<TNL::Devices::Cuda>(params))

    return EXIT_FAILURE;

#endif

  return 0;

  return EXIT_SUCCESS;

}