Added profiling executable.

CMakeLists.txt

@@ -264,6 +264,9 @@ target_link_libraries(test_differentiation gmshdg)
 add_executable(test_differentiation_gpu tests/test_differentiation_gpu.cpp)
 target_link_libraries(test_differentiation_gpu gmshdg)
 
+add_executable(profile_differentiation_gpu tests/profile_differentiation_gpu.cpp)
+target_link_libraries(profile_differentiation_gpu gmshdg)
+
 #add_executable(test test.cpp test_basics.cpp test_mass.cpp test_differentiation.cpp test_differentiation_gpu.cpp test_lifting.cpp ${COMMON_SOURCES})
 #target_link_libraries(test ${LINK_LIBS})
 

tests/profile_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 profile_differentiation(int geometric_order, int approximation_order)
+{
+    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()); 
+    };
+
+        double h = 0.04;
+        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;
+
+    return 0;
+}
+
+int main(int argc, const char *argv[])
+{
+    gmsh::initialize();
+    gmsh::option::setNumber("General.Terminal", 0);
+    gmsh::option::setNumber("Mesh.Algorithm", 1);
+    gmsh::option::setNumber("Mesh.Algorithm3D", 1);
+
+    if (argc != 3)
+    {
+        std::cout << argv[0] << "<gorder> <aorder>" << std::endl;
+        return 1;
+    }
+
+    auto go = atoi(argv[1]);
+    auto ao = atoi(argv[2]);
+
+    profile_differentiation(go, ao);
+
+    return 0;
+}
\ No newline at end of file