kokkos_sample.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include <filesystem>
#include <cmath>
#include <chrono>
#include <Kokkos_Core.hpp>
#include <cstdlib>

#include "FDTD_kokkos.h"
#include "FDTD_PML_kokkos.h"

#define __PML_TEST__

using namespace FDTD_kokkos;


void spherical_wave(int n, int it, const std::string base_path = "../../PlotScript/") {
    CurrentParameters cur_param {
        8,
        4,
        0.2
    };

    double T = cur_param.period;
    double Tx = cur_param.period_x;
    double Ty = cur_param.period_y;
    double Tz = cur_param.period_z;

    cur_param.iterations = static_cast<int>(static_cast<double>(cur_param.period) / cur_param.dt);

    std::function<double(double, double, double, double)> cur_func =
        [T, Tx, Ty, Tz](double x, double y, double z, double t) {
        return sin(2.0 * FDTD_const::PI * t / T)
            * pow(cos(2.0 * FDTD_const::PI * x / Tx), 2.0)
            * pow(cos(2.0 * FDTD_const::PI * y / Ty), 2.0)
            * pow(cos(2.0 * FDTD_const::PI * z / Tz), 2.0);
    };

    double d = FDTD_const::C;
    double boundary = static_cast<double>(n) / 2.0 * d;

    Parameters params {
        n,          // Ni
        n,          // Nj
        n,          // Nk
        -boundary,  // x_min
        boundary,   // x_max
        -boundary,  // y_min
        boundary,   // y_max
        -boundary,  // z_min
        boundary,   // z_max
        d,          // dx
        d,          // dy
        d           // dz
    };

    FDTD_kokkos::FDTD method(params, cur_param.dt);

    int cur_time = std::min(cur_param.iterations, it);

    int start_i = static_cast<int>(floor((-Tx / 4.0 - params.ax) / params.dx));
    int start_j = static_cast<int>(floor((-Ty / 4.0 - params.ay) / params.dy));
    int start_k = static_cast<int>(floor((-Tz / 4.0 - params.az) / params.dz));

    int max_i = static_cast<int>(floor((Tx / 4.0 - params.ax) / params.dx));
    int max_j = static_cast<int>(floor((Ty / 4.0 - params.ay) / params.dy));
    int max_k = static_cast<int>(floor((Tz / 4.0 - params.az) / params.dz));

    auto start = std::chrono::high_resolution_clock::now();
    for (int t = 0; t < cur_time; t++) {
        for (int i = start_i; i < max_i; ++i) {
            for (int j = start_j; j < max_j; ++j) {
                for (int k = start_k; k < max_k; ++k) {
                    int index = i + j * params.Ni + k * params.Ni * params.Nj;
                    double value = cur_func(static_cast<double>(i) * params.dx,
                                            static_cast<double>(j) * params.dy,
                                            static_cast<double>(k) * params.dz,
                                            static_cast<double>(t + 1) * cur_param.dt);

                    method.get_field(Component::JX)[index] = value;
                    method.get_field(Component::JY)[index] = value;
                    method.get_field(Component::JZ)[index] = value;
                }
            }
        }
        method.update_fields();
    }
    method.zeroed_currents();
    for (int t = cur_time; t < it; t++) {
        method.update_fields();
    }
    auto end = std::chrono::high_resolution_clock::now();

    std::chrono::duration<double> elapsed = end - start;
    std::cout << "Execution time: " << elapsed.count() << " s" << std::endl;

#ifdef __PML_TEST__
    FDTD_kokkos::FDTD_PML pml_method(params, cur_param.dt, 0.2);

    auto start_pml = std::chrono::high_resolution_clock::now();
    for (int t = 0; t < cur_time; t++) {
        for (int i = start_i; i < max_i; ++i) {
            for (int j = start_j; j < max_j; ++j) {
                for (int k = start_k; k < max_k; ++k) {
                    int index = i + j * params.Ni + k * params.Ni * params.Nj;
                    double value = cur_func(static_cast<double>(i) * params.dx,
                                            static_cast<double>(j) * params.dy,
                                            static_cast<double>(k) * params.dz,
                                            static_cast<double>(t + 1) * cur_param.dt);

                    pml_method.get_field(Component::JX)[index] = value;
                    pml_method.get_field(Component::JY)[index] = value;
                    pml_method.get_field(Component::JZ)[index] = value;
                }
            }
        }
        pml_method.update_fields();
    }
    pml_method.zeroed_currents();
    for (int t = cur_time; t < it; t++) {
        pml_method.update_fields();
    }
    auto end_pml = std::chrono::high_resolution_clock::now();

    std::chrono::duration<double> elapsed_pml = end_pml - start_pml;
    std::cout << "Execution time (PML): " << elapsed_pml.count() << " s" << std::endl;
#endif //__PML_TEST__

    int k = params.Nk/2;
    for (int j = params.Nj/2 - 5; j < params.Nj/2 + 5; j++) {
        for (int i = params.Ni/2 - 5; i < params.Ni/2 + 5; i++) {
            int index = i + j * params.Ni + k * params.Ni * params.Nj;
            std::cout << std::setw(12) << std::fixed << std::setprecision(5) 
                  << method.get_field(Component::EX)[index];
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

#ifdef __PML_TEST__
    std::cout << "PML: \n" << std::endl;
    for (int j = params.Nj/2 - 5; j < params.Nj/2 + 5; j++) {
        for (int i = params.Ni/2 - 5; i < params.Ni/2 + 5; i++) {
            int index = i + j * params.Ni + k * params.Ni * params.Nj;
            std::cout << std::setw(12) << std::fixed << std::setprecision(5) 
                  << pml_method.get_field(Component::EX)[index];
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;
#endif //__PML_TEST__
}

int main(int argc, char* argv[]) {
    Kokkos::initialize(argc, argv); {
        if (Kokkos::hwloc::available())
            std::cout << "hwloc available" << std::endl;
        std::ostringstream msg;
        msg << "hwloc( NUMA[" << Kokkos::hwloc::get_available_numa_count() << "] x CORE["
          << Kokkos::hwloc::get_available_cores_per_numa() << "] x HT["
          << Kokkos::hwloc::get_available_threads_per_core() << "] )" << std::endl;
        Kokkos::print_configuration(msg);
        std::cout << msg.str();


        std::ifstream source_fin;
        std::vector<char*> arguments(argv, argv + argc);

        if (argc == 1) {
            int N = 32;
            int Iterations = 100;
            spherical_wave(N, Iterations, "../../");
        }
        else if (argc == 3) {
            int N = std::atoi(arguments[1]);
            int Iterations = std::atoi(arguments[2]);
            spherical_wave(N, Iterations, "");
        }
        else {
            std::cout << "ERROR: Incorrect number of parameters" << std::endl;
            exit(1);
        }
    }
    Kokkos::finalize();
    return 0;
}