diff --git a/CHANGELOG.md b/CHANGELOG.md
index 97a54c1ab..6891e9a22 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
index 886cf9a80..680ea17b1 100644
--- a/benchmarks/CMakeLists.txt
+++ b/benchmarks/CMakeLists.txt
@@ -28,3 +28,4 @@ endfunction()
 
 add_subdirectory(fields)
 add_subdirectory(finiteVolume/cellCentred/operator)
+add_subdirectory(finiteVolume/cellCentred/interpolation)
diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/CMakeLists.txt b/benchmarks/finiteVolume/cellCentred/interpolation/CMakeLists.txt
new file mode 100644
index 000000000..c82f7b0aa
--- /dev/null
+++ b/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)
diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp
new file mode 100644
index 000000000..af5dfd859
--- /dev/null
+++ b/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)); };
+    }
+}
+
+}
diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp
new file mode 100644
index 000000000..1e3b53dfe
--- /dev/null
+++ b/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)); };
+    }
+}
+
+}
diff --git a/benchmarks/finiteVolume/cellCentred/operator/gaussGreenGrad.cpp b/benchmarks/finiteVolume/cellCentred/operator/gaussGreenGrad.cpp
new file mode 100644
index 000000000..3d0d2e522
--- /dev/null
+++ b/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)); };
+    }
+}
diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp
similarity index 53%
rename from test/finiteVolume/cellCentred/operator/divOperator.cpp
rename to benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp
index 241acd149..78b3bc797 100644
--- a/test/finiteVolume/cellCentred/operator/divOperator.cpp
+++ b/benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp
@@ -1,20 +1,19 @@
 // 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 {}),
@@ -22,29 +21,22 @@ TEST_CASE("DivOperator")
     );
 
     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)); };
     }
 }
diff --git a/include/NeoFOAM/core/primitives/label.hpp b/include/NeoFOAM/core/primitives/label.hpp
index 41323e81b..279fb5c2b 100644
--- a/include/NeoFOAM/core/primitives/label.hpp
+++ b/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;
+};
+
 }
diff --git a/include/NeoFOAM/core/primitives/scalar.hpp b/include/NeoFOAM/core/primitives/scalar.hpp
index 47e172263..f6f8f8b63 100644
--- a/include/NeoFOAM/core/primitives/scalar.hpp
+++ b/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
diff --git a/include/NeoFOAM/core/primitives/traits.hpp b/include/NeoFOAM/core/primitives/traits.hpp
new file mode 100644
index 000000000..7ec060d56
--- /dev/null
+++ b/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;
+};
+
+}
diff --git a/include/NeoFOAM/core/primitives/vector.hpp b/include/NeoFOAM/core/primitives/vector.hpp
index 1843adefd..a7b6833fd 100644
--- a/include/NeoFOAM/core/primitives/vector.hpp
+++ b/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
diff --git a/include/NeoFOAM/dsl/explicit.hpp b/include/NeoFOAM/dsl/explicit.hpp
index e2fe0297e..6fa5274d8 100644
--- a/include/NeoFOAM/dsl/explicit.hpp
+++ b/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));
 }
 
 
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp
index e74de98d5..382c54ceb 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp
+++ b/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:
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
index a0b4334b0..e4acbd2c4 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp
@@ -18,8 +18,11 @@
 namespace NeoFOAM::finiteVolume::cellCentred
 {
 
+/* @class SurfaceInterpolationFactor
+**
+*/
 class SurfaceInterpolationFactory :
-    public NeoFOAM::RuntimeSelectionFactory<
+    public RuntimeSelectionFactory<
         SurfaceInterpolationFactory,
         Parameters<const Executor&, const UnstructuredMesh&, Input>>
 {
@@ -47,13 +50,16 @@ class SurfaceInterpolationFactory :
 
     virtual ~SurfaceInterpolationFactory() {} // Virtual destructor
 
-    virtual void
-    interpolate(const VolumeField<scalar>& volField, ScalarSurfaceField& surfaceField) const = 0;
+    virtual void interpolate(const VolumeField<scalar>& src, SurfaceField<scalar>& dst) const = 0;
 
     virtual void interpolate(
-        const ScalarSurfaceField& faceFlux,
-        const VolumeField<scalar>& volField,
-        ScalarSurfaceField& surfaceField
+        const ScalarSurfaceField& flux, const VolumeField<scalar>& src, SurfaceField<scalar>& dst
+    ) const = 0;
+
+    virtual void interpolate(const VolumeField<Vector>& src, SurfaceField<Vector>& dst) const = 0;
+
+    virtual void interpolate(
+        const ScalarSurfaceField& flux, const VolumeField<Vector>& src, SurfaceField<Vector>& dst
     ) const = 0;
 
     // Pure virtual function for cloning
@@ -91,39 +97,50 @@ class SurfaceInterpolation
           interpolationKernel_(SurfaceInterpolationFactory::create(exec, mesh, input)) {};
 
 
-    void interpolate(const VolumeField<scalar>& volField, ScalarSurfaceField& surfaceField) const
+    void interpolate(const VolumeField<scalar>& src, ScalarSurfaceField& dst) const
     {
-        interpolationKernel_->interpolate(volField, surfaceField);
+        interpolationKernel_->interpolate(src, dst);
     }
 
-    ScalarSurfaceField interpolate(const VolumeField<scalar>& volField) const
+    void interpolate(const VolumeField<Vector>& src, SurfaceField<Vector>& dst) const
     {
-        std::string nameInterpolated = "interpolated_" + volField.name;
-        ScalarSurfaceField surfaceField(
-            exec_, nameInterpolated, mesh_, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh_)
-        );
-        interpolate(surfaceField, volField);
-        return surfaceField;
+        interpolationKernel_->interpolate(src, dst);
     }
 
+
     void interpolate(
-        const ScalarSurfaceField& faceFlux,
-        const VolumeField<scalar>& volField,
-        ScalarSurfaceField& surfaceField
+        const ScalarSurfaceField& flux, const VolumeField<scalar>& src, ScalarSurfaceField& dst
     ) const
     {
-        interpolationKernel_->interpolate(faceFlux, volField, surfaceField);
+        interpolationKernel_->interpolate(flux, src, dst);
+    }
+
+    void interpolate(
+        const ScalarSurfaceField& flux, const VolumeField<Vector>& src, SurfaceField<Vector>& dst
+    ) const
+    {
+        interpolationKernel_->interpolate(flux, src, dst);
+    }
+
+    ScalarSurfaceField interpolate(const VolumeField<scalar>& src) const
+    {
+        std::string nameInterpolated = "interpolated_" + src.name;
+        ScalarSurfaceField dst(
+            exec_, nameInterpolated, mesh_, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh_)
+        );
+        interpolate(dst, src);
+        return dst;
     }
 
     ScalarSurfaceField
-    interpolate(const ScalarSurfaceField& faceFlux, const VolumeField<scalar>& volField) const
+    interpolate(const ScalarSurfaceField& flux, const VolumeField<scalar>& src) const
     {
-        std::string name = "interpolated_" + volField.name;
-        ScalarSurfaceField surfaceField(
+        std::string name = "interpolated_" + src.name;
+        ScalarSurfaceField dst(
             exec_, name, mesh_, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh_)
         );
-        interpolate(faceFlux, volField, surfaceField);
-        return surfaceField;
+        interpolate(flux, src, dst);
+        return dst;
     }
 
 private:
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp
index aede0152b..1ba6e2678 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp
@@ -29,13 +29,16 @@ class Upwind : public SurfaceInterpolationFactory::Register<Upwind>
 
     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>& flux, const VolumeField<scalar>& src, SurfaceField<scalar>& dst
+    ) const override;
 
     void interpolate(
-        const SurfaceField<scalar>& faceFlux,
-        const VolumeField<scalar>& volField,
-        SurfaceField<scalar>& surfaceField
+        const SurfaceField<scalar>& flux, const VolumeField<Vector>& src, SurfaceField<Vector>& dst
     ) const override;
 
     std::unique_ptr<SurfaceInterpolationFactory> clone() const override;
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
index 2ed0fb16f..f88f53a34 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
@@ -1,5 +1,5 @@
 // SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2023 NeoFOAM authors
+// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors
 
 #pragma once
 
@@ -41,6 +41,9 @@ class DivOperatorFactory :
 
     virtual ~DivOperatorFactory() {} // Virtual destructor
 
+    // NOTE currently simple overloading is used here, because templating the virtual function
+    // does not work and we cant template the entire class because the static create function
+    // cannot access keyExistsOrError and table anymore.
     virtual void
     div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) = 0;
@@ -48,9 +51,19 @@ class DivOperatorFactory :
     virtual void
     div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) = 0;
 
+    virtual void
+    div(VolumeField<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
+    ) = 0;
+
+    virtual void
+    div(Field<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) = 0;
+
     virtual VolumeField<scalar>
     div(const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) = 0;
 
+    virtual VolumeField<Vector>
+    div(const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) = 0;
+
     // Pure virtual function for cloning
     virtual std::unique_ptr<DivOperatorFactory> clone() const = 0;
 
@@ -61,14 +74,15 @@ class DivOperatorFactory :
     const UnstructuredMesh& mesh_;
 };
 
-class DivOperator : public dsl::OperatorMixin<VolumeField<scalar>>
+template<typename ValueType>
+class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>
 {
 
 public:
 
     // copy constructor
     DivOperator(const DivOperator& divOp)
-        : dsl::OperatorMixin<VolumeField<scalar>>(divOp.exec_, divOp.field_, divOp.type_),
+        : dsl::OperatorMixin<VolumeField<ValueType>>(divOp.exec_, divOp.field_, divOp.type_),
           faceFlux_(divOp.faceFlux_),
           divOperatorStrategy_(
               divOp.divOperatorStrategy_ ? divOp.divOperatorStrategy_->clone() : nullptr
@@ -77,16 +91,17 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<scalar>>
     DivOperator(
         dsl::Operator::Type termType,
         const SurfaceField<scalar>& faceFlux,
-        VolumeField<scalar>& phi,
+        VolumeField<ValueType>& phi,
         Input input
     )
-        : dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), phi, termType), faceFlux_(faceFlux),
-          divOperatorStrategy_(DivOperatorFactory::create(exec_, phi.mesh(), input)) {};
+        : dsl::OperatorMixin<VolumeField<ValueType>>(phi.exec(), phi, termType),
+          faceFlux_(faceFlux),
+          divOperatorStrategy_(DivOperatorFactory::create(phi.exec(), phi.mesh(), input)) {};
 
     DivOperator(
         dsl::Operator::Type termType,
         const SurfaceField<scalar>& faceFlux,
-        VolumeField<scalar>& phi,
+        VolumeField<ValueType>& phi,
         std::unique_ptr<DivOperatorFactory> divOperatorStrategy
     )
         : dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), phi, termType), faceFlux_(faceFlux),
@@ -106,32 +121,35 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<scalar>>
             NF_ERROR_EXIT("DivOperatorStrategy not initialized");
         }
         NeoFOAM::Field<NeoFOAM::scalar> tmpsource(source.exec(), source.size(), 0.0);
-        divOperatorStrategy_->div(tmpsource, faceFlux_, field_);
+        divOperatorStrategy_->div(tmpsource, faceFlux_, this->getField());
         source += tmpsource;
     }
 
-    void div(Field<scalar>& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, getField()); }
+    void div(Field<ValueType>& divPhi)
+    {
+        divOperatorStrategy_->div(divPhi, faceFlux_, this->getField());
+    }
 
-    void div(VolumeField<scalar>& divPhi)
+    void div(VolumeField<ValueType>& divPhi)
     {
-        divOperatorStrategy_->div(divPhi, faceFlux_, getField());
+        divOperatorStrategy_->div(divPhi, faceFlux_, this->getField());
     }
 
 
     void build(const Input& input)
     {
-        const UnstructuredMesh& mesh = field_.mesh();
+        const UnstructuredMesh& mesh = this->getField().mesh();
         if (std::holds_alternative<NeoFOAM::Dictionary>(input))
         {
             auto dict = std::get<NeoFOAM::Dictionary>(input);
-            std::string schemeName = "div(" + faceFlux_.name + "," + field_.name + ")";
+            std::string schemeName = "div(" + faceFlux_.name + "," + this->getField().name + ")";
             auto tokens = dict.subDict("divSchemes").get<NeoFOAM::TokenList>(schemeName);
-            divOperatorStrategy_ = DivOperatorFactory::create(exec(), mesh, tokens);
+            divOperatorStrategy_ = DivOperatorFactory::create(this->exec(), mesh, tokens);
         }
         else
         {
             auto tokens = std::get<NeoFOAM::TokenList>(input);
-            divOperatorStrategy_ = DivOperatorFactory::create(exec(), mesh, tokens);
+            divOperatorStrategy_ = DivOperatorFactory::create(this->exec(), mesh, tokens);
         }
     }
 
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
index 6d0fdd9d7..9e0a2fa54 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
@@ -24,17 +24,28 @@ class GaussGreenDiv : public DivOperatorFactory::Register<GaussGreenDiv>
 
     GaussGreenDiv(const Executor& exec, const UnstructuredMesh& mesh, const Input& inputs);
 
-    void
+    virtual void
     div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) override;
 
-    void
+    virtual void
     div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
     ) override;
 
-    VolumeField<scalar>
+    virtual VolumeField<scalar>
     div(const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) override;
 
+    virtual void
+    div(VolumeField<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
+    ) override;
+
+    virtual void
+    div(Field<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
+    ) override;
+
+    virtual VolumeField<Vector>
+    div(const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) override;
+
     std::unique_ptr<DivOperatorFactory> clone() const override;
 
 private:
diff --git a/scripts/catch2json.py b/scripts/catch2json.py
index 9980b0937..4fefe5c1b 100644
--- a/scripts/catch2json.py
+++ b/scripts/catch2json.py
@@ -15,7 +15,12 @@ def parse_xml_dict(d):
         data = [data]
     records = []
     for cases in data:
-        test_case = cases["@name"]
+        test_case_raw_name = cases["@name"].split(" - ")
+        test_case = test_case_raw_name[0]
+        if len(test_case_raw_name) > 1:
+            test_type = test_case_raw_name[-1]
+        else:
+            test_type = ""
         for d in cases["Section"]:
             size = d["@name"]
             res = {}
@@ -33,6 +38,7 @@ def parse_xml_dict(d):
                     continue
             res["size"] = size
             res["test_case"] = test_case
+            res["test_type"] = test_type
             records.append(res)
     return records
 
diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp
index 9fc257852..37c7dddd0 100644
--- a/src/finiteVolume/cellCentred/interpolation/linear.cpp
+++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp
@@ -9,39 +9,47 @@
 namespace NeoFOAM::finiteVolume::cellCentred
 {
 
+/* @brief computional kernel to perform a linear interpolation
+** from a source volumeField to a surface field. It performs an interpolation
+** of the form
+**
+** d_f = w_f * s_O + ( 1 - w_f ) * s_N
+**
+**@param src the input field
+**@param weights weights for the interpolation
+**@param dst the target field
+*/
+template<typename ValueType>
 void computeLinearInterpolation(
-    const VolumeField<scalar>& volField,
-    const std::shared_ptr<GeometryScheme> geometryScheme,
-    SurfaceField<scalar>& surfaceField
+    const VolumeField<ValueType>& src,
+    const SurfaceField<scalar>& weights,
+    SurfaceField<ValueType>& dst
 )
 {
-    const UnstructuredMesh& mesh = surfaceField.mesh();
-    const auto& exec = surfaceField.exec();
-    auto sfield = surfaceField.internalField().span();
-    const NeoFOAM::labelField& owner = mesh.faceOwner();
-    const NeoFOAM::labelField& neighbour = mesh.faceNeighbour();
-
-    const auto sWeight = geometryScheme->weights().internalField().span();
-    const auto sVolField = volField.internalField().span();
-    const auto sBField = volField.boundaryField().value().span();
-    const auto sOwner = owner.span();
-    const auto sNeighbour = neighbour.span();
-    size_t nInternalFaces = mesh.nInternalFaces();
+    const auto exec = dst.exec();
+    auto dstS = dst.internalField().span();
+    const auto [srcS, weightS, ownerS, neighS, boundS] = spans(
+        src.internalField(),
+        weights.internalField(),
+        dst.mesh().faceOwner(),
+        dst.mesh().faceNeighbour(),
+        src.boundaryField().value()
+    );
+    size_t nInternalFaces = dst.mesh().nInternalFaces();
 
     NeoFOAM::parallelFor(
         exec,
-        {0, sfield.size()},
+        {0, dstS.size()},
         KOKKOS_LAMBDA(const size_t facei) {
-            size_t own = static_cast<size_t>(sOwner[facei]);
-            size_t nei = static_cast<size_t>(sNeighbour[facei]);
             if (facei < nInternalFaces)
             {
-                sfield[facei] =
-                    sWeight[facei] * sVolField[own] + (1 - sWeight[facei]) * sVolField[nei];
+                size_t own = static_cast<size_t>(ownerS[facei]);
+                size_t nei = static_cast<size_t>(neighS[facei]);
+                dstS[facei] = weightS[facei] * srcS[own] + (1 - weightS[facei]) * srcS[nei];
             }
             else
             {
-                sfield[facei] = sWeight[facei] * sBField[facei - nInternalFaces];
+                dstS[facei] = weightS[facei] * boundS[facei - nInternalFaces];
             }
         }
     );
@@ -55,21 +63,35 @@ Linear::Linear(const Executor& exec, const UnstructuredMesh& mesh)
     : SurfaceInterpolationFactory::Register<Linear>(exec, mesh),
       geometryScheme_(GeometryScheme::readOrCreate(mesh)) {};
 
-void Linear::interpolate(const VolumeField<scalar>& volField, SurfaceField<scalar>& surfaceField)
-    const
+void Linear::interpolate(const VolumeField<scalar>& src, SurfaceField<scalar>& dst) const
+{
+    computeLinearInterpolation(src, geometryScheme_->weights(), dst);
+}
+
+void Linear::interpolate(const VolumeField<Vector>& src, SurfaceField<Vector>& dst) const
 {
-    computeLinearInterpolation(volField, geometryScheme_, surfaceField);
+    computeLinearInterpolation(src, geometryScheme_->weights(), dst);
 }
 
 void Linear::interpolate(
-    [[maybe_unused]] const SurfaceField<scalar>& faceFlux,
-    const VolumeField<scalar>& volField,
-    SurfaceField<scalar>& surfaceField
+    [[maybe_unused]] const SurfaceField<scalar>& flux,
+    const VolumeField<Vector>& src,
+    SurfaceField<Vector>& dst
 ) const
 {
-    interpolate(volField, surfaceField);
+    interpolate(src, dst);
 }
 
+void Linear::interpolate(
+    [[maybe_unused]] const SurfaceField<scalar>& flux,
+    const VolumeField<scalar>& src,
+    SurfaceField<scalar>& dst
+) const
+{
+    interpolate(src, dst);
+}
+
+
 std::unique_ptr<SurfaceInterpolationFactory> Linear::clone() const
 {
     return std::make_unique<Linear>(*this);
diff --git a/src/finiteVolume/cellCentred/interpolation/upwind.cpp b/src/finiteVolume/cellCentred/interpolation/upwind.cpp
index e0d493a8c..85e7e1de2 100644
--- a/src/finiteVolume/cellCentred/interpolation/upwind.cpp
+++ b/src/finiteVolume/cellCentred/interpolation/upwind.cpp
@@ -9,47 +9,57 @@
 namespace NeoFOAM::finiteVolume::cellCentred
 {
 
+/* @brief computional kernel to perform an upwind interpolation
+** from a source volumeField to a surface field. It performs an interpolation
+** of the form
+**
+** d_f = w_f * s_O + ( 1 - w_f ) * s_N where w_1 is either 0,1 depending on the
+** direction of the flux
+**
+**@param src the input field
+**@param weights weights for the interpolation
+**@param dst the target field
+*/
+template<typename ValueType>
 void computeUpwindInterpolation(
-    const SurfaceField<scalar>& faceFlux,
-    const VolumeField<scalar>& volField,
-    const std::shared_ptr<GeometryScheme> geometryScheme,
-    SurfaceField<scalar>& surfaceField
+    const VolumeField<ValueType>& src,
+    const SurfaceField<scalar>& flux,
+    const SurfaceField<scalar>& weights,
+    SurfaceField<ValueType>& dst
 )
 {
-    const UnstructuredMesh& mesh = surfaceField.mesh();
-    const auto& exec = surfaceField.exec();
-
-    auto sfield = surfaceField.internalField().span();
-    const NeoFOAM::labelField& owner = mesh.faceOwner();
-    const NeoFOAM::labelField& neighbour = mesh.faceNeighbour();
-    const auto sWeight = geometryScheme->weights().internalField().span();
-    const auto sFaceFlux = faceFlux.internalField().span();
-    const auto sVolField = volField.internalField().span();
-    const auto sBField = volField.boundaryField().value().span();
-    const auto sOwner = owner.span();
-    const auto sNeighbour = neighbour.span();
-    size_t nInternalFaces = mesh.nInternalFaces();
+    const auto exec = dst.exec();
+    auto dstS = dst.internalField().span();
+    const auto [srcS, weightS, ownerS, neighS, boundS, fluxS] = spans(
+        src.internalField(),
+        weights.internalField(),
+        dst.mesh().faceOwner(),
+        dst.mesh().faceNeighbour(),
+        src.boundaryField().value(),
+        flux.internalField()
+    );
+    size_t nInternalFaces = dst.mesh().nInternalFaces();
 
     NeoFOAM::parallelFor(
         exec,
-        {0, sfield.size()},
+        {0, dstS.size()},
         KOKKOS_LAMBDA(const size_t facei) {
             if (facei < nInternalFaces)
             {
-                if (sFaceFlux[facei] >= 0)
+                if (fluxS[facei] >= 0)
                 {
-                    size_t own = static_cast<size_t>(sOwner[facei]);
-                    sfield[facei] = sVolField[own];
+                    size_t own = static_cast<size_t>(ownerS[facei]);
+                    dstS[facei] = srcS[own];
                 }
                 else
                 {
-                    size_t nei = static_cast<size_t>(sNeighbour[facei]);
-                    sfield[facei] = sVolField[nei];
+                    size_t nei = static_cast<size_t>(neighS[facei]);
+                    dstS[facei] = srcS[nei];
                 }
             }
             else
             {
-                sfield[facei] = sWeight[facei] * sBField[facei - nInternalFaces];
+                dstS[facei] = weightS[facei] * boundS[facei - nInternalFaces];
             }
         }
     );
@@ -60,20 +70,31 @@ Upwind::Upwind(const Executor& exec, const UnstructuredMesh& mesh, [[maybe_unuse
       geometryScheme_(GeometryScheme::readOrCreate(mesh)) {};
 
 void Upwind::interpolate(
-    [[maybe_unused]] const VolumeField<scalar>& volField,
-    [[maybe_unused]] SurfaceField<scalar>& surfaceField
+    [[maybe_unused]] const VolumeField<scalar>& src, [[maybe_unused]] SurfaceField<scalar>& dst
+) const
+{
+    NF_ERROR_EXIT("limited scheme require a faceFlux");
+}
+
+void Upwind::interpolate(
+    const SurfaceField<scalar>& flux, const VolumeField<scalar>& src, SurfaceField<scalar>& dst
+) const
+{
+    computeUpwindInterpolation(src, flux, geometryScheme_->weights(), dst);
+}
+
+void Upwind::interpolate(
+    [[maybe_unused]] const VolumeField<Vector>& src, [[maybe_unused]] SurfaceField<Vector>& dst
 ) const
 {
     NF_ERROR_EXIT("limited scheme require a faceFlux");
 }
 
 void Upwind::interpolate(
-    const SurfaceField<scalar>& faceFlux,
-    const VolumeField<scalar>& volField,
-    SurfaceField<scalar>& surfaceField
+    const SurfaceField<scalar>& flux, const VolumeField<Vector>& src, SurfaceField<Vector>& dst
 ) const
 {
-    computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField);
+    computeUpwindInterpolation(src, flux, geometryScheme_->weights(), dst);
 }
 
 std::unique_ptr<SurfaceInterpolationFactory> Upwind::clone() const
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
index 890bb3779..9ba87ffd0 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
@@ -8,51 +8,53 @@ namespace NeoFOAM::finiteVolume::cellCentred
 {
 
 
+/* @brief free standing function implementation of the divergence operator
+** ie computes 1/V \sum_f S_f \cdot \phi_f
+** where S_f is the face normal flux of a given face
+**  phi_f is the face interpolate value
+**
+**
+** @param faceFlux
+** @param neighbour - mapping from face id to neighbour cell id
+** @param owner - mapping from face id to owner cell id
+** @param faceCells - mapping from boundary face id to owner cell id
+*/
+template<typename ValueType>
 void computeDiv(
-    const SurfaceField<scalar>& faceFlux,
-    const VolumeField<scalar>& phi,
-    const SurfaceInterpolation& surfInterp,
-    Field<scalar>& divPhi
+    const Executor& exec,
+    size_t nInternalFaces,
+    size_t nBoundaryFaces,
+    std::span<const int> neighbour,
+    std::span<const int> owner,
+    std::span<const int> faceCells,
+    std::span<const scalar> faceFlux,
+    std::span<const ValueType> phiF,
+    std::span<const scalar> V,
+    std::span<ValueType> res
 )
 {
-    const UnstructuredMesh& mesh = phi.mesh();
-    const auto exec = phi.exec();
-    SurfaceField<scalar> phif(
-        exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh)
-    );
-    fill(phif.internalField(), 0.0);
-    const auto surfFaceCells = mesh.boundaryMesh().faceCells().span();
-    surfInterp.interpolate(faceFlux, phi, phif);
-
-    auto surfDivPhi = divPhi.span();
-
-    const auto surfPhif = phif.internalField().span();
-    const auto surfOwner = mesh.faceOwner().span();
-    const auto surfNeighbour = mesh.faceNeighbour().span();
-    const auto surfFaceFlux = faceFlux.internalField().span();
-    size_t nInternalFaces = mesh.nInternalFaces();
-    const auto surfV = mesh.cellVolumes().span();
-
+    size_t nCells {V.size()};
     // check if the executor is GPU
     if (std::holds_alternative<SerialExecutor>(exec))
     {
         for (size_t i = 0; i < nInternalFaces; i++)
         {
-            scalar flux = surfFaceFlux[i] * surfPhif[i];
-            surfDivPhi[static_cast<size_t>(surfOwner[i])] += flux;
-            surfDivPhi[static_cast<size_t>(surfNeighbour[i])] -= flux;
+            ValueType flux = faceFlux[i] * phiF[i];
+            res[static_cast<size_t>(owner[i])] += flux;
+            res[static_cast<size_t>(neighbour[i])] -= flux;
         }
 
-        for (size_t i = nInternalFaces; i < surfPhif.size(); i++)
+        for (size_t i = nInternalFaces; i < nInternalFaces + nBoundaryFaces; i++)
         {
-            auto own = static_cast<size_t>(surfFaceCells[i - nInternalFaces]);
-            scalar valueOwn = surfFaceFlux[i] * surfPhif[i];
-            surfDivPhi[own] += valueOwn;
+            auto own = static_cast<size_t>(faceCells[i - nInternalFaces]);
+            ValueType valueOwn = faceFlux[i] * phiF[i];
+            res[own] += valueOwn;
         }
 
-        for (size_t celli = 0; celli < mesh.nCells(); celli++)
+        // TODO does it make sense to store invVol and multiply?
+        for (size_t celli = 0; celli < nCells; celli++)
         {
-            surfDivPhi[celli] *= 1 / surfV[celli];
+            res[celli] *= 1 / V[celli];
         }
     }
     else
@@ -61,39 +63,62 @@ void computeDiv(
             exec,
             {0, nInternalFaces},
             KOKKOS_LAMBDA(const size_t i) {
-                scalar flux = surfFaceFlux[i] * surfPhif[i];
-                Kokkos::atomic_add(&surfDivPhi[static_cast<size_t>(surfOwner[i])], flux);
-                Kokkos::atomic_sub(&surfDivPhi[static_cast<size_t>(surfNeighbour[i])], flux);
+                ValueType flux = faceFlux[i] * phiF[i];
+                Kokkos::atomic_add(&res[static_cast<size_t>(owner[i])], flux);
+                Kokkos::atomic_sub(&res[static_cast<size_t>(neighbour[i])], flux);
             }
         );
 
         parallelFor(
             exec,
-            {nInternalFaces, surfPhif.size()},
+            {nInternalFaces, nInternalFaces + nBoundaryFaces},
             KOKKOS_LAMBDA(const size_t i) {
-                auto own = static_cast<size_t>(surfFaceCells[i - nInternalFaces]);
-                scalar valueOwn = surfFaceFlux[i] * surfPhif[i];
-                Kokkos::atomic_add(&surfDivPhi[own], valueOwn);
+                auto own = static_cast<size_t>(faceCells[i - nInternalFaces]);
+                ValueType valueOwn = faceFlux[i] * phiF[i];
+                Kokkos::atomic_add(&res[own], valueOwn);
             }
         );
 
         parallelFor(
-            exec,
-            {0, mesh.nCells()},
-            KOKKOS_LAMBDA(const size_t celli) { surfDivPhi[celli] *= 1 / surfV[celli]; }
+            exec, {0, nCells}, KOKKOS_LAMBDA(const size_t celli) { res[celli] *= 1 / V[celli]; }
         );
     }
 }
 
+
+template<typename ValueType>
 void computeDiv(
     const SurfaceField<scalar>& faceFlux,
-    const VolumeField<scalar>& phi,
+    const VolumeField<ValueType>& phi,
     const SurfaceInterpolation& surfInterp,
-    VolumeField<scalar>& divPhi
+    Field<ValueType>& divPhi
 )
 {
-    Field<scalar>& divPhiField = divPhi.internalField();
-    computeDiv(faceFlux, phi, surfInterp, divPhiField);
+    const UnstructuredMesh& mesh = phi.mesh();
+    const auto exec = phi.exec();
+    SurfaceField<ValueType> phif(
+        exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<ValueType>>(mesh)
+    );
+    // fill(phif.internalField(), NeoFOAM::zero<ValueType>::value);
+    surfInterp.interpolate(faceFlux, phi, phif);
+
+    // FIXME: currently we just copy the boundary values over
+    phif.boundaryField().value() = phi.boundaryField().value();
+
+    size_t nInternalFaces = mesh.nInternalFaces();
+    size_t nBoundaryFaces = mesh.nBoundaryFaces();
+    computeDiv<ValueType>(
+        exec,
+        nInternalFaces,
+        nBoundaryFaces,
+        mesh.faceNeighbour().span(),
+        mesh.faceOwner().span(),
+        mesh.boundaryMesh().faceCells().span(),
+        faceFlux.internalField().span(),
+        phif.internalField().span(),
+        mesh.cellVolumes().span(),
+        divPhi.span()
+    );
 }
 
 GaussGreenDiv::GaussGreenDiv(
@@ -106,14 +131,14 @@ void GaussGreenDiv::div(
     VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
 )
 {
-    computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi);
+    computeDiv<scalar>(faceFlux, phi, surfaceInterpolation_, divPhi.internalField());
 };
 
 void GaussGreenDiv::div(
     Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
 )
 {
-    computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi);
+    computeDiv<scalar>(faceFlux, phi, surfaceInterpolation_, divPhi);
 };
 
 VolumeField<scalar>
@@ -123,7 +148,33 @@ GaussGreenDiv::div(const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& ph
     VolumeField<scalar> divPhi(
         exec_, name, mesh_, createCalculatedBCs<VolumeBoundary<scalar>>(mesh_)
     );
-    computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi);
+    computeDiv<scalar>(faceFlux, phi, surfaceInterpolation_, divPhi.internalField());
+    return divPhi;
+};
+
+
+void GaussGreenDiv::div(
+    VolumeField<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
+)
+{
+    computeDiv<Vector>(faceFlux, phi, surfaceInterpolation_, divPhi.internalField());
+};
+
+void GaussGreenDiv::div(
+    Field<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
+)
+{
+    computeDiv<Vector>(faceFlux, phi, surfaceInterpolation_, divPhi);
+};
+
+VolumeField<Vector>
+GaussGreenDiv::div(const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi)
+{
+    std::string name = "div(" + faceFlux.name + "," + phi.name + ")";
+    VolumeField<Vector> divPhi(
+        exec_, name, mesh_, createCalculatedBCs<VolumeBoundary<Vector>>(mesh_)
+    );
+    computeDiv<Vector>(faceFlux, phi, surfaceInterpolation_, divPhi.internalField());
     return divPhi;
 };
 
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
index 319ae16c9..5a7954e04 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
@@ -8,35 +8,43 @@
 namespace NeoFOAM::finiteVolume::cellCentred
 {
 
+/* @brief free standing function implementation of the explicit gradient operator
+** ie computes \sum_f \phi_f
+**
+** @param[in] in - Field on which the gradient should be computed
+** @param[in,out] out - Field to hold the result
+*/
 void computeGrad(
-    const VolumeField<scalar>& phi,
-    const SurfaceInterpolation& surfInterp,
-    VolumeField<Vector>& gradPhi
+    const VolumeField<scalar>& in, const SurfaceInterpolation& surfInterp, VolumeField<Vector>& out
 )
 {
-    const UnstructuredMesh& mesh = gradPhi.mesh();
-    const auto exec = gradPhi.exec();
+    const UnstructuredMesh& mesh = out.mesh();
+    const auto exec = out.exec();
     SurfaceField<scalar> phif(
         exec, "phif", mesh, createCalculatedBCs<SurfaceBoundary<scalar>>(mesh)
     );
-    const auto surfFaceCells = mesh.boundaryMesh().faceCells().span();
-    const auto sBSf = mesh.boundaryMesh().sf().span();
-    auto surfGradPhi = gradPhi.internalField().span();
-
-    surfInterp.interpolate(phi, phif);
-    const auto surfPhif = phif.internalField().span();
-    const auto surfOwner = mesh.faceOwner().span();
-    const auto surfNeighbour = mesh.faceNeighbour().span();
-    const auto sSf = mesh.faceAreas().span();
+    surfInterp.interpolate(in, phif);
+
+    auto surfGradPhi = out.internalField().span();
+
+    const auto [surfFaceCells, sBSf, surfPhif, surfOwner, surfNeighbour, faceAreaS, surfV] = spans(
+        mesh.boundaryMesh().faceCells(),
+        mesh.boundaryMesh().sf(),
+        phif.internalField(),
+        mesh.faceOwner(),
+        mesh.faceNeighbour(),
+        mesh.faceAreas(),
+        mesh.cellVolumes()
+    );
+
     size_t nInternalFaces = mesh.nInternalFaces();
 
-    const auto surfV = mesh.cellVolumes().span();
 
     if (std::holds_alternative<SerialExecutor>(exec))
     {
         for (size_t i = 0; i < nInternalFaces; i++)
         {
-            Vector flux = sSf[i] * surfPhif[i];
+            Vector flux = faceAreaS[i] * surfPhif[i];
             surfGradPhi[static_cast<size_t>(surfOwner[i])] += flux;
             surfGradPhi[static_cast<size_t>(surfNeighbour[i])] -= flux;
         }
@@ -59,7 +67,7 @@ void computeGrad(
             exec,
             {0, nInternalFaces},
             KOKKOS_LAMBDA(const size_t i) {
-                Vector flux = sSf[i] * surfPhif[i];
+                Vector flux = faceAreaS[i] * surfPhif[i];
                 Kokkos::atomic_add(&surfGradPhi[static_cast<size_t>(surfOwner[i])], flux);
                 Kokkos::atomic_sub(&surfGradPhi[static_cast<size_t>(surfNeighbour[i])], flux);
             }
@@ -70,7 +78,7 @@ void computeGrad(
             {nInternalFaces, surfPhif.size()},
             KOKKOS_LAMBDA(const size_t i) {
                 size_t own = static_cast<size_t>(surfFaceCells[i - nInternalFaces]);
-                Vector valueOwn = sSf[i] * surfPhif[i];
+                Vector valueOwn = faceAreaS[i] * surfPhif[i];
                 Kokkos::atomic_add(&surfGradPhi[own], valueOwn);
             }
         );
diff --git a/test/core/primitives/CMakeLists.txt b/test/core/primitives/CMakeLists.txt
index ecd61c0b2..49ea4f15d 100644
--- a/test/core/primitives/CMakeLists.txt
+++ b/test/core/primitives/CMakeLists.txt
@@ -1,4 +1,5 @@
 # SPDX-License-Identifier: Unlicense
 # SPDX-FileCopyrightText: 2024 NeoFOAM authors
 
+neofoam_unit_test(scalar)
 neofoam_unit_test(vector)
diff --git a/test/core/primitives/scalar.cpp b/test/core/primitives/scalar.cpp
new file mode 100644
index 000000000..4b38084e9
--- /dev/null
+++ b/test/core/primitives/scalar.cpp
@@ -0,0 +1,35 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 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/core.hpp"
+
+TEST_CASE("Primitives")
+{
+    SECTION("Scalar", "[Traits]")
+    {
+        auto one = NeoFOAM::one<NeoFOAM::scalar>::value;
+
+        REQUIRE(one == 1.0);
+
+        auto zero = NeoFOAM::zero<NeoFOAM::scalar>::value;
+
+        REQUIRE(zero == 0.0);
+    }
+
+    SECTION("LocalIdx", "[Traits]")
+    {
+        auto one = NeoFOAM::one<NeoFOAM::localIdx>::value;
+
+        REQUIRE(one == 1);
+
+        auto zero = NeoFOAM::zero<NeoFOAM::localIdx>::value;
+
+        REQUIRE(zero == 0);
+    }
+}
diff --git a/test/core/primitives/vector.cpp b/test/core/primitives/vector.cpp
index 8fa3e86ca..0a07ee68d 100644
--- a/test/core/primitives/vector.cpp
+++ b/test/core/primitives/vector.cpp
@@ -11,6 +11,7 @@
 
 TEST_CASE("Primitives")
 {
+
     SECTION("Vector")
     {
         SECTION("CPU")
@@ -48,4 +49,20 @@ TEST_CASE("Primitives")
             REQUIRE((a + 2 * a + a) == d);
         }
     }
+
+    SECTION("Vector", "[Traits]")
+    {
+
+        auto one = NeoFOAM::one<NeoFOAM::Vector>::value;
+
+        REQUIRE(one(0) == 1.0);
+        REQUIRE(one(1) == 1.0);
+        REQUIRE(one(2) == 1.0);
+
+        auto zero = NeoFOAM::zero<NeoFOAM::Vector>::value;
+
+        REQUIRE(zero(0) == 0.0);
+        REQUIRE(zero(1) == 0.0);
+        REQUIRE(zero(2) == 0.0);
+    }
 }
diff --git a/test/finiteVolume/cellCentred/interpolation/CMakeLists.txt b/test/finiteVolume/cellCentred/interpolation/CMakeLists.txt
index 89f7c887a..b17cd174d 100644
--- a/test/finiteVolume/cellCentred/interpolation/CMakeLists.txt
+++ b/test/finiteVolume/cellCentred/interpolation/CMakeLists.txt
@@ -2,4 +2,5 @@
 # SPDX-FileCopyrightText: 2024 NeoFOAM authors
 
 neofoam_unit_test(linear)
+neofoam_unit_test(upwind)
 neofoam_unit_test(surfaceInterpolation)
diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp
index 64180c5d3..449c38cdd 100644
--- a/test/finiteVolume/cellCentred/interpolation/linear.cpp
+++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp
@@ -6,6 +6,7 @@
 #include <catch2/catch_session.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/generators/catch_generators_all.hpp>
+#include <catch2/catch_template_test_macros.hpp>
 
 #include "NeoFOAM/NeoFOAM.hpp"
 
@@ -13,7 +14,13 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation;
 using NeoFOAM::finiteVolume::cellCentred::VolumeField;
 using NeoFOAM::finiteVolume::cellCentred::SurfaceField;
 
-TEST_CASE("linear")
+namespace NeoFOAM
+{
+
+template<typename T>
+using I = std::initializer_list<T>;
+
+TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector)
 {
     NeoFOAM::Executor exec = GENERATE(
         NeoFOAM::Executor(NeoFOAM::SerialExecutor {}),
@@ -22,10 +29,40 @@ TEST_CASE("linear")
     );
 
     std::string execName = std::visit([](auto e) { return e.name(); }, exec);
-    auto mesh = NeoFOAM::createSingleCellMesh(exec);
-    NeoFOAM::Input input = NeoFOAM::TokenList({std::string("linear")});
+    auto mesh = create1DUniformMesh(exec, 10);
+    Input input = TokenList({std::string("linear")});
     auto linear = SurfaceInterpolation(exec, mesh, input);
+    std::vector<fvcc::VolumeBoundary<TestType>> vbcs {};
+    std::vector<fvcc::SurfaceBoundary<TestType>> sbcs {};
+    for (auto patchi : I<size_t> {0, 1})
+    {
+        Dictionary dict;
+        dict.insert("type", std::string("fixedValue"));
+        dict.insert("fixedValue", one<TestType>::value);
+        sbcs.push_back(fvcc::SurfaceBoundary<TestType>(mesh, dict, patchi));
+        vbcs.push_back(fvcc::VolumeBoundary<TestType>(mesh, dict, patchi));
+    }
+
+    auto in = VolumeField<TestType>(exec, "in", mesh, vbcs);
+    auto out = SurfaceField<TestType>(exec, "out", mesh, sbcs);
 
-    auto in = VolumeField<NeoFOAM::scalar>(exec, "in", mesh, {});
-    auto out = SurfaceField<NeoFOAM::scalar>(exec, "out", mesh, {});
+    fill(in.internalField(), one<TestType>::value);
+    in.correctBoundaryConditions();
+
+    linear.interpolate(in, out);
+    out.correctBoundaryConditions();
+
+    auto outHost = out.internalField().copyToHost();
+    auto nInternal = mesh.nInternalFaces();
+    auto nBoundary = mesh.nBoundaryFaces();
+    for (int i = 0; i < nInternal; i++)
+    {
+        REQUIRE(outHost[i] == one<TestType>::value);
+    }
+
+    for (int i = nInternal; i < nInternal + nBoundary; i++)
+    {
+        REQUIRE(outHost[i] == one<TestType>::value);
+    }
+}
 }
diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp
new file mode 100644
index 000000000..69f039382
--- /dev/null
+++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp
@@ -0,0 +1,68 @@
+// 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 <catch2/catch_session.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators_all.hpp>
+#include <catch2/catch_template_test_macros.hpp>
+
+#include "NeoFOAM/NeoFOAM.hpp"
+
+using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation;
+using NeoFOAM::finiteVolume::cellCentred::VolumeField;
+using NeoFOAM::finiteVolume::cellCentred::SurfaceField;
+
+namespace NeoFOAM
+{
+
+template<typename T>
+using I = std::initializer_list<T>;
+
+TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector)
+{
+    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);
+    auto mesh = create1DUniformMesh(exec, 10);
+    Input input = TokenList({std::string("upwind")});
+    auto upwind = SurfaceInterpolation(exec, mesh, input);
+    std::vector<fvcc::SurfaceBoundary<TestType>> bcs {};
+    for (auto patchi : I<size_t> {0, 1})
+    {
+        Dictionary dict;
+        dict.insert("type", std::string("fixedValue"));
+        dict.insert("fixedValue", one<TestType>::value);
+        bcs.push_back(fvcc::SurfaceBoundary<TestType>(mesh, dict, patchi));
+    }
+
+    auto in = VolumeField<TestType>(exec, "in", mesh, {});
+    auto flux = SurfaceField<scalar>(exec, "flux", mesh, {});
+    auto out = SurfaceField<TestType>(exec, "out", mesh, bcs);
+
+    fill(flux.internalField(), one<scalar>::value);
+    fill(in.internalField(), one<TestType>::value);
+
+    upwind.interpolate(flux, in, out);
+    out.correctBoundaryConditions();
+
+    auto outHost = out.internalField().copyToHost();
+    auto nInternal = mesh.nInternalFaces();
+    auto nBoundary = mesh.nBoundaryFaces();
+    for (int i = 0; i < nInternal; i++)
+    {
+        REQUIRE(outHost[i] == one<TestType>::value);
+    }
+
+    for (int i = nInternal; i < nInternal + nBoundary; i++)
+    {
+        REQUIRE(outHost[i] == one<TestType>::value);
+    }
+}
+
+}
diff --git a/test/finiteVolume/cellCentred/operator/CMakeLists.txt b/test/finiteVolume/cellCentred/operator/CMakeLists.txt
index 8d58bb29d..3b26396f4 100644
--- a/test/finiteVolume/cellCentred/operator/CMakeLists.txt
+++ b/test/finiteVolume/cellCentred/operator/CMakeLists.txt
@@ -1,4 +1,4 @@
 # SPDX-License-Identifier: Unlicense
 # SPDX-FileCopyrightText: 2024 NeoFOAM authors
 
-neofoam_unit_test(divOperator)
+neofoam_unit_test(gaussGreenDiv)
diff --git a/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp b/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp
new file mode 100644
index 000000000..a9ffe8101
--- /dev/null
+++ b/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp
@@ -0,0 +1,78 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2024 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 <catch2/catch_template_test_macros.hpp>
+
+#include "NeoFOAM/NeoFOAM.hpp"
+
+namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
+
+using Operator = NeoFOAM::dsl::Operator;
+
+namespace NeoFOAM
+{
+
+TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector)
+{
+    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);
+    auto mesh = create1DUniformMesh(exec, 10);
+    auto surfaceBCs = fvcc::createCalculatedBCs<fvcc::SurfaceBoundary<scalar>>(mesh);
+
+    // compute corresponding uniform faceFlux
+    // TODO this should be handled outside of the unit test
+    fvcc::SurfaceField<scalar> faceFlux(exec, "sf", mesh, surfaceBCs);
+    fill(faceFlux.internalField(), 1.0);
+    auto boundFaceFlux = faceFlux.internalField().span();
+    // face on the left side has different orientation
+    parallelFor(
+        exec,
+        {mesh.nCells() - 1, mesh.nCells()},
+        KOKKOS_LAMBDA(const size_t i) { boundFaceFlux[i] = -1.0; }
+    );
+
+    auto volumeBCs = fvcc::createCalculatedBCs<fvcc::VolumeBoundary<TestType>>(mesh);
+    fvcc::VolumeField<TestType> phi(exec, "sf", mesh, volumeBCs);
+    fill(phi.internalField(), one<TestType>::value);
+    fill(phi.boundaryField().value(), one<TestType>::value);
+    phi.correctBoundaryConditions();
+
+    auto result = Field<TestType>(exec, phi.size());
+    fill(result, zero<TestType>::value);
+
+    SECTION("Construct from Token" + execName)
+    {
+        Input input = TokenList({std::string("Gauss"), std::string("linear")});
+        fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input);
+    }
+
+    SECTION("Construct from Dictionary" + execName)
+    {
+        Input input = Dictionary(
+            {{std::string("DivOperator"), std::string("Gauss")},
+             {std::string("surfaceInterpolation"), std::string("linear")}}
+        );
+        auto op = fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input);
+        op.div(result);
+
+        // divergence of a uniform field should be zero
+        auto outHost = result.copyToHost();
+        for (int i = 0; i < result.size(); i++)
+        {
+            std::cout << "outHost[" << i << "]" << outHost[i] << "\n";
+            REQUIRE(outHost[i] == zero<TestType>::value);
+        }
+    }
+}
+
+}