Implementation of div operator for Field<Vector>

CHANGELOG.md

@@ -1,5 +1,6 @@
 # Version 0.2.0 (unreleased)
 ## Features
+- add vectorField implementation  []()
 - 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

benchmarks/CMakeLists.txt

@@ -28,3 +28,4 @@ endfunction()
 
 add_subdirectory(fields)
 add_subdirectory(finiteVolume/cellCentred/operator)
+add_subdirectory(finiteVolume/cellCentred/interpolation)

benchmarks/finiteVolume/cellCentred/interpolation/CMakeLists.txt

@@ -0,0 +1,5 @@
+# SPDX-License-Identifier: Unlicense
+# SPDX-FileCopyrightText: 2025 NeoFOAM authors
+
+neofoam_benchmark(linear)
+neofoam_benchmark(upwind)

benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2025 NeoFOAM authors
+
+#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
+                            // a custom main
+
+#include "NeoFOAM/NeoFOAM.hpp"
+#include "../../../catch_main.hpp"
+
+#include <catch2/catch_template_test_macros.hpp>
+
+using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation;
+using NeoFOAM::finiteVolume::cellCentred::VolumeField;
+using NeoFOAM::finiteVolume::cellCentred::SurfaceField;
+
+namespace NeoFOAM
+{
+
+TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector)
+{
+    auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20);
+
+    NeoFOAM::Executor exec = GENERATE(
+        NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
+    );
+
+    std::string execName = std::visit([](auto e) { return e.name(); }, exec);
+    UnstructuredMesh mesh = create1DUniformMesh(exec, size);
+    auto surfaceBCs = fvcc::createCalculatedBCs<fvcc::SurfaceBoundary<TestType>>(mesh);
+    Input input = TokenList({std::string("linear")});
+    auto linear = SurfaceInterpolation(exec, mesh, input);
+
+    auto in = VolumeField<TestType>(exec, "in", mesh, {});
+    auto out = SurfaceField<TestType>(exec, "out", mesh, surfaceBCs);
+
+    fill(in.internalField(), one<TestType>::value);
+
+    // capture the value of size as section name
+    DYNAMIC_SECTION("" << size)
+    {
+        BENCHMARK(std::string(execName)) { return (linear.interpolate(in, out)); };
+    }
+}
+
+}

benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp

@@ -0,0 +1,50 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2025 NeoFOAM authors
+
+#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
+                            // a custom main
+
+#include "NeoFOAM/NeoFOAM.hpp"
+#include "../../../catch_main.hpp"
+
+#include <catch2/catch_template_test_macros.hpp>
+
+using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation;
+using NeoFOAM::finiteVolume::cellCentred::VolumeField;
+using NeoFOAM::finiteVolume::cellCentred::SurfaceField;
+using NeoFOAM::Input;
+
+namespace NeoFOAM
+{
+
+TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector)
+{
+    auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20);
+
+    NeoFOAM::Executor exec = GENERATE(
+        NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
+    );
+
+    std::string execName = std::visit([](auto e) { return e.name(); }, exec);
+    UnstructuredMesh mesh = create1DUniformMesh(exec, size);
+    auto surfaceBCs = fvcc::createCalculatedBCs<fvcc::SurfaceBoundary<TestType>>(mesh);
+    Input input = TokenList({std::string("upwind")});
+    auto upwind = SurfaceInterpolation(exec, mesh, input);
+
+    auto in = VolumeField<TestType>(exec, "in", mesh, {});
+    auto flux = SurfaceField<scalar>(exec, "flux", mesh, {});
+    auto out = SurfaceField<TestType>(exec, "out", mesh, surfaceBCs);
+
+    fill(flux.internalField(), one<scalar>::value);
+    fill(in.internalField(), one<TestType>::value);
+
+    // capture the value of size as section name
+    DYNAMIC_SECTION("" << size)
+    {
+        BENCHMARK(std::string(execName)) { return (upwind.interpolate(flux, in, out)); };
+    }
+}
+
+}

benchmarks/finiteVolume/cellCentred/operator/gaussGreenGrad.cpp

