Implement Field memory region backend infrastructure

CMakeLists.txt

@@ -3,12 +3,14 @@ project(Redesign CXX)
 
 option(NEKTAR_USE_CUDA "Enable CUDA support" FALSE)
 
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
 set(SRC LibUtilities/ErrorUtil.cpp Operators/Operator.cpp Operators/OperatorBwdTrans.cpp Operators/BwdTrans/BwdTransImpl.cpp)
 
 if (NEKTAR_USE_CUDA)
     enable_language(CUDA)
     add_definitions(-DNEKTAR_USE_CUDA)
-    set(SRC ${SRC} MemoryRegionCUDA.cu)
+    set(SRC ${SRC} MemoryRegionCUDA.cu Operators/BwdTrans/BwdTransCUDA.cu)
 endif()
 
 add_library(Operators SHARED ${SRC})
@@ -17,4 +19,4 @@ target_include_directories(Operators PRIVATE "${CMAKE_SOURCE_DIR}")
 set(TEST_SRC main.cpp)
 add_executable(main ${TEST_SRC})
 target_link_libraries(main Operators)
-target_include_directories(main PRIVATE "${CMAKE_SOURCE_DIR}")
\ No newline at end of file
+target_include_directories(main PRIVATE "${CMAKE_SOURCE_DIR}")

Field.hpp

 #pragma once
 
 #include <array>
+#include <memory>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
-#include "MemoryRegion.hpp"
 #include "MemoryRegionCPU.hpp"
 
