Update/matrix assembly

CMakeLists.txt

@@ -36,7 +36,10 @@ option(NEOFOAM_DEVEL_TOOLS "Add development tools to the build system" OFF)
 option(NEOFOAM_BUILD_TESTS "Build the unit tests" OFF)
 option(NEOFOAM_BUILD_BENCHMARKS "Build benchmarks" OFF)
 option(NEOFOAM_BUILD_DOC "Build documentation" OFF)
+
 option(NEOFOAM_WITH_SUNDIALS "Build NeoFOAM with Sundials support [currently required]" ON)
+option(NEOFOAM_WITH_GINKGO "Enable using Ginkgo as linear algebra backend." ON)
+option(NEOFOAM_WITH_PETSC "Enable using PETSc as linear algebra backend." OFF)
 
 option(NEOFOAM_ENABLE_SANITIZE_ADDRESS "Enable address sanitizer" OFF)
 option(NEOFOAM_ENABLE_SANITIZE_LEAK "Enable leak sanitizer" OFF)

cmake/CxxThirdParty.cmake

@@ -108,6 +108,36 @@ cpmaddpackage(
   3.2.0
   SYSTEM)
 
+if(${NEOFOAM_WITH_GINKGO})
+  set(GINKGO_BUILD_TESTS
+      OFF
+      CACHE INTERNAL "")
+  set(GINKGO_BUILD_BENCHMARKS
+      OFF
+      CACHE INTERNAL "")
+  set(GINKGO_BUILD_EXAMPLES
+      OFF
+      CACHE INTERNAL "")
+  set(GINKGO_BUILD_EXAMPLES
+      OFF
+      CACHE INTERNAL "")
+  cpmaddpackage(
+    NAME
+    Ginkgo
+    GITHUB_REPOSITORY
+    ginkgo-project/ginkgo
+    GIT_TAG
+    develop
+    VERSION
+    1.9.0
+    SYSTEM)
+endif()
+
+if(${NEOFOAM_WITH_PETSC})
+  find_package(PkgConfig REQUIRED)
+  pkg_search_module(PETSc REQUIRED IMPORTED_TARGET PETSc)
+endif()
+
 if(NEOFOAM_BUILD_TESTS OR NEOFOAM_BUILD_BENCHMARKS)
   cpmaddpackage(
     NAME
@@ -118,3 +148,35 @@ if(NEOFOAM_BUILD_TESTS OR NEOFOAM_BUILD_BENCHMARKS)
     3.4.0
     SYSTEM)
 endif()
+
+if(${NEOFOAM_WITH_SUNDIALS})
+  set(SUNDIALS_OPTIONS
+      "BUILD_TESTING OFF"
+      "EXAMPLES_INSTALL OFF"
+      "BUILD_STATIC_LIBS OFF"
+      "BUILD_SHARED_LIBS ON"
+      "BUILD_ARKODE ON"
+      "BUILD_CVODE OFF"
+      "BUILD_CVODES OFF"
+      "BUILD_IDA OFF"
+      "BUILD_IDAS OFF"
+      "BUILD_KINSOL OFF"
+      "BUILD_CPODES OFF")
+  if(Kokkos_ENABLE_CUDA)
+    set(SUNDIALS_CUDA_OPTIONS "ENABLE_CUDA ON" "SUNDIALS_BUILD_KOKKOS ON")
+  else()
+    set(SUNDIALS_CUDA_OPTIONS "ENABLE_CUDA OFF" "SUNDIALS_BUILD_KOKKOS ON")
+  endif()
+
+  cpmaddpackage(
+    NAME
+    sundials
+    GITHUB_REPOSITORY
+    LLNL/sundials
+    VERSION
+    7.1.1
+    OPTIONS
+    ${SUNDIALS_OPTIONS}
+    ${SUNDIALS_CUDA_OPTIONS}
+    SYSTEM)
+endif()

include/CMakeLists.txt