@@ -0,0 +1,41 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+
+#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
+                            // a custom main
+
+#include "NeoFOAM/NeoFOAM.hpp"
+#include "../../../catch_main.hpp"
+
+using Operator = NeoFOAM::dsl::Operator;
+
+TEST_CASE("DivOperator::div", "[bench]")
+{
+    auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20);
+
+    NeoFOAM::Executor exec = GENERATE(
+        NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
+        NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
+    );
+
+    std::string execName = std::visit([](auto e) { return e.name(); }, exec);
+    NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size);
+    auto surfaceBCs = fvcc::createCalculatedBCs<fvcc::SurfaceBoundary<NeoFOAM::scalar>>(mesh);
+    fvcc::SurfaceField<NeoFOAM::scalar> faceFlux(exec, "sf", mesh, surfaceBCs);
+    NeoFOAM::fill(faceFlux.internalField(), 1.0);
+
+    auto volumeBCs = fvcc::createCalculatedBCs<fvcc::VolumeBoundary<NeoFOAM::scalar>>(mesh);
+    fvcc::VolumeField<NeoFOAM::scalar> phi(exec, "vf", mesh, volumeBCs);
+    fvcc::VolumeField<NeoFOAM::scalar> divPhi(exec, "divPhi", mesh, volumeBCs);
+    NeoFOAM::fill(phi.internalField(), 1.0);
+
+    // capture the value of size as section name
+    DYNAMIC_SECTION("" << size)
+    {
+        NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")});
+        auto op = fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input);
+
+        BENCHMARK(std::string(execName)) { return (op.div(divPhi)); };
+    }
+}

test/finiteVolume/cellCentred/operator/divOperator.cpp → benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp

@@ -1,50 +1,42 @@
 // SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2024 NeoFOAM authors
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
 
 #define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create
                             // a custom main
-#include <catch2/catch_session.hpp>
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/generators/catch_generators_all.hpp>
 
 #include "NeoFOAM/NeoFOAM.hpp"
-
-namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
+#include "../../../catch_main.hpp"
 
 using Operator = NeoFOAM::dsl::Operator;
 
-TEST_CASE("DivOperator")
+
+TEST_CASE("DivOperator::grad", "[bench]")
 {
+    auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20);
+
     NeoFOAM::Executor exec = GENERATE(
         NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
         NeoFOAM::Executor(NeoFOAM::CPUExecutor {}),
         NeoFOAM::Executor(NeoFOAM::GPUExecutor {})
     );
 
     std::string execName = std::visit([](auto e) { return e.name(); }, exec);
-    // TODO take 1d mesh
-    NeoFOAM::UnstructuredMesh mesh = NeoFOAM::createSingleCellMesh(exec);
+    NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size);
     auto surfaceBCs = fvcc::createCalculatedBCs<fvcc::SurfaceBoundary<NeoFOAM::scalar>>(mesh);
-
     fvcc::SurfaceField<NeoFOAM::scalar> faceFlux(exec, "sf", mesh, surfaceBCs);
     NeoFOAM::fill(faceFlux.internalField(), 1.0);
 
     auto volumeBCs = fvcc::createCalculatedBCs<fvcc::VolumeBoundary<NeoFOAM::scalar>>(mesh);
-    fvcc::VolumeField<NeoFOAM::scalar> phi(exec, "sf", mesh, volumeBCs);
+    fvcc::VolumeField<NeoFOAM::scalar> phi(exec, "vf", mesh, volumeBCs);
+    fvcc::VolumeField<NeoFOAM::scalar> divPhi(exec, "divPhi", mesh, volumeBCs);
     NeoFOAM::fill(phi.internalField(), 1.0);
 
-    SECTION("Construct from Token" + execName)
+    // capture the value of size as section name
+    DYNAMIC_SECTION("" << size)
     {
         NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")});
-        fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input);
-    }
+        auto op = fvcc::GradOperator(Operator::Type::Explicit, faceFlux, phi, input);
 
-    SECTION("Construct from Dictionary" + execName)
-    {
-        NeoFOAM::Input input = NeoFOAM::Dictionary(
-            {{std::string("DivOperator"), std::string("Gauss")},
-             {std::string("surfaceInterpolation"), std::string("linear")}}
-        );
-        fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input);
+        BENCHMARK(std::string(execName)) { return (op.grad(divPhi)); };
     }
 }

include/NeoFOAM/core/primitives/label.hpp

@@ -4,6 +4,9 @@
 
 #include <cstdint>
 
+#include "NeoFOAM/core/primitives/traits.hpp"
+
+
 namespace NeoFOAM
 {
 #ifdef NEOFOAM_DP_LABEL
@@ -16,4 +19,19 @@ using localIdx = uint32_t;
 using globalIdx = uint64_t;
 using size_t = std::size_t;
 using mpi_label_t = int;
+
+// traits for vector
+template<>
+struct one<localIdx>
+{
+    static const inline localIdx value = 1;
+};
+
+
+template<>
+struct zero<localIdx>
+{
+    static const inline localIdx value = 0;
+};
+
 }

include/NeoFOAM/core/primitives/scalar.hpp

@@ -2,6 +2,8 @@
 // SPDX-FileCopyrightText: 2023 NeoFOAM authors
 #pragma once
 
+#include "NeoFOAM/core/primitives/traits.hpp"
+
 // TODO this needs to be implemented in the corresponding cmake file
 namespace NeoFOAM
 {
@@ -13,4 +15,18 @@ typedef float scalar;
 
 constexpr scalar ROOTVSMALL = 1e-18;
 
+// traits for vector
+template<>
+struct one<scalar>
+{
+    static const inline scalar value = 1.0;
+};
+
+
+template<>
+struct zero<scalar>
+{
+    static const inline scalar value = 0.0;
+};
+
 } // namespace NeoFOAM

include/NeoFOAM/core/primitives/traits.hpp

@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+
+#pragma once
+
+namespace NeoFOAM
+{
+
+template<typename T>
+struct one
+{
+    static const T value;
+};
+
+template<typename T>
+struct zero
+{
+    static const T value;
+};
+
+}

include/NeoFOAM/core/primitives/vector.hpp

@@ -7,6 +7,8 @@
 
 #include "NeoFOAM/core/primitives/scalar.hpp"
 #include "NeoFOAM/core/primitives/label.hpp"
+#include "NeoFOAM/core/primitives/traits.hpp"
+
 
 namespace NeoFOAM
 {
@@ -153,4 +155,19 @@ scalar mag(const Vector& vec) { return sqrt(vec[0] * vec[0] + vec[1] * vec[1] +
 
 std::ostream& operator<<(std::ostream& out, const Vector& vec);
 
+
+// traits for vector
+template<>
+struct one<Vector>
+{
+    static const inline Vector value = {1.0, 1.0, 1.0};
+};
+
+
+template<>
+struct zero<Vector>
+{
+    static const inline Vector value = {0.0, 0.0, 0.0};
+};
+
 } // namespace NeoFOAM

include/NeoFOAM/dsl/explicit.hpp

@@ -18,10 +18,9 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
 namespace NeoFOAM::dsl::exp
 {
 
-Operator
-div(const fvcc::SurfaceField<NeoFOAM::scalar>& faceFlux, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+Operator div(const fvcc::SurfaceField<scalar>& faceFlux, fvcc::VolumeField<scalar>& phi)
 {
-    return Operator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi));
+    return Operator(fvcc::DivOperator<scalar>(dsl::Operator::Type::Explicit, faceFlux, phi));
 }
 
 

include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp

@@ -33,15 +33,19 @@ class Linear : public SurfaceInterpolationFactory::Register<Linear>
 
     static std::string schema() { return "none"; }
 
-    void interpolate(const VolumeField<scalar>& volField, SurfaceField<scalar>& surfaceField)
-        const override;
+    void interpolate(const VolumeField<scalar>& src, SurfaceField<scalar>& dst) const override;
+
+    void interpolate(const VolumeField<Vector>& src, SurfaceField<Vector>& dst) const override;
 
     void interpolate(
-        const SurfaceField<scalar>& faceFlux,
-        const VolumeField<scalar>& volField,
-        SurfaceField<scalar>& surfaceField
+        const SurfaceField<scalar>& flux, const VolumeField<scalar>& src, SurfaceField<scalar>& dst
     ) const override;
 
+    void interpolate(
+        const SurfaceField<scalar>& flux, const VolumeField<Vector>& src, SurfaceField<Vector>& dst
+    ) const override;
+
+
     std::unique_ptr<SurfaceInterpolationFactory> clone() const override;
 
 private: