Merge branch 'stack/implicitOperators' into stack/vectorFields

.gitignore

@@ -12,6 +12,10 @@
 *.bak[0-9][0-9]
 \#*\#
 
+# clangd
+.cache
+compile_commands.json
+
 # File-browser settings - anywhere
 .directory
 

CHANGELOG.md

@@ -1,6 +1,9 @@
 # Version 0.2.0 (unreleased)
 ## Features
-- add vectorField implementation  []()
+- add basic vectorField implementation  [#260](https://github.com/exasim-project/NeoFOAM/pull/260)
+- add temporal operators inside expressions and the ability to solve linear system [#259 ](https://github.com/exasim-project/NeoFOAM/pull/259)
+- Add basic Ginkgo solver interface [#250](https://github.com/exasim-project/NeoFOAM/pull/250)
+- support for implicit operators in the DSL [#246](https://github.com/exasim-project/NeoFOAM/pull/246)
 - add segmentedField to represent vector of vectors [#202](https://github.com/exasim-project/NeoFOAM/pull/202)
 - Adds a minimal implementation linear algebra functionality [#219](https://github.com/exasim-project/NeoFOAM/pull/219)
 ## Fixes

CMakeLists.txt

@@ -36,7 +36,9 @@ 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_ENABLE_SANITIZE_ADDRESS "Enable address sanitizer" OFF)
 option(NEOFOAM_ENABLE_SANITIZE_LEAK "Enable leak sanitizer" OFF)

CMakePresets.json

@@ -11,7 +11,8 @@
       "hidden": true,
       "binaryDir": "${sourceDir}/build/${presetName}",
       "cacheVariables": {
-        "Kokkos_ENABLE_SERIAL": "ON"
+        "Kokkos_ENABLE_SERIAL": "ON",
+        "GINKGO_BUILD_PAPI_SDE": "OFF"
       },
       "generator": "Ninja"
     },

cmake/CxxThirdParty.cmake

@@ -108,6 +108,32 @@ 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_MPI
+      OFF
+      CACHE INTERNAL "")
+  set(GINKGO_ENABLE_HALF
+      OFF
+      CACHE INTERNAL "")
+  cpmaddpackage(
+    NAME
+    Ginkgo
+    GITHUB_REPOSITORY
+    ginkgo-project/ginkgo
+    GIT_TAG
+    develop
+    SYSTEM)
+endif()
+
 if(NEOFOAM_BUILD_TESTS OR NEOFOAM_BUILD_BENCHMARKS)
   cpmaddpackage(
     NAME

include/CMakeLists.txt

@@ -13,6 +13,14 @@ 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)
+  # set_property( TARGET Ginkgo::ginkgo APPEND PROPERTY INTERFACE_COMPILE_DEFINITIONS )
+else()
+  target_compile_definitions(NeoFOAM_public_api INTERFACE NF_WITH_GINKGO=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,10 @@
 // SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors
 
 #include "core.hpp"
 #include "dsl.hpp"
 #include "timeIntegration.hpp"
 #include "fields.hpp"
+#include "linearAlgebra.hpp"
 
 #include "finiteVolume/cellCentred.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/core/parallelAlgorithms.hpp

@@ -124,7 +124,7 @@ void parallelReduce(
     const NeoFOAM::Executor& exec, std::pair<size_t, size_t> range, Kernel kernel, T& value
 )
 {
-    return std::visit([&](const auto& e) { return parallelReduce(e, range, kernel, value); }, exec);
+    std::visit([&](const auto& e) { parallelReduce(e, range, kernel, value); }, exec);
 }
 
 
@@ -159,9 +159,7 @@ void parallelReduce(
 template<typename ValueType, typename Kernel, typename T>
 void parallelReduce(Field<ValueType>& field, Kernel kernel, T& value)
 {
-    return std::visit(
-        [&](const auto& e) { return parallelReduce(e, field, kernel, value); }, field.exec()
-    );
+    std::visit([&](const auto& e) { parallelReduce(e, field, kernel, value); }, field.exec());
 }
 
 template<typename Executor, typename Kernel>
@@ -203,9 +201,7 @@ void parallelScan(
     ReturnType& returnValue
 )
 {
-    return std::visit(
-        [&](const auto& e) { return parallelScan(e, range, kernel, returnValue); }, exec
-    );
+    std::visit([&](const auto& e) { parallelScan(e, range, kernel, returnValue); }, exec);
 }
 
 

include/NeoFOAM/dsl.hpp

@@ -6,5 +6,6 @@
 #include "dsl/explicit.hpp"
 #include "dsl/expression.hpp"
 #include "dsl/implicit.hpp"
-#include "dsl/operator.hpp"
+#include "dsl/spatialOperator.hpp"
+#include "dsl/temporalOperator.hpp"
 #include "dsl/solver.hpp"

include/NeoFOAM/dsl/ddt.hpp

@@ -5,7 +5,7 @@
 #pragma once
 
 #include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/dsl/operator.hpp"
+#include "NeoFOAM/dsl/spatialOperator.hpp"
 
 namespace NeoFOAM::dsl::temporal
 {
@@ -16,17 +16,21 @@ class Ddt : public OperatorMixin<FieldType>
 
 public:
 
-    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Temporal)
+    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Implicit)
     {}
 
     std::string getName() const { return "TimeOperator"; }
 
-    void explicitOperation([[maybe_unused]] Field<scalar>& source, [[maybe_unused]] scalar scale)
+    void explicitOperation(
+        [[maybe_unused]] Field<scalar>& source,
+        [[maybe_unused]] scalar t,
+        [[maybe_unused]] scalar dt
+    )
     {
         NF_ERROR_EXIT("Not implemented");
     }
 
-    void implicitOperation([[maybe_unused]] Field<scalar>& phi)
+    void implicitOperation( [[maybe_unused]] la::LinearSystem<scalar, localIdx>& ls,[[maybe_unused]] scalar t, [[maybe_unused]] scalar dt)
     {
         NF_ERROR_EXIT("Not implemented");
     }

include/NeoFOAM/dsl/explicit.hpp

@@ -5,22 +5,31 @@
 #pragma once
 
 #include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/dsl/operator.hpp"
+#include "NeoFOAM/dsl/spatialOperator.hpp"
 #include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
 #include "NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp"
 
 // TODO we should get rid of this include since it includes details
 // from a general implementation
 #include "NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp"
 
 namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
 
 namespace NeoFOAM::dsl::exp
 {
 
-Operator div(const fvcc::SurfaceField<scalar>& faceFlux, fvcc::VolumeField<scalar>& phi)
+
+SpatialOperator
+div(const fvcc::SurfaceField<NeoFOAM::scalar>& faceFlux, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return SpatialOperator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi));
+}
+
+SpatialOperator
+Source(fvcc::VolumeField<NeoFOAM::scalar>& coeff, fvcc::VolumeField<NeoFOAM::scalar>& phi)
 {
-    return Operator(fvcc::DivOperator<scalar>(dsl::Operator::Type::Explicit, faceFlux, phi));
+    return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Explicit, coeff, phi));
 }