test_differentiation_gpu.cpp

#include <iostream>
#include <iomanip>

#include <unistd.h>

#include "test.h"
#include "gmsh_io.h"
#include "sgr.hpp"
#include "entity_data.h"
#include "kernels_cpu.h"
#include "timecounter.h"

#include "kernels_gpu.h"

using namespace sgr;

static void
make_geometry(int order, double mesh_h)
{
    gm::add("difftest");
    gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
    gmo::synchronize();

    gvp_t vp;
    gm::getEntities(vp);
    gmm::setSize(vp, mesh_h);
}

int test_differentiation_convergence(int geometric_order, int approximation_order)
{
    std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
    std::vector<double> errors;

    std::cout << cyanfg << "Testing geometric order " << geometric_order;
    std::cout << ", approximation order = " << approximation_order << nofg;
    std::cout << std::endl;

    auto f = [](const point_3d& pt) -> double {
        return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z()); 
    };
    auto df_dx = [](const point_3d& pt) -> double {
        return M_PI*std::cos(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z()); 
    };
    auto df_dy = [](const point_3d& pt) -> double {
        return M_PI*std::sin(M_PI*pt.x())*std::cos(M_PI*pt.y())*std::sin(M_PI*pt.z()); 
    };
    auto df_dz = [](const point_3d& pt) -> double {
        return M_PI*std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::cos(M_PI*pt.z()); 
    };

    std::vector<double> errs_x;
    std::vector<double> errs_y;
    std::vector<double> errs_z;

    for (size_t i = 0; i < sizes.size(); i++)
    {
        double h = sizes[i];
        make_geometry(0,h);
        model mod(geometric_order, approximation_order);

        auto& e0 = mod.entity_at(0);
        e0.sort_by_orientation();
        vecxd Pf_cpu = e0.project(f);

        /* Make CPU entity data */
        entity_data_cpu ed;
        e0.populate_entity_data(ed);
        /* Derive GPU entity data */
        entity_data_gpu edg(ed);

        /* Prepare I/O vectors and call kernel */
        device_vector<double> Pf_gpu(Pf_cpu.data(), Pf_cpu.size());
        device_vector<double> df_dx_gpu(Pf_cpu.size());
        device_vector<double> df_dy_gpu(Pf_cpu.size());
        device_vector<double> df_dz_gpu(Pf_cpu.size());

        timecounter_gpu tc;
        tc.tic();
        gpu_compute_field_derivatives(edg, Pf_gpu.data(), df_dx_gpu.data(),
            df_dy_gpu.data(), df_dz_gpu.data());
        double time = tc.toc();

        vecxd df_dx_cpu = vecxd::Zero( Pf_cpu.size() );
        vecxd df_dy_cpu = vecxd::Zero( Pf_cpu.size() );
        vecxd df_dz_cpu = vecxd::Zero( Pf_cpu.size() );

        df_dx_gpu.copyout( df_dx_cpu.data() );
        df_dy_gpu.copyout( df_dy_cpu.data() );
        df_dz_gpu.copyout( df_dz_cpu.data() );

        if (geometric_order == 1)
        {
            std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
            double flops = 21*ed.num_bf*ed.num_bf*e0.num_cells();
            std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
        }
        else
        {
            std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
            double flops = ((21*ed.num_bf+6)*ed.num_bf*ed.num_qp + 3*(2*ed.num_bf-1)*ed.num_bf)*e0.num_cells();
            std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
        }

        vecxd Pdf_dx_e0 = e0.project(df_dx);
        vecxd Pdf_dy_e0 = e0.project(df_dy);
        vecxd Pdf_dz_e0 = e0.project(df_dz);

        double err_x = 0.0;
        double err_y = 0.0;
        double err_z = 0.0;

        for (size_t iT = 0; iT < e0.num_cells(); iT++)
        {
            auto& pe = e0.cell(iT);
            auto& re = e0.cell_refelem(pe);
            matxd mass = e0.mass_matrix(iT);
            auto num_bf = re.num_basis_functions();
            vecxd diff_x = Pdf_dx_e0.segment(iT*num_bf, num_bf) - df_dx_cpu.segment(iT*num_bf, num_bf); 
            vecxd diff_y = Pdf_dy_e0.segment(iT*num_bf, num_bf) - df_dy_cpu.segment(iT*num_bf, num_bf); 
            vecxd diff_z = Pdf_dz_e0.segment(iT*num_bf, num_bf) - df_dz_cpu.segment(iT*num_bf, num_bf); 
        
            err_x += diff_x.dot(mass*diff_x);
            err_y += diff_y.dot(mass*diff_y);
            err_z += diff_z.dot(mass*diff_z);
        }
        std::cout << "Errors: " << std::sqrt(err_x) << " " << std::sqrt(err_y);
        std::cout << " " << std::sqrt(err_z) << std::endl;
    
        errs_x.push_back( std::sqrt(err_x) );
        errs_y.push_back( std::sqrt(err_y) );
        errs_z.push_back( std::sqrt(err_z) );
    }

    double rate_x = 0.0;
    double rate_y = 0.0;
    double rate_z = 0.0;

    std::cout << Byellowfg << "rate df/dx  rate df/dy  rate df/dz" << reset << std::endl;
    for (size_t i = 1; i < sizes.size(); i++)
    {
        std::cout << (rate_x = std::log2(errs_x[i-1]/errs_x[i]) ) << " ";
        std::cout << (rate_y = std::log2(errs_y[i-1]/errs_y[i]) ) << " ";
        std::cout << (rate_z = std::log2(errs_z[i-1]/errs_z[i]) ) << std::endl;
    }

    COMPARE_VALUES_ABSOLUTE("df/dx", rate_x, double(approximation_order), 0.2);
    COMPARE_VALUES_ABSOLUTE("df/dy", rate_y, double(approximation_order), 0.2);
    COMPARE_VALUES_ABSOLUTE("df/dz", rate_z, double(approximation_order), 0.2);

    return 0;
}

int main(void)
{
    gmsh::initialize();
    gmsh::option::setNumber("General.Terminal", 0);
    gmsh::option::setNumber("Mesh.Algorithm", 1);
    gmsh::option::setNumber("Mesh.Algorithm3D", 1);


    int failed_tests = 0;

    std::cout << Bmagentafg << " *** TESTING: DIFFERENTIATION ***" << reset << std::endl;
    for (size_t go = 1; go < 2; go++)
        for (size_t ao = go; ao < 7; ao++)
            failed_tests += test_differentiation_convergence(go, ao);

    return failed_tests;
}