-#ifdef NEKTAR_USE_CUDA
-#include "MemoryRegionCUDA.hpp"
-#endif
-
 enum class ShapeType
 {
     eQuadrilateral,
@@ -32,14 +29,34 @@ enum class FieldState
 
 static constexpr FieldState DefaultState = FieldState::Phys;
 
-template <typename TType = double, FieldState TState = DefaultState,
-          typename TBackend = DefaultBackend>
-class Field
+template <typename TType = double, FieldState TState = DefaultState> class Field
 {
 public:
-    Field(std::vector<BlockAttributes> &blocks, int num_components = 1)
-        : block_attributes(blocks)
+    Field(const Field &) = delete;
+    virtual ~Field()     = default;
+
+    Field(Field &&rhs)
+        : m_storage(std::move(rhs.m_storage)),
+          block_attributes(std::move(rhs.block_attributes)),
+          component_names(std::move(component_names))
     {
+    }
+
+    Field &operator=(Field &&rhs)
+    {
+        m_storage        = std::move(rhs.m_storage);
+        block_attributes = std::move(rhs.block_attributes);
+        component_names  = std::move(rhs.component_names);
+
+        return *this;
+    }
+
+    template <template <typename> class TMemoryRegion = MemoryRegionCPU>
+    static Field<TType, TState> create(std::vector<BlockAttributes> &blocks,
+                                       int num_components = 1)
+    {
+        auto field = Field(blocks, num_components);
+
         size_t storage_size = 0;
         for (int i = 0; i < blocks.size(); ++i)
         {
@@ -50,15 +67,30 @@ public:
             storage_size += blockSize;
         }
 
-        m_storage =
-            MemoryRegion<TType, TBackend>(storage_size * num_components);
+        field.m_storage = std::make_unique<TMemoryRegion<TType>>(
+            storage_size * num_components);
+
+        return field;
     }
-    Field(const Field &) = delete;
-    ~Field()             = default;
 
-    MemoryRegion<TType, TBackend> &GetStorage()
+    template <template <typename> class TMemoryRegion = MemoryRegionCPU>
+    TMemoryRegion<TType> &GetStorage()
     {
-        return m_storage;
+        static_assert(std::is_base_of<MemoryRegionCPU<TType>,
+                                      TMemoryRegion<TType>>::value,
+                      "TMemoryRegion must derive MemoryRegionCPU<TType>");
+
+        try
+        {
+            return dynamic_cast<TMemoryRegion<TType> &>(*m_storage);
+        }
+        catch (const std::bad_cast &e)
+        {
+            throw std::runtime_error(
+                "Failed to cast memory storage from type " +
+                std::string(typeid(*m_storage).name()) + " to type " +
+                std::string(typeid(TMemoryRegion<TType>).name()));
+        }
     }
 
     size_t GetNumComponents()
@@ -67,7 +99,12 @@ public:
     }
 
 private:
-    MemoryRegion<TType, TBackend> m_storage;
+    Field(std::vector<BlockAttributes> &blocks, int num_components = 1)
+        : block_attributes(blocks)
+    {
+    }
+
+    std::unique_ptr<MemoryRegionCPU<TType>> m_storage;
     // std::vector<MemoryView<TType>> m_views;
     std::vector<BlockAttributes> block_attributes;
     std::vector<std::string> component_names;

MemoryRegion.hpp

-#pragma once
-
-struct BackendCPU;
-struct BackendCUDA;
-
-#ifdef NEKTAR_USE_CUDA
-using DefaultBackend = BackendCUDA;
-#else
-using DefaultBackend = BackendCPU;
-#endif
-
-template <typename tData, typename tBackend = DefaultBackend>
-class MemoryRegion;

MemoryRegionCPU.hpp

 #pragma once
-#include "MemoryRegion.hpp"
 
-template <typename tData> class MemoryRegion<tData, BackendCPU>
+#include <cstdio>
+
+template <typename TData> class MemoryRegionCPU
 {
 public:
-    MemoryRegion()                        = default;
-    MemoryRegion(const MemoryRegion &rhs) = delete;
+    MemoryRegionCPU()                           = delete;
+    MemoryRegionCPU(const MemoryRegionCPU &rhs) = delete;
 
-    MemoryRegion(MemoryRegion &&rhs)
+    MemoryRegionCPU(MemoryRegionCPU &&rhs)
     {
         m_host     = rhs.m_host;
         m_size     = rhs.m_size;
         rhs.m_host = nullptr;
     }
 
-    MemoryRegion(size_t n)
+    MemoryRegionCPU(size_t n)
     {
-        m_host = new tData[n];
+        m_host = new TData[n];
         m_size = n;
     }
 
-    virtual ~MemoryRegion()
+    virtual ~MemoryRegionCPU()
     {
         if (m_host != nullptr)
         {
@@ -29,23 +30,24 @@ public:
         }
     }
 
-    void operator=(MemoryRegion &&rhs)
+    void operator=(MemoryRegionCPU &&rhs)
     {
         m_host     = rhs.m_host;
         m_size     = rhs.m_size;
         rhs.m_host = nullptr;
     }
 
-    virtual double *GetPtr()
+    virtual TData *GetCPUPtr()
     {
         return m_host;
     }
 
-    size_t size() {
+    size_t size()
+    {
         return m_size;
     }
 
 protected:
-    double *m_host = nullptr;
-    size_t m_size  = 0;
+    TData *m_host = nullptr;
+    size_t m_size = 0;
 };

MemoryRegionCUDA.cu

@@ -3,40 +3,40 @@
 
 #include "MemoryRegionCUDA.hpp"
 
-template<typename tData>
-MemoryRegion<tData, BackendCUDA>::MemoryRegion(size_t n)
-    : MemoryRegion<tData, BackendCPU>(n)
+template <typename TData>
+MemoryRegionCUDA<TData>::MemoryRegionCUDA(size_t n) : MemoryRegionCPU<TData>(n)
 {
-    //m_host = new tData[n];
-    cudaMalloc((void **)&m_device, sizeof(tData) * n);
+    // m_host = new TData[n];
+    cudaMalloc((void **)&m_device, sizeof(TData) * n);
     m_size = n;
-    std::cout << "HEllo" << std::endl;
 }
 
-template<typename tData>
-MemoryRegion<tData, BackendCUDA>::~MemoryRegion() {
-    if (m_device != nullptr) {
+template <typename TData> MemoryRegionCUDA<TData>::~MemoryRegionCUDA()
+{
+    if (m_device != nullptr)
+    {
         cudaFree(m_device);
         m_device = nullptr; // lol
     }
 }
 
-template<typename tData>
-void MemoryRegion<tData, BackendCUDA>::HostToDevice()
+template <typename TData> void MemoryRegionCUDA<TData>::HostToDevice()
 {
-    cudaMemcpy(m_device, this->m_host, m_size*sizeof(tData), cudaMemcpyHostToDevice);
+    cudaMemcpy(m_device, this->m_host, m_size * sizeof(TData),
+               cudaMemcpyHostToDevice);
     m_ondevice = true;
 }
 
-template<typename tData>
-void MemoryRegion<tData, BackendCUDA>::DeviceToHost()
+template <typename TData> void MemoryRegionCUDA<TData>::DeviceToHost()
 {
-    cudaMemcpy(this->m_host, m_device, m_size*sizeof(tData), cudaMemcpyDeviceToHost);
+    cudaMemcpy(this->m_host, m_device, m_size * sizeof(TData),
+               cudaMemcpyDeviceToHost);
     m_ondevice = false;
 }
 
-template MemoryRegion<double, BackendCUDA>::MemoryRegion(size_t n);
-template MemoryRegion<double, BackendCUDA>::~MemoryRegion();
-template void MemoryRegion<double, BackendCUDA>::HostToDevice();
-template void MemoryRegion<double, BackendCUDA>::DeviceToHost();
-
+template class MemoryRegionCUDA<double>;
+//
+// template MemoryRegionCUDA<double>::MemoryRegionCUDA(size_t n);
+// template MemoryRegionCUDA<double>::~MemoryRegionCUDA();
+// template void MemoryRegionCUDA<double>::HostToDevice();
+// template void MemoryRegionCUDA<double>::DeviceToHost();

MemoryRegionCUDA.hpp

 #pragma once
-#include "MemoryRegion.hpp"
+
+#include <utility>
+
 #include "MemoryRegionCPU.hpp"
 
-template<typename tData>
-class MemoryRegion<tData, BackendCUDA> : 
-        public MemoryRegion<tData, BackendCPU>
+template <typename TData> class MemoryRegionCUDA : public MemoryRegionCPU<TData>
 {
 public:
-    MemoryRegion() : MemoryRegion<tData, BackendCPU>() {}
-    MemoryRegion(const MemoryRegion &rhs) = delete;
-    MemoryRegion(MemoryRegion &&rhs) : MemoryRegion<tData, BackendCPU>(rhs) {
+    MemoryRegionCUDA(const MemoryRegionCUDA<TData> &rhs) = delete;
+    MemoryRegionCUDA(MemoryRegionCUDA &&rhs)
+        : MemoryRegionCPU<TData>(std::move(rhs))
+    {
         // C++ is less retarded
-        m_device = rhs.m_device;
-        m_size = rhs.m_size;
+        m_device     = rhs.m_device;
+        m_size       = rhs.m_size;
         rhs.m_device = nullptr;
-        rhs.m_size = 0;
+        rhs.m_size   = 0;
     }
-    MemoryRegion(size_t n);
-    virtual ~MemoryRegion() override;
+    MemoryRegionCUDA(size_t n);
+    virtual ~MemoryRegionCUDA() override;
 
-    void operator=(MemoryRegion &&rhs) {
-        MemoryRegion<tData, BackendCPU>::operator=(std::move(rhs));
-        m_device = rhs.m_device;
-        m_size = rhs.m_size;
+    void operator=(MemoryRegionCUDA &&rhs)
+    {
+        MemoryRegionCPU<TData>::operator=(std::move(rhs));
+        m_device     = rhs.m_device;
+        m_size       = rhs.m_size;
         rhs.m_device = nullptr;
-        rhs.m_size = 0;
+        rhs.m_size   = 0;
     }
 
-    virtual double *GetPtr() override
+    virtual TData *GetCPUPtr() override
     {
-        return m_ondevice ? m_device : this->m_host;
+        if (m_ondevice)
+        {
+            DeviceToHost();
+        }
+
+        return this->m_host;
+    }
+
+    TData *GetGPUPtr()
+    {
+        if (!m_ondevice)
+        {
+            HostToDevice();
+        }
+
+        return m_device;
     }
 
     void HostToDevice();
@@ -41,10 +58,7 @@ public:
     }
 
 protected:
-    // double *m_host = nullptr;
-    double *m_device = nullptr;
-    size_t m_size = 0;
+    TData *m_device = nullptr;
+    size_t m_size   = 0;
     bool m_ondevice = false;
-
 };
-

Operators/BwdTrans/BwdTransCUDA.cu

+#include "BwdTransCUDA.hpp"
+
+namespace Nektar::Operators::detail
+{
+template <>
+std::string OperatorBwdTransImpl<double, ImplCUDA>::className =
+    GetOperatorFactory<double>().RegisterCreatorFunction(
+        "BwdTransCUDA", OperatorBwdTransImpl<double, ImplCUDA>::instantiate,
+        "...");
+}

Operators/BwdTrans/BwdTransCUDA.hpp

+#include "MemoryRegionCUDA.hpp"
 #include "Operators/OperatorBwdTrans.hpp"
 
 namespace Nektar::Operators::detail
@@ -11,15 +12,17 @@ public:
     void apply(Field<TData, FieldState::Coeff> &in,
                Field<TData, FieldState::Phys> &out) override
     {
+        TData *x = in.template GetStorage<MemoryRegionCUDA>().GetGPUPtr();
+        TData *y = out.template GetStorage<MemoryRegionCUDA>().GetGPUPtr();
         std::cout << "Op bwd trans CUDA\n";
     }
 
     static std::unique_ptr<Operator<TData>> instantiate()
     {
-        return std::make_unique<OperatorBwdTransImpl<TData, ImplMatFree>>();
+        return std::make_unique<OperatorBwdTransImpl<TData, ImplCUDA>>();
     }
 
     static std::string className;
 };
 
-}
\ No newline at end of file
+} // namespace Nektar::Operators::detail

Operators/BwdTrans/BwdTransMatFree.hpp

@@ -11,10 +11,11 @@ public:
     void apply(Field<TData, FieldState::Coeff> &in,
                Field<TData, FieldState::Phys> &out) override
     {
-        TData* x = in.GetStorage().GetPtr();
-        TData* y = out.GetStorage().GetPtr();
+        TData *x       = in.GetStorage().GetCPUPtr();
+        TData *y       = out.GetStorage().GetCPUPtr();
         const size_t n = in.GetStorage().size();
-        for (size_t i = 0; i < n; ++i) {
+        for (size_t i = 0; i < n; ++i)
+        {
             y[i] = 2.0 * x[i];
         }
         std::cout << "Op bwd trans mat free\n";
@@ -28,4 +29,4 @@ public:
     static std::string className;
 };
 
-}
\ No newline at end of file
+} // namespace Nektar::Operators::detail

main.cpp

@@ -4,6 +4,10 @@
 #include "Field.hpp"
 #include "Operators/OperatorBwdTrans.hpp"
 
+#ifdef NEKTAR_USE_CUDA
+#include "MemoryRegionCUDA.hpp"
+#endif
+
 using namespace Nektar::Operators;
 
 int main()
@@ -12,24 +16,49 @@ int main()
         {ShapeType::eQuadrilateral, {4, 4, 1}, 10},
         {ShapeType::eTriangle, {4, 4, 1}, 20}};
 