@@ -13,6 +13,21 @@ set_property(
   APPEND
   PROPERTY INTERFACE_INCLUDE_DIRECTORIES "${CMAKE_CURRENT_BINARY_DIR}")
 
+if(NEOFOAM_WITH_GINKGO)
+  target_link_libraries(NeoFOAM_public_api INTERFACE Ginkgo::ginkgo)
+  target_compile_definitions(NeoFOAM_public_api INTERFACE NF_WITH_GINKGO=1 KOKKOS_ENABLE_SERIAL=1)
+  # set_property( TARGET Ginkgo::ginkgo APPEND PROPERTY INTERFACE_COMPILE_DEFINITIONS )
+else()
+  target_compile_definitions(NeoFOAM_public_api INTERFACE NF_WITH_GINKGO=0)
+endif()
+
+if(NEOFOAM_WITH_PETSC)
+  target_link_libraries(NeoFOAM_public_api INTERFACE PkgConfig::PETSc)
+  target_compile_definitions(NeoFOAM_public_api INTERFACE NF_WITH_PETSC=1)
+else()
+  target_compile_definitions(NeoFOAM_public_api INTERFACE NF_WITH_PETSC=0)
+endif()
+
 if(${CMAKE_BUILD_TYPE} MATCHES Debug)
   target_compile_definitions(NeoFOAM_public_api INTERFACE NF_DEBUG)
 endif()

include/NeoFOAM/NeoFOAM.hpp

@@ -1,9 +1,12 @@
 // SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
 
 #include "core.hpp"
 #include "dsl.hpp"
 #include "timeIntegration.hpp"
 #include "fields.hpp"
 
 #include "finiteVolume/cellCentred.hpp"
+
+#include "NeoFOAM/linearAlgebra/ginkgo.hpp"
+#include "linearAlgebra/petsc.hpp"

include/NeoFOAM/core/executor/CPUExecutor.hpp

@@ -58,6 +58,8 @@ class CPUExecutor
     void free(void* ptr) const noexcept { Kokkos::kokkos_free<exec>(ptr); };
 
     std::string name() const { return "CPUExecutor"; };
+
+    exec underlyingExec() const { return exec {}; }
 };
 
 } // namespace NeoFOAM

include/NeoFOAM/core/executor/GPUExecutor.hpp

@@ -59,6 +59,8 @@ class GPUExecutor
     void free(void* ptr) const noexcept { Kokkos::kokkos_free<exec>(ptr); }
 
     std::string name() const { return "GPUExecutor"; };
+
+    exec underlyingExec() const { return exec {}; }
 };
 
 } // namespace NeoFOAM

include/NeoFOAM/core/executor/serialExecutor.hpp

@@ -57,6 +57,8 @@ class SerialExecutor
     void free(void* ptr) const noexcept { Kokkos::kokkos_free<exec>(ptr); };
 
     std::string name() const { return "SerialExecutor"; };
+
+    exec underlyingExec() const { return exec {}; }
 };
 
 } // namespace NeoFOAM

include/NeoFOAM/linearAlgebra/ginkgo.hpp

