test_lifting_gpu.cpp

#include <iostream>
#include <iomanip>
#include <sstream>
#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 "silo_output.hpp"
#include "kernels_gpu.h"

using namespace sgr;

static void
make_geometry(int order, double mesh_h)
{
    gm::add("difftest");

    std::vector<std::pair<int,int>> objects;

    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.5) )
        );
    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.5, 0.5, 0.5, 0.5) )
        );

    std::vector<std::pair<int, int>> tools;
    gmsh::vectorpair odt;
    std::vector<gmsh::vectorpair> odtm;
    gmo::fragment(objects, tools, odt, odtm);

    gmo::synchronize();

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

int test_lifting(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 << reset;
    std::cout << std::endl;

    /* Test field */
    auto Fx = [](const point_3d& pt) -> double {
        return std::sin(M_PI*pt.x()); 
    };

    auto Fy = [](const point_3d& pt) -> double {
        return std::sin(M_PI*pt.y());
    };

    auto Fz = [](const point_3d& pt) -> double {
        return std::sin(M_PI*pt.z());
    };

    auto F = [&](const point_3d& pt) -> vecxd {
        vec3d Fret;
        Fret(0) = Fx(pt);
        Fret(1) = Fy(pt);
        Fret(2) = Fz(pt);
        return Fret;
    };

    /* Expected divergence */
    auto div_F = [](const point_3d& pt) -> double {
        auto dFx_dx = M_PI*std::cos(M_PI*pt.x());
        auto dFy_dy = M_PI*std::cos(M_PI*pt.y());
        auto dFz_dz = M_PI*std::cos(M_PI*pt.z());
        return dFx_dx + dFy_dy + dFz_dz; 
    };

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

#ifdef WRITE_TEST_OUTPUTS
        std::stringstream ss;
        ss << "lift_go_" << geometric_order << "_ao_" << approximation_order;
        ss << "_seq_" << sz << ".silo";

        silo silodb;
        silodb.create_db(ss.str());
        silodb.import_mesh_from_gmsh();
        silodb.write_mesh();
#endif

        auto model_num_dofs = mod.num_dofs();
        auto model_num_fluxes = mod.num_fluxes();

        vecxd Pdiv_F    = vecxd::Zero(model_num_dofs);
        vecxd PFdotn    = vecxd::Zero(model_num_fluxes);
        vecxd LiftF     = vecxd::Zero(model_num_dofs);

        std::vector<entity_data_gpu> edgs;

        for (auto& e : mod)
        {
            e.project(div_F, Pdiv_F);
            entity_data_cpu ed;
            e.populate_entity_data(ed);

            vecxd PFx = e.project(Fx);
            vecxd PFy = e.project(Fy);
            vecxd PFz = e.project(Fz);

            for (size_t iT = 0; iT < e.num_cells(); iT++)
            {
                for (size_t iF = 0; iF < 4; iF++)
                {
                    auto face_det = ed.face_determinants[4*iT+iF];
                    vec3d n = ed.normals.row(4*iT+iF);
                    for (size_t k = 0; k < ed.num_fluxes; k++)
                    {
                        auto base = e.flux_base();
                        auto ofs = (4*iT+iF)*ed.num_fluxes + k;
                        auto Fx = PFx( ed.fluxdofs_mine[ofs] );
                        auto Fy = PFy( ed.fluxdofs_mine[ofs] );
                        auto Fz = PFz( ed.fluxdofs_mine[ofs] );
                        PFdotn(base+ofs) = face_det*Fx*n(0) +
                                           face_det*Fy*n(1) +
                                           face_det*Fz*n(2);
                    }
                }
            }
            entity_data_gpu edg(ed);
            edgs.push_back( std::move(edg) );
        }

        texture_allocator<double> PFdotn_gpu(PFdotn.data(), PFdotn.size());
        device_vector<double> LiftF_gpu(LiftF.data(), LiftF.size());

        for (auto& edg : edgs)
        {
            timecounter_gpu tc;
            tc.tic();
            gpu_compute_flux_lifting(edg, PFdotn_gpu.get_texture(), LiftF_gpu.data());
            double time = tc.toc();

            auto num_cells = edg.num_all_elems;
            if (geometric_order == 1)
            {
                std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
                double flops = 3*(edg.num_bf)*4*edg.num_fluxes*num_cells;
                std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
            }
            else
            {
                //std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
                //double flops = ((21*edg.num_bf+6)*edg.num_bf*edg.num_qp + 3*(2*edg.num_bf-1)*edg.num_bf)*num_cells;
                //std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
            }
        }
        
        LiftF_gpu.copyout(LiftF.data());
        
#ifdef WRITE_TEST_OUTPUTS
        std::vector<double> var_ediv_F( mod.num_cells() );   
        std::vector<double> var_div_F( mod.num_cells() );
        std::vector<double> var_lift( mod.num_cells() );
        std::vector<double> var_vol( mod.num_cells() );
#endif

        /* Compute (∇∙F,g)_T = (F∙n,g)_∂T - (F,∇g)_T on each element T
         * to verify that the lifting computation is correct */

        double err = 0.0;
        for (auto& e : mod)
        {
            for (size_t iT = 0; iT < e.num_cells(); iT++)
            {
                auto& pe = e.cell(iT);
                auto& re = e.cell_refelem(pe);
                auto num_bf = re.num_basis_functions();
                matxd mass = matxd::Zero(num_bf, num_bf);
                vecxd vol = vecxd::Zero(num_bf);

                const auto pqps = pe.integration_points();
                for(size_t iQp = 0; iQp < pqps.size(); iQp++)
                {
                    const auto& pqp = pqps[iQp];
                    const vecxd phi = re.basis_functions(iQp);
                    const matxd dphi = re.basis_gradients(iQp);
                    const mat3d J = pe.jacobian(iQp);
                    mass += pqp.weight() * phi * phi.transpose();
                    vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
                }

                auto ofs = e.cell_global_dof_offset(iT);
                vecxd excpected_div_F = Pdiv_F.segment(ofs, num_bf);
                vecxd Plf = LiftF.segment(ofs, num_bf);
                vecxd Pvol = mass.ldlt().solve(vol);


                vecxd computed_div_F = Plf - Pvol;
                vecxd diff = computed_div_F - excpected_div_F;
                double loc_err = diff.dot(mass*diff);
                err += loc_err;

                auto bar = pe.barycenter();
                vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
#ifdef WRITE_TEST_OUTPUTS
                auto gi = e.cell_global_index_by_gmsh(iT);
                var_ediv_F[gi] = div_F(bar);
                var_div_F[gi] = computed_div_F.dot(phi_bar);
                var_lift[gi] = Plf.dot(phi_bar);
                var_vol[gi] = Pvol.dot(phi_bar);
#endif
            }
        }

        errors.push_back( std::sqrt(err) );

#ifdef WRITE_TEST_OUTPUTS
        silodb.write_zonal_variable("expected_div_F", var_ediv_F);
        silodb.write_zonal_variable("div_F", var_div_F);
        silodb.write_zonal_variable("lift", var_lift);
        silodb.write_zonal_variable("vol", var_vol);
#endif
        std::cout << "Error: " << std::sqrt(err) << std::endl;
    }

    double rate = 0.0;

    std::cout << Byellowfg << "rate" << reset << std::endl;
    for (size_t i = 1; i < sizes.size(); i++)
        std::cout << (rate = std::log2(errors[i-1]/errors[i]) ) << std::endl;

    COMPARE_VALUES_ABSOLUTE("df/dx", rate, 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: LIFTING ***" << reset << std::endl;
    for (size_t ao = 1; ao < 7; ao++)
            test_lifting(1, ao);

    gmsh::finalize();

    return failed_tests;
}