-    Field<double, FieldState::Coeff> in(blocks);
-    Field<double, FieldState::Phys> out(blocks);
+    auto in  = Field<double, FieldState::Coeff>::create(blocks);
+    auto out = Field<double, FieldState::Phys>::create(blocks);
 
     GetOperatorFactory<double>().PrintAvailableClasses();
 
-    double* x = in.GetStorage().GetPtr();
+    double *x      = in.GetStorage().GetCPUPtr();
     const size_t n = in.GetStorage().size();
-    for (size_t i = 0; i < n; ++i) {
+    for (size_t i = 0; i < n; ++i)
+    {
         x[i] = i;
     }
 
-    BwdTrans<>::create()->apply(in, out);
-   // BwdTrans<>::create("SumFac")->apply(in, out);
-    
-    double* y = out.GetStorage().GetPtr();
-    for (size_t i = 0; i < n; ++i) {
+    BwdTrans<>::create("")->apply(in, out);
+    // BwdTrans<>::create("SumFac")->apply(in, out);
+
+    double *y = out.GetStorage().GetCPUPtr();
+    for (size_t i = 0; i < n; ++i)
+    {
         std::cout << y[i] << std::endl;
     }
 
+#ifdef NEKTAR_USE_CUDA
+
+    // Test CUDA call
+    in  = Field<double, FieldState::Coeff>::create<MemoryRegionCUDA>(blocks);
+    out = Field<double, FieldState::Phys>::create<MemoryRegionCUDA>(blocks);
+    BwdTrans<>::create("CUDA")->apply(in, out);
+
+    // Using a CPU operator with a CUDA mem region should also work
+    BwdTrans<>::create("MatFree")->apply(in, out);
+
+    try
+    {
+        // It does not work the other way around though
+        in  = Field<double, FieldState::Coeff>::create(blocks);
+        out = Field<double, FieldState::Phys>::create(blocks);
+        BwdTrans<>::create("CUDA")->apply(in, out);
+    }
+    catch (std::exception const &e)
+    {
+        std::cout << e.what() << std::endl;
+    }
+#endif
+
     return 0;
 }