@@ -0,0 +1,159 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2024 NeoFOAM authors
+
+#pragma once
+
+#if NF_WITH_GINKGO
+
+#include <ginkgo/ginkgo.hpp>
+#include <ginkgo/extensions/kokkos.hpp>
+
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/linearAlgebra/utilities.hpp"
+
+
+namespace NeoFOAM::la::ginkgo
+{
+
+std::shared_ptr<gko::Executor> getGkoExecutor(Executor exec)
+{
+    return std::visit(
+        [](auto concreteExec)
+        { return gko::ext::kokkos::create_executor(concreteExec.underlyingExec()); },
+        exec
+    );
+}
+
+template<typename ValueType, typename IndexType>
+class Matrix;
+
+template<typename ValueType, typename IndexType = int32_t>
+class MatrixBuilder
+{
+    friend class Matrix<ValueType, IndexType>;
+
+public:
+
+    struct SymbolicAssembly : CompatibleWithAnyExecutor
+    {
+        SymbolicAssembly(MatrixBuilder& mat)
+            : rowIdxs(mat.data.get_row_idxs()), colIdxs(mat.data.get_col_idxs())
+        {}
+
+        KOKKOS_FUNCTION void insert(size_t nnzId, MatrixCoordinate<IndexType> coordinate) const
+        {
+            rowIdxs[nnzId] = coordinate.row;
+            colIdxs[nnzId] = coordinate.col;
+        }
+
+    private:
+
+        IndexType* rowIdxs;
+        IndexType* colIdxs;
+    };
+
+    struct NumericAssembly : CompatibleWithAnyExecutor
+    {
+        NumericAssembly(MatrixBuilder& mat) : values(mat.data.get_values()) {}
+
+        KOKKOS_FUNCTION void insert(size_t nnzId, ValueType value) const { values[nnzId] = value; }
+
+    private:
+
+        ValueType* values;
+    };
+
+    struct Assembly : CompatibleWithAnyExecutor
+    {
+        Assembly(MatrixBuilder& mat) : symAssembly(mat), numAssembly(mat) {}
+
+        KOKKOS_FUNCTION void insert(size_t nnzId, MatrixEntry<ValueType, IndexType> entry) const
+        {
+            symAssembly.insert(nnzId, {entry.row, entry.col});
+            numAssembly.insert(nnzId, entry.value);
+        }
+
+    private:
+
+        SymbolicAssembly symAssembly;
+        NumericAssembly numAssembly;
+    };
+
+    MatrixBuilder(Executor exec, Dim dim, size_t nnzEstimation)
+        : exec_(exec),
+          data(getGkoExecutor(exec), gko::dim<2>(dim.numRows, dim.numCols), nnzEstimation)
+    {
+        auto size = data.get_size();
+        auto arrays = data.empty_out();
+        arrays.row_idxs.fill(0);
+        arrays.col_idxs.fill(0);
+        arrays.values.fill(ValueType {});
+        data = gko::device_matrix_data<ValueType, IndexType>(
+            data.get_executor(),
+            size,
+            std::move(arrays.row_idxs),
+            std::move(arrays.col_idxs),
+            std::move(arrays.values)
+        );
+    }
+
+    SymbolicAssembly startSymbolicAssembly() { return {*this}; }
+
+    NumericAssembly startNumericAssembly(SymbolicAssembly&& assembly) { return {*this}; }
+
+    Assembly startAssembly() { return {*this}; }
+
+    Executor getExecutor() const { return exec_; }
+
+private:
+
+    Executor exec_;
+    gko::device_matrix_data<ValueType, IndexType> data;
+};
+
+
+template<typename ValueType, typename IndexType = int32_t>
+class Matrix
+{
+public:
+
+    Matrix(MatrixBuilder<ValueType, IndexType>&& builder)
+        : exec_(builder.exec_),
+          mtx(gko::matrix::Coo<ValueType, IndexType>::create(builder.data.get_executor()))
+    {
+        builder.data.sum_duplicates();
+        builder.data.remove_zeros();
+        mtx->read(std::move(builder.data));
+    }
+
+    void apply(const Field<ValueType>& in, Field<ValueType>& out)
+    {
+        auto wrapField = [](auto& field)
+        {
+            auto data = const_cast<ValueType*>(field.data());
+            auto gkoExec = getGkoExecutor(field.exec());
+            return gko::matrix::Dense<ValueType>::create(
+                gkoExec,
+                gko::dim<2>(field.size(), 1),
+                gko::make_array_view(gkoExec, field.size(), data),
+                1
+            );
+        };
+        mtx->apply(wrapField(in), wrapField(out));
+    }
+
+    Executor getExecutor() const { return exec_; }
+
+    auto* getUnderlyingData() { return mtx.get(); }
+
+private:
+
+    Executor exec_;
+
+    std::shared_ptr<gko::matrix::Coo<ValueType, IndexType>> mtx;
+};
+
+
+}
+
+#endif