Ginkgo config solver

cmake/CxxThirdParty.cmake

@@ -109,28 +109,22 @@ cpmaddpackage(
   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_MPI
-      OFF
-      CACHE INTERNAL "")
-  set(GINKGO_ENABLE_HALF
-      OFF
-      CACHE INTERNAL "")
   cpmaddpackage(
     NAME
     Ginkgo
+    VERSION
+    1.10.0
     GITHUB_REPOSITORY
     ginkgo-project/ginkgo
     GIT_TAG
     develop
+    OPTIONS
+    "GINKGO_BUILD_TESTS OFF"
+    "GINKGO_BUILD_BENCHMARKS OFF"
+    "GINKGO_BUILD_EXAMPLES OFF"
+    "GINKGO_BUILD_MPI ${NEOFOAM_ENABLE_MPI}"
+    "GINKGO_BUILD_CUDA ${Kokkos_ENABLE_CUDA}"
+    "GINKGO_BUILD_HIP ${Kokkos_ENABLE_HIP}"
     SYSTEM)
 endif()
 

include/NeoFOAM/core/dictionary.hpp

@@ -95,7 +95,7 @@ class Dictionary
         catch (const std::bad_any_cast& e)
         {
             logBadAnyCast<T>(e, key, data_);
-            throw e;
+            throw;
         }
     }
 
@@ -117,7 +117,7 @@ class Dictionary
         catch (const std::bad_any_cast& e)
         {
             logBadAnyCast<T>(e, key, data_);
-            throw e;
+            throw;
         }
     }
 

include/NeoFOAM/dsl/solver.hpp

@@ -112,7 +112,7 @@
     [[maybe_unused]] const Dictionary& fvSolution
 )
 {
     // FIXME:
     if (exp.temporalOperators().size() == 0 && exp.spatialOperators().size() == 0)
     {
         NF_ERROR_EXIT("No temporal or implicit terms to solve.");
@@ -132,8 +132,7 @@
         using ValueType = typename FieldType::ElementType;
         auto ls = ginkgoMatrix(exp, solution);
 
-
-        NeoFOAM::la::ginkgo::BiCGStab<ValueType> solver(solution.exec(), fvSolution);
+        NeoFOAM::la::ginkgo::Solver<ValueType> solver(solution.exec(), fvSolution);
         solver.solve(ls, solution.internalField());
     }
 }

include/NeoFOAM/linearAlgebra/ginkgo.hpp

@@ -17,98 +17,42 @@
 namespace NeoFOAM::la::ginkgo
 {
 
+gko::config::pnode parse(const Dictionary& dict);
+
 template<typename ValueType>
-class GkoSolverBase
+class Solver
 {
 
-private:
-
-    std::shared_ptr<const gko::Executor> gkoExec_;
-    Dictionary solverDict_;
-
-    virtual std::shared_ptr<gko::LinOp> solverGen(
-        std::shared_ptr<const gko::Executor> exec,
-        std::shared_ptr<const gko::LinOp> mtx,
-        size_t maxIter,
-        float relTol
-    ) = 0;
-
-protected:
+public:
 
-    GkoSolverBase(Executor exec, Dictionary solverDict)
-        : gkoExec_(getGkoExecutor(exec)), solverDict_(solverDict)
+    Solver(Executor exec, Dictionary solverConfig)
+        : gkoExec_(getGkoExecutor(exec)), config_(parse(solverConfig)),
+          factory_(
+              gko::config::parse(
+                  config_, gko::config::registry(), gko::config::make_type_descriptor<ValueType>()
+              )
+                  .on(gkoExec_)
+          )
     {}
 
-public:
-
-    virtual void solve(LinearSystem<ValueType, int>& sys, Field<ValueType>& x)
+    void solve(LinearSystem<ValueType, int>& sys, Field<ValueType>& x)
     {
         size_t nrows = sys.rhs().size();
 
-        auto solver = solverGen(
-            gkoExec_,
-            detail::createGkoMtx(gkoExec_, sys),
-            size_t(solverDict_.get<int>("maxIters")),
-            float(solverDict_.get<double>("relTol"))
-        );
+        auto solver = factory_->generate(detail::createGkoMtx(gkoExec_, sys));
 
         auto rhs = detail::createGkoDense(gkoExec_, sys.rhs().data(), nrows);
         auto gkoX = detail::createGkoDense(gkoExec_, x.data(), nrows);
         solver->apply(rhs, gkoX);
     }
-};
-
-
-template<typename ValueType>
-class CG : public GkoSolverBase<ValueType>
-{
-
-    virtual std::shared_ptr<gko::LinOp> solverGen(
-        std::shared_ptr<const gko::Executor> exec,
-        std::shared_ptr<const gko::LinOp> mtx,
-        size_t maxIter,
-        float relTol
-    ) override
-    {
-        auto fact =
-            gko::solver::Bicgstab<ValueType>::build()
-                .with_criteria(
-                    gko::stop::Iteration::build().with_max_iters(maxIter),
-                    gko::stop::ResidualNorm<ValueType>::build().with_reduction_factor(relTol)
-                )
-                .on(exec);
-        return fact->generate(mtx);
-    }
 
-public:
+private:
 
-    CG(Executor exec, Dictionary solverDict) : GkoSolverBase<ValueType>(exec, solverDict) {}
+    std::shared_ptr<const gko::Executor> gkoExec_;
+    gko::config::pnode config_;
+    std::shared_ptr<const gko::LinOpFactory> factory_;
 };
 
-template<typename ValueType>
-class BiCGStab : public GkoSolverBase<ValueType>
-{
-    virtual std::shared_ptr<gko::LinOp> solverGen(
-        std::shared_ptr<const gko::Executor> exec,
-        std::shared_ptr<const gko::LinOp> mtx,
-        size_t maxIter,
-        float relTol
-    )
-    {
-        auto fact =
-            gko::solver::Bicgstab<ValueType>::build()
-                .with_criteria(
-                    gko::stop::Iteration::build().with_max_iters(maxIter),
-                    gko::stop::ResidualNorm<ValueType>::build().with_reduction_factor(relTol)
-                )
-                .on(exec);
-        return fact->generate(mtx);
-    }
-
-public:
-
-    BiCGStab(Executor exec, Dictionary solverDict) : GkoSolverBase<ValueType>(exec, solverDict) {}
-};
 
 }
 

include/NeoFOAM/linearAlgebra/utilities.hpp

@@ -18,19 +18,33 @@ std::shared_ptr<gko::Executor> getGkoExecutor(Executor exec);
 
 namespace detail
 {
+
+template<typename T>
+gko::array<T> createGkoArray(std::shared_ptr<const gko::Executor> exec, std::span<T> values)
+{
+    return gko::make_array_view(exec, values.size(), values.data());
+}
+
+template<typename T>
+gko::detail::const_array_view<T>
+createGkoArray(std::shared_ptr<const gko::Executor> exec, std::span<const T> values)
+{
+    return gko::make_const_array_view(exec, values.size(), values.data());
+}
+
 template<typename ValueType, typename IndexType>
 std::shared_ptr<gko::matrix::Csr<ValueType, IndexType>>
 createGkoMtx(std::shared_ptr<const gko::Executor> exec, LinearSystem<ValueType, IndexType>& sys)
 {
     auto& mtx = sys.matrix();
     size_t nrows = sys.rhs().size();
 
-    auto valuesView = gko::array<scalar>::view(exec, mtx.values().size(), mtx.values().data());
-    auto colIdxView = gko::array<int>::view(exec, mtx.colIdxs().size(), mtx.colIdxs().data());
-    auto rowPtrView = gko::array<int>::view(exec, mtx.rowPtrs().size(), mtx.rowPtrs().data());
-
-    return gko::share(gko::matrix::Csr<scalar, int>::create(
-        exec, gko::dim<2> {nrows, nrows}, valuesView, colIdxView, rowPtrView
+    return gko::share(gko::matrix::Csr<ValueType, IndexType>::create(
+        exec,
+        gko::dim<2> {nrows, nrows},
+        createGkoArray(exec, mtx.values()),
+        createGkoArray(exec, mtx.colIdxs()),
+        createGkoArray(exec, mtx.rowPtrs())
     ));
 }
 
@@ -39,7 +53,7 @@ std::shared_ptr<gko::matrix::Dense<ValueType>>
 createGkoDense(std::shared_ptr<const gko::Executor> exec, ValueType* ptr, size_t size)
 {
     return gko::share(gko::matrix::Dense<ValueType>::create(
-        exec, gko::dim<2> {size, 1}, gko::array<scalar>::view(exec, size, ptr), 1
+        exec, gko::dim<2> {size, 1}, createGkoArray(exec, std::span {ptr, size}), 1
     ));
 }
 }

include/NeoFOAM/timeIntegration/backwardEuler.hpp

@@ -57,7 +57,7 @@ class BackwardEuler :
         auto ginkgoLs = NeoFOAM::dsl::ginkgoMatrix(ls, solutionField);
 
 
-        NeoFOAM::la::ginkgo::BiCGStab<ValueType> solver(solutionField.exec(), this->solutionDict_);
+        NeoFOAM::la::ginkgo::Solver<ValueType> solver(solutionField.exec(), this->solutionDict_);
         solver.solve(ginkgoLs, solutionField.internalField());
 
         // check if executor is GPU

src/CMakeLists.txt

@@ -28,6 +28,7 @@ target_sources(
           "executor/CPUExecutor.cpp"
           "executor/GPUExecutor.cpp"
           "executor/serialExecutor.cpp"
+          "linearAlgebra/ginkgo.cpp"
           "linearAlgebra/utilities.cpp"
           "mesh/unstructured/boundaryMesh.cpp"
           "mesh/unstructured/unstructuredMesh.cpp"

src/linearAlgebra/ginkgo.cpp

@@ -0,0 +1,66 @@
+// SPDX-FileCopyrightText: 2025 NeoFOAM authors
+//
+// SPDX-License-Identifier: MIT
+
+#include "NeoFOAM/linearAlgebra/ginkgo.hpp"
+
+gko::config::pnode NeoFOAM::la::ginkgo::parse(const Dictionary& dict)
+{
+    auto parseData = [&](auto key)
+    {
+        auto parseAny = [&](auto blueprint)
+        {
+            using value_type = decltype(blueprint);
+            if (dict[key].type() == typeid(value_type))
+            {
+                return gko::config::pnode(dict.get<value_type>(key));
+            }
+            else
+            {
+                return gko::config::pnode();
+            }
+        };
+
+        if (auto node = parseAny(std::string()))
+        {
+            return node;
+        }
+        if (auto node = parseAny(static_cast<const char*>(nullptr)))
+        {
+            return node;
+        }
+        if (auto node = parseAny(int {}))
+        {
+            return node;
+        }
+        if (auto node = parseAny((unsigned int) {}))
+        {
+            return node;
+        }
+        if (auto node = parseAny(double {}))
+        {
+            return node;
+        }
+        if (auto node = parseAny(float {}))
+        {
+            return node;
+        }
+
+        NF_THROW("Dictionary key " + key + " has unsupported type: " + dict[key].type().name());
+    };
+    gko::config::pnode::map_type result;
+    for (const auto& key : dict.keys())
+    {
+        gko::config::pnode node;
+        if (dict.isDict(key))
+        {
+            node = parse(dict.subDict(key));
+        }
+        else
+        {
+            node = parseData(key);
+        }
+        result.emplace(key, node);
+    }
+    return gko::config::pnode {result};
+}

test/linearAlgebra/CMakeLists.txt

@@ -1,5 +1,6 @@
 # SPDX-License-Identifier: Unlicense
 # SPDX-FileCopyrightText: 2024 NeoFOAM authors
 
+neofoam_unit_test(ginkgo)
 neofoam_unit_test(linearSystem)
 neofoam_unit_test(linearAlgebra)