From fd67a4d1bfe6d27cc02fe658187b86830e1f5b0e Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 11:36:12 +0100 Subject: [PATCH 01/35] update CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) 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 From 57a5ea253a7e05552d2a7ebd4999d28ff5b627a5 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 09:41:12 +0100 Subject: [PATCH 02/35] add basic vector traits --- include/NeoFOAM/core/primitives/traits.hpp | 20 ++++++++++++++++++++ include/NeoFOAM/core/primitives/vector.hpp | 17 +++++++++++++++++ test/core/primitives/vector.cpp | 16 ++++++++++++++++ 3 files changed, 53 insertions(+) create mode 100644 include/NeoFOAM/core/primitives/traits.hpp diff --git a/include/NeoFOAM/core/primitives/traits.hpp b/include/NeoFOAM/core/primitives/traits.hpp new file mode 100644 index 000000000..f7ee449b0 --- /dev/null +++ b/include/NeoFOAM/core/primitives/traits.hpp @@ -0,0 +1,20 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023 NeoFOAM authors + + +namespace NeoFOAM +{ +template +struct one +{ + static const T value = 1; +}; + +template +struct zero +{ + static const T value = 0; +}; + + +} 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 +{ + static const inline Vector value = {1.0, 1.0, 1.0}; +}; + + +template<> +struct zero +{ + static const inline Vector value = {0.0, 0.0, 0.0}; +}; + } // namespace NeoFOAM diff --git a/test/core/primitives/vector.cpp b/test/core/primitives/vector.cpp index 8fa3e86ca..1b9bde122 100644 --- a/test/core/primitives/vector.cpp +++ b/test/core/primitives/vector.cpp @@ -48,4 +48,20 @@ TEST_CASE("Primitives") REQUIRE((a + 2 * a + a) == d); } } + + SECTION("Traits") + { + + auto one = NeoFOAM::one::value; + + REQUIRE(one(0) == 1.0); + REQUIRE(one(1) == 1.0); + REQUIRE(one(2) == 1.0); + + auto zero = NeoFOAM::zero::value; + + REQUIRE(zero(0) == 0.0); + REQUIRE(zero(1) == 0.0); + REQUIRE(zero(2) == 0.0); + } } From 29e75f188eed4338f3dc2b44f300a20e57f00348 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 10:03:04 +0100 Subject: [PATCH 03/35] add basic scalar traits --- include/NeoFOAM/core/primitives/label.hpp | 18 ++++++ include/NeoFOAM/core/primitives/scalar.hpp | 16 ++++++ include/NeoFOAM/core/primitives/traits.hpp | 7 ++- test/core/primitives/CMakeLists.txt | 1 + test/core/primitives/scalar.cpp | 67 ++++++++++++++++++++++ test/core/primitives/vector.cpp | 3 +- 6 files changed, 108 insertions(+), 4 deletions(-) create mode 100644 test/core/primitives/scalar.cpp 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 +#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 +{ + static const inline localIdx value = 1; +}; + + +template<> +struct zero +{ + 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 +{ + static const inline scalar value = 1.0; +}; + + +template<> +struct zero +{ + 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 index f7ee449b0..7ec060d56 100644 --- a/include/NeoFOAM/core/primitives/traits.hpp +++ b/include/NeoFOAM/core/primitives/traits.hpp @@ -1,20 +1,21 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors +#pragma once namespace NeoFOAM { + template struct one { - static const T value = 1; + static const T value; }; template struct zero { - static const T value = 0; + static const T value; }; - } 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..1b9bde122 --- /dev/null +++ b/test/core/primitives/scalar.cpp @@ -0,0 +1,67 @@ +// 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 +#include +#include + +#include "NeoFOAM/core.hpp" + +TEST_CASE("Primitives") +{ + SECTION("Vector") + { + SECTION("CPU") + { + NeoFOAM::Vector a(1.0, 2.0, 3.0); + REQUIRE(a(0) == 1.0); + REQUIRE(a(1) == 2.0); + REQUIRE(a(2) == 3.0); + + NeoFOAM::Vector b(1.0, 2.0, 3.0); + REQUIRE(a == b); + + NeoFOAM::Vector c(2.0, 4.0, 6.0); + + REQUIRE(a + b == c); + + REQUIRE((a - b) == NeoFOAM::Vector(0.0, 0.0, 0.0)); + + a += b; + REQUIRE(a == c); + + a -= b; + REQUIRE(a == b); + a *= 2; + REQUIRE(a == c); + a = b; + + REQUIRE(a == b); + + NeoFOAM::Vector d(4.0, 8.0, 12.0); + REQUIRE((a + a + a + a) == d); + REQUIRE((4 * a) == d); + REQUIRE((a * 4) == d); + REQUIRE((a + 3 * a) == d); + REQUIRE((a + 2 * a + a) == d); + } + } + + SECTION("Traits") + { + + auto one = NeoFOAM::one::value; + + REQUIRE(one(0) == 1.0); + REQUIRE(one(1) == 1.0); + REQUIRE(one(2) == 1.0); + + auto zero = NeoFOAM::zero::value; + + REQUIRE(zero(0) == 0.0); + REQUIRE(zero(1) == 0.0); + REQUIRE(zero(2) == 0.0); + } +} diff --git a/test/core/primitives/vector.cpp b/test/core/primitives/vector.cpp index 1b9bde122..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") @@ -49,7 +50,7 @@ TEST_CASE("Primitives") } } - SECTION("Traits") + SECTION("Vector", "[Traits]") { auto one = NeoFOAM::one::value; From be5df76b42a3b806e560b94f8d513be00203a54e Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 13:09:03 +0100 Subject: [PATCH 04/35] clean-up --- test/core/primitives/scalar.cpp | 52 +++++++-------------------------- 1 file changed, 10 insertions(+), 42 deletions(-) diff --git a/test/core/primitives/scalar.cpp b/test/core/primitives/scalar.cpp index 1b9bde122..4b38084e9 100644 --- a/test/core/primitives/scalar.cpp +++ b/test/core/primitives/scalar.cpp @@ -11,57 +11,25 @@ TEST_CASE("Primitives") { - SECTION("Vector") + SECTION("Scalar", "[Traits]") { - SECTION("CPU") - { - NeoFOAM::Vector a(1.0, 2.0, 3.0); - REQUIRE(a(0) == 1.0); - REQUIRE(a(1) == 2.0); - REQUIRE(a(2) == 3.0); + auto one = NeoFOAM::one::value; - NeoFOAM::Vector b(1.0, 2.0, 3.0); - REQUIRE(a == b); + REQUIRE(one == 1.0); - NeoFOAM::Vector c(2.0, 4.0, 6.0); + auto zero = NeoFOAM::zero::value; - REQUIRE(a + b == c); - - REQUIRE((a - b) == NeoFOAM::Vector(0.0, 0.0, 0.0)); - - a += b; - REQUIRE(a == c); - - a -= b; - REQUIRE(a == b); - a *= 2; - REQUIRE(a == c); - a = b; - - REQUIRE(a == b); - - NeoFOAM::Vector d(4.0, 8.0, 12.0); - REQUIRE((a + a + a + a) == d); - REQUIRE((4 * a) == d); - REQUIRE((a * 4) == d); - REQUIRE((a + 3 * a) == d); - REQUIRE((a + 2 * a + a) == d); - } + REQUIRE(zero == 0.0); } - SECTION("Traits") + SECTION("LocalIdx", "[Traits]") { + auto one = NeoFOAM::one::value; - auto one = NeoFOAM::one::value; - - REQUIRE(one(0) == 1.0); - REQUIRE(one(1) == 1.0); - REQUIRE(one(2) == 1.0); + REQUIRE(one == 1); - auto zero = NeoFOAM::zero::value; + auto zero = NeoFOAM::zero::value; - REQUIRE(zero(0) == 0.0); - REQUIRE(zero(1) == 0.0); - REQUIRE(zero(2) == 0.0); + REQUIRE(zero == 0); } } From 816679f403598fc3fc3d232575c01dced2813777 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 13:05:04 +0100 Subject: [PATCH 05/35] wip add vector field interpolation --- .../cellCentred/interpolation/linear.hpp | 10 +++++++++ .../interpolation/surfaceInterpolation.hpp | 18 +++++++++++++--- .../cellCentred/interpolation/upwind.hpp | 9 ++++++++ .../cellCentred/interpolation/linear.cpp | 20 +++++++++++++++++- .../cellCentred/interpolation/upwind.cpp | 21 ++++++++++++++++++- .../cellCentred/interpolation/linear.cpp | 4 +++- 6 files changed, 76 insertions(+), 6 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp index e74de98d5..28bb97e6a 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp @@ -42,6 +42,16 @@ class Linear : public SurfaceInterpolationFactory::Register SurfaceField& surfaceField ) const override; + void interpolate(const VolumeField& volField, SurfaceField& surfaceField) + const override; + + void interpolate( + const SurfaceField& faceFlux, + const VolumeField& volField, + SurfaceField& surfaceField + ) const override; + + std::unique_ptr clone() const override; private: diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp index a0b4334b0..2efaa33d1 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> { @@ -48,12 +51,21 @@ class SurfaceInterpolationFactory : virtual ~SurfaceInterpolationFactory() {} // Virtual destructor virtual void - interpolate(const VolumeField& volField, ScalarSurfaceField& surfaceField) const = 0; + interpolate(const VolumeField& volField, SurfaceField& surfaceField) const = 0; virtual void interpolate( const ScalarSurfaceField& faceFlux, const VolumeField& volField, - ScalarSurfaceField& surfaceField + SurfaceField& surfaceField + ) const = 0; + + virtual void + interpolate(const VolumeField& volField, SurfaceField& surfaceField) const = 0; + + virtual void interpolate( + const ScalarSurfaceField& faceFlux, + const VolumeField& volField, + SurfaceField& surfaceField ) const = 0; // Pure virtual function for cloning diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp index aede0152b..d67603cd5 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp @@ -38,6 +38,15 @@ class Upwind : public SurfaceInterpolationFactory::Register SurfaceField& surfaceField ) const override; + void interpolate(const VolumeField& volField, SurfaceField& surfaceField) + const override; + + void interpolate( + const SurfaceField& faceFlux, + const VolumeField& volField, + SurfaceField& surfaceField + ) const override; + std::unique_ptr clone() const override; private: diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index 9fc257852..b0f6893f7 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -58,7 +58,8 @@ Linear::Linear(const Executor& exec, const UnstructuredMesh& mesh) void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) const { - computeLinearInterpolation(volField, geometryScheme_, surfaceField); + // FIXME: + // computeLinearInterpolation(volField, geometryScheme_, surfaceField); } void Linear::interpolate( @@ -70,6 +71,23 @@ void Linear::interpolate( interpolate(volField, surfaceField); } +void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) + const +{ + // FIXME: + // computeLinearInterpolation(volField, geometryScheme_, surfaceField); +} + +void Linear::interpolate( + [[maybe_unused]] const SurfaceField& faceFlux, + const VolumeField& volField, + SurfaceField& surfaceField +) const +{ + interpolate(volField, surfaceField); +} + + std::unique_ptr Linear::clone() const { return std::make_unique(*this); diff --git a/src/finiteVolume/cellCentred/interpolation/upwind.cpp b/src/finiteVolume/cellCentred/interpolation/upwind.cpp index e0d493a8c..d875fd28b 100644 --- a/src/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/src/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -73,7 +73,26 @@ void Upwind::interpolate( SurfaceField& surfaceField ) const { - computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField); + // FIXME: + // computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField); +} + +void Upwind::interpolate( + [[maybe_unused]] const VolumeField& volField, + [[maybe_unused]] SurfaceField& surfaceField +) const +{ + NF_ERROR_EXIT("limited scheme require a faceFlux"); +} + +void Upwind::interpolate( + const SurfaceField& faceFlux, + const VolumeField& volField, + SurfaceField& surfaceField +) const +{ + // FIXME: + // computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField); } std::unique_ptr Upwind::clone() const diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 64180c5d3..8f8c1bcbf 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -22,10 +22,12 @@ TEST_CASE("linear") ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - auto mesh = NeoFOAM::createSingleCellMesh(exec); + auto mesh = NeoFOAM::create1DUniformMesh(exec, 10); NeoFOAM::Input input = NeoFOAM::TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, {}); + + linear.interpolate(in, out); } From 1a102e2f007954f66316dafa1a5b9da54e2c4c50 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 22:14:45 +0100 Subject: [PATCH 06/35] template interpolation test --- .../interpolation/surfaceInterpolation.hpp | 6 ++++++ .../cellCentred/interpolation/linear.cpp | 13 ++++++++++--- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp index 2efaa33d1..682e053a0 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp @@ -108,6 +108,12 @@ class SurfaceInterpolation interpolationKernel_->interpolate(volField, surfaceField); } + void interpolate(const VolumeField& volField, SurfaceField& surfaceField) const + { + interpolationKernel_->interpolate(volField, surfaceField); + } + + ScalarSurfaceField interpolate(const VolumeField& volField) const { std::string nameInterpolated = "interpolated_" + volField.name; diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 8f8c1bcbf..44fd3e82e 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include "NeoFOAM/NeoFOAM.hpp" @@ -13,7 +14,7 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; -TEST_CASE("linear") +TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) { NeoFOAM::Executor exec = GENERATE( NeoFOAM::Executor(NeoFOAM::SerialExecutor {}), @@ -26,8 +27,14 @@ TEST_CASE("linear") NeoFOAM::Input input = NeoFOAM::TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); - auto in = VolumeField(exec, "in", mesh, {}); - auto out = SurfaceField(exec, "out", mesh, {}); + auto in = VolumeField(exec, "in", mesh, {}); + auto out = SurfaceField(exec, "out", mesh, {}); + + // FIXME add fill + // fill(in.internalField(), 1); linear.interpolate(in, out); + + // FIXME add fill + // REQUIRE(out.internalField()[0] == one::value); } From c8f96e48b90098b5064bc83c096e633a8bbd3c3a Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 14:07:39 +0100 Subject: [PATCH 07/35] make it work with scalar and vector, test fails --- .../cellCentred/interpolation/linear.cpp | 15 +++++++++------ .../cellCentred/interpolation/linear.cpp | 9 +++++---- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index b0f6893f7..ed2236174 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -9,10 +9,11 @@ namespace NeoFOAM::finiteVolume::cellCentred { +template void computeLinearInterpolation( - const VolumeField& volField, + const VolumeField& volField, const std::shared_ptr geometryScheme, - SurfaceField& surfaceField + SurfaceField& surfaceField ) { const UnstructuredMesh& mesh = surfaceField.mesh(); @@ -36,8 +37,12 @@ void computeLinearInterpolation( size_t nei = static_cast(sNeighbour[facei]); if (facei < nInternalFaces) { + std::cout << __FILE__ << " sWeight" << sWeight[facei] << "\n"; + std::cout << __FILE__ << " sVolfield[own]" << sVolField[own] << "\n"; + std::cout << __FILE__ << " sVolfield[nei]" << sVolField[nei] << "\n"; sfield[facei] = sWeight[facei] * sVolField[own] + (1 - sWeight[facei]) * sVolField[nei]; + std::cout << __FILE__ << " sfield[facei]" << sfield[facei] << "\n"; } else { @@ -58,8 +63,7 @@ Linear::Linear(const Executor& exec, const UnstructuredMesh& mesh) void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) const { - // FIXME: - // computeLinearInterpolation(volField, geometryScheme_, surfaceField); + computeLinearInterpolation(volField, geometryScheme_, surfaceField); } void Linear::interpolate( @@ -74,8 +78,7 @@ void Linear::interpolate( void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) const { - // FIXME: - // computeLinearInterpolation(volField, geometryScheme_, surfaceField); + computeLinearInterpolation(volField, geometryScheme_, surfaceField); } void Linear::interpolate( diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 44fd3e82e..b0597b142 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -30,11 +30,12 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, {}); - // FIXME add fill - // fill(in.internalField(), 1); + fill(in.internalField(), NeoFOAM::one::value); linear.interpolate(in, out); - // FIXME add fill - // REQUIRE(out.internalField()[0] == one::value); + for (int i = 0; i < out.internalField().size(); i++) + { + REQUIRE(out.internalField()[i] == NeoFOAM::one::value); + } } From 305d3329b8d16c85dc68590771e765cad3acb643 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 14:43:14 +0100 Subject: [PATCH 08/35] fix failing test, bcs were missing --- .../cellCentred/interpolation/linear.cpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index b0597b142..6c1e482e1 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -14,6 +14,9 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; +template +using I = std::initializer_list; + TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) { NeoFOAM::Executor exec = GENERATE( @@ -26,13 +29,22 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) auto mesh = NeoFOAM::create1DUniformMesh(exec, 10); NeoFOAM::Input input = NeoFOAM::TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); + std::vector> bcs {}; + for (auto patchi : I {0, 1}) + { + NeoFOAM::Dictionary dict; + dict.insert("type", std::string("fixedValue")); + dict.insert("fixedValue", NeoFOAM::one::value); + bcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); + } auto in = VolumeField(exec, "in", mesh, {}); - auto out = SurfaceField(exec, "out", mesh, {}); + auto out = SurfaceField(exec, "out", mesh, bcs); fill(in.internalField(), NeoFOAM::one::value); linear.interpolate(in, out); + out.correctBoundaryConditions(); for (int i = 0; i < out.internalField().size(); i++) { From 6f4e157f31ed8cf9f45f3e57af6fab3bb76481fb Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 16:13:01 +0100 Subject: [PATCH 09/35] add linear benchmark --- .../cellCentred/interpolation/linear.cpp | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp new file mode 100644 index 000000000..2f4d62c00 --- /dev/null +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -0,0 +1,39 @@ +// 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 NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; +using NeoFOAM::finiteVolume::cellCentred::VolumeField; +using NeoFOAM::finiteVolume::cellCentred::SurfaceField; + +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>(mesh); + + auto in = VolumeField(exec, "in", mesh, {}); + auto out = SurfaceField(exec, "out", mesh, surfaceBCs); + + fill(in.internalField(), NeoFOAM::one::value); + + // capture the value of size as section name + DYNAMIC_SECTION("" << size) + { + + BENCHMARK(std::string(execName)) {return (linear.interpolate(in, out); }; + } +} From 1fc3625d5a9233e500a52d7c9ee008daabaf2ed1 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 16:13:29 +0100 Subject: [PATCH 10/35] fixup! add linear benchmark --- benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp index 2f4d62c00..c551d6355 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -34,6 +34,6 @@ TEST_CASE("DivOperator::div", "[bench]") DYNAMIC_SECTION("" << size) { - BENCHMARK(std::string(execName)) {return (linear.interpolate(in, out); }; + BENCHMARK(std::string(execName)) { return (linear.interpolate(in, out)); }; } } From dac7479b3fa8cd9c334624e905a1d4a2174aa87c Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 17:42:00 +0100 Subject: [PATCH 11/35] wip implement upwind --- .../cellCentred/interpolation/linear.hpp | 16 ++-- .../interpolation/surfaceInterpolation.hpp | 45 ++++++----- .../cellCentred/interpolation/upwind.hpp | 16 ++-- .../cellCentred/interpolation/linear.cpp | 81 +++++++++---------- .../cellCentred/interpolation/upwind.cpp | 76 ++++++++--------- .../cellCentred/interpolation/CMakeLists.txt | 1 + .../cellCentred/interpolation/linear.cpp | 1 + .../cellCentred/interpolation/upwind.cpp | 56 +++++++++++++ 8 files changed, 171 insertions(+), 121 deletions(-) create mode 100644 test/finiteVolume/cellCentred/interpolation/upwind.cpp diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp index 28bb97e6a..382c54ceb 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp @@ -33,22 +33,16 @@ class Linear : public SurfaceInterpolationFactory::Register static std::string schema() { return "none"; } - void interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const override; + void interpolate(const VolumeField& src, SurfaceField& dst) const override; + + void interpolate(const VolumeField& src, SurfaceField& dst) const override; void interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const override; - void interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const override; - void interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const override; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp index 682e053a0..01028121d 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp @@ -103,45 +103,50 @@ class SurfaceInterpolation interpolationKernel_(SurfaceInterpolationFactory::create(exec, mesh, input)) {}; - void interpolate(const VolumeField& volField, ScalarSurfaceField& surfaceField) const + void interpolate(const VolumeField& src, ScalarSurfaceField& dst) const { - interpolationKernel_->interpolate(volField, surfaceField); + interpolationKernel_->interpolate(src, dst); } - void interpolate(const VolumeField& volField, SurfaceField& surfaceField) const + void interpolate(const VolumeField& src, SurfaceField& dst) const { - interpolationKernel_->interpolate(volField, surfaceField); + interpolationKernel_->interpolate(src, dst); } - ScalarSurfaceField interpolate(const VolumeField& volField) const + void interpolate( + const ScalarSurfaceField& flux, const VolumeField& src, ScalarSurfaceField& dst + ) const { - std::string nameInterpolated = "interpolated_" + volField.name; - ScalarSurfaceField surfaceField( - exec_, nameInterpolated, mesh_, createCalculatedBCs>(mesh_) - ); - interpolate(surfaceField, volField); - return surfaceField; + interpolationKernel_->interpolate(flux, src, dst); } void interpolate( - const ScalarSurfaceField& faceFlux, - const VolumeField& volField, - ScalarSurfaceField& surfaceField + const ScalarSurfaceField& flux, const VolumeField& src, VolumeField& dst ) const { - interpolationKernel_->interpolate(faceFlux, volField, surfaceField); + interpolationKernel_->interpolate(flux, src, dst); + } + + ScalarSurfaceField interpolate(const VolumeField& src) const + { + std::string nameInterpolated = "interpolated_" + volField.name; + ScalarSurfaceField dst( + exec_, nameInterpolated, mesh_, createCalculatedBCs>(mesh_) + ); + interpolate(dst, src); + return dst; } ScalarSurfaceField - interpolate(const ScalarSurfaceField& faceFlux, const VolumeField& volField) const + interpolate(const ScalarSurfaceField& flux, const VolumeField& src) const { - std::string name = "interpolated_" + volField.name; - ScalarSurfaceField surfaceField( + std::string name = "interpolated_" + src.name; + ScalarSurfaceField dst( exec_, name, mesh_, createCalculatedBCs>(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 d67603cd5..1ba6e2678 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp @@ -29,22 +29,16 @@ class Upwind : public SurfaceInterpolationFactory::Register static std::string schema() { return "none"; } - void interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const override; + void interpolate(const VolumeField& src, SurfaceField& dst) const override; + + void interpolate(const VolumeField& src, SurfaceField& dst) const override; void interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const override; - void interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const override; - void interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const override; std::unique_ptr clone() const override; diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index ed2236174..2097f7aae 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -9,44 +9,45 @@ 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 void computeLinearInterpolation( - const VolumeField& volField, - const std::shared_ptr geometryScheme, - SurfaceField& surfaceField + const VolumeField& src, + const SurfaceField& weights, + SurfaceField& 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 = src.internalField().span(); + const auto weightS = weights.internalField().span(); + const auto ownerS = dst.mesh().faceOwner().span(); + const auto neighS = dst.mesh().faceNeighbour().span(); + const auto boundS = src.boundaryField().value().span(); + 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(sOwner[facei]); - size_t nei = static_cast(sNeighbour[facei]); + size_t own = static_cast(ownerS[facei]); + size_t nei = static_cast(neighS[facei]); if (facei < nInternalFaces) { - std::cout << __FILE__ << " sWeight" << sWeight[facei] << "\n"; - std::cout << __FILE__ << " sVolfield[own]" << sVolField[own] << "\n"; - std::cout << __FILE__ << " sVolfield[nei]" << sVolField[nei] << "\n"; - sfield[facei] = - sWeight[facei] * sVolField[own] + (1 - sWeight[facei]) * sVolField[nei]; - std::cout << __FILE__ << " sfield[facei]" << sfield[facei] << "\n"; + 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]; } } ); @@ -60,34 +61,32 @@ Linear::Linear(const Executor& exec, const UnstructuredMesh& mesh) : SurfaceInterpolationFactory::Register(exec, mesh), geometryScheme_(GeometryScheme::readOrCreate(mesh)) {}; -void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const +void Linear::interpolate(const VolumeField& src, SurfaceField& dst) const { - computeLinearInterpolation(volField, geometryScheme_, surfaceField); + computeLinearInterpolation(src, geometryScheme_->weights(), dst); } -void Linear::interpolate( - [[maybe_unused]] const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField -) const +void Linear::interpolate(const VolumeField& src, SurfaceField& dst) const { - interpolate(volField, surfaceField); + computeLinearInterpolation(src, geometryScheme_->weights(), dst); } -void Linear::interpolate(const VolumeField& volField, SurfaceField& surfaceField) - const +void Linear::interpolate( + [[maybe_unused]] const SurfaceField& flux, + const VolumeField& src, + SurfaceField& dst +) const { - computeLinearInterpolation(volField, geometryScheme_, surfaceField); + interpolate(src, dst); } void Linear::interpolate( - [[maybe_unused]] const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + [[maybe_unused]] const SurfaceField& flux, + const VolumeField& src, + SurfaceField& dst ) const { - interpolate(volField, surfaceField); + interpolate(src, dst); } diff --git a/src/finiteVolume/cellCentred/interpolation/upwind.cpp b/src/finiteVolume/cellCentred/interpolation/upwind.cpp index d875fd28b..8c655ff6e 100644 --- a/src/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/src/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -9,47 +9,55 @@ 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 void computeUpwindInterpolation( - const SurfaceField& faceFlux, - const VolumeField& volField, - const std::shared_ptr geometryScheme, - SurfaceField& surfaceField + const VolumeField& src, + const SurfaceField& flux, + const SurfaceField& weights, + SurfaceField& 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 = src.internalField().span(); + const auto weightS = weights.internalField().span(); + const auto ownerS = dst.mesh().faceOwner().span(); + const auto neighS = dst.mesh().faceNeighbour().span(); + const auto boundS = src.boundaryField().value().span(); + const auto fluxS = flux.internalField().span(); + 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(sOwner[facei]); - sfield[facei] = sVolField[own]; + size_t own = static_cast(ownerS[facei]); + dstS[facei] = srcS[own]; } else { - size_t nei = static_cast(sNeighbour[facei]); - sfield[facei] = sVolField[nei]; + size_t nei = static_cast(neighS[facei]); + dstS[facei] = srcS[nei]; } } else { - sfield[facei] = sWeight[facei] * sBField[facei - nInternalFaces]; + dstS[facei] = weightS[facei] * boundS[facei - nInternalFaces]; } } ); @@ -60,39 +68,31 @@ Upwind::Upwind(const Executor& exec, const UnstructuredMesh& mesh, [[maybe_unuse geometryScheme_(GeometryScheme::readOrCreate(mesh)) {}; void Upwind::interpolate( - [[maybe_unused]] const VolumeField& volField, - [[maybe_unused]] SurfaceField& surfaceField + [[maybe_unused]] const VolumeField& src, [[maybe_unused]] SurfaceField& dst ) const { NF_ERROR_EXIT("limited scheme require a faceFlux"); } void Upwind::interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const { - // FIXME: - // computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField); + computeUpwindInterpolation(src, flux, geometryScheme_->weights(), dst); } void Upwind::interpolate( - [[maybe_unused]] const VolumeField& volField, - [[maybe_unused]] SurfaceField& surfaceField + [[maybe_unused]] const VolumeField& src, [[maybe_unused]] SurfaceField& dst ) const { NF_ERROR_EXIT("limited scheme require a faceFlux"); } void Upwind::interpolate( - const SurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const SurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const { - // FIXME: - // computeUpwindInterpolation(faceFlux, volField, geometryScheme_, surfaceField); + computeUpwindInterpolation(src, flux, geometryScheme_->weights(), dst); } std::unique_ptr Upwind::clone() const 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 6c1e482e1..40bac4766 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -46,6 +46,7 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) linear.interpolate(in, out); out.correctBoundaryConditions(); + // FIXME: needs to copy back to host for (int i = 0; i < out.internalField().size(); i++) { REQUIRE(out.internalField()[i] == NeoFOAM::one::value); diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp new file mode 100644 index 000000000..8a432f1ee --- /dev/null +++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -0,0 +1,56 @@ +// 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 +#include +#include +#include + +#include "NeoFOAM/NeoFOAM.hpp" + +using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; +using NeoFOAM::finiteVolume::cellCentred::VolumeField; +using NeoFOAM::finiteVolume::cellCentred::SurfaceField; +using NeoFOAM::Input; + +template +using I = std::initializer_list; + +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 = NeoFOAM::create1DUniformMesh(exec, 10); + Input input = NeoFOAM::TokenList({std::string("upwind")}); + auto upwind = SurfaceInterpolation(exec, mesh, input); + std::vector> bcs {}; + for (auto patchi : I {0, 1}) + { + NeoFOAM::Dictionary dict; + dict.insert("type", std::string("fixedValue")); + dict.insert("fixedValue", NeoFOAM::one::value); + bcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); + } + + auto in = VolumeField(exec, "in", mesh, {}); + auto flux = SurfaceField(exec, "flux", mesh, {}); + auto out = SurfaceField(exec, "out", mesh, bcs); + + fill(flux.internalField(), NeoFOAM::one::value); + fill(in.internalField(), NeoFOAM::one::value); + + upwind.interpolate(flux, in, out); + out.correctBoundaryConditions(); + + for (int i = 0; i < out.internalField().size(); i++) + { + REQUIRE(out.internalField()[i] == NeoFOAM::one::value); + } +} From cdcc44e285214ebbb9390d540e7840262ed08b89 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 19:29:35 +0100 Subject: [PATCH 12/35] fixup! compilation issue, more rename --- .../interpolation/surfaceInterpolation.hpp | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp index 01028121d..e4acbd2c4 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp @@ -50,22 +50,16 @@ class SurfaceInterpolationFactory : virtual ~SurfaceInterpolationFactory() {} // Virtual destructor - virtual void - interpolate(const VolumeField& volField, SurfaceField& surfaceField) const = 0; + virtual void interpolate(const VolumeField& src, SurfaceField& dst) const = 0; virtual void interpolate( - const ScalarSurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const ScalarSurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const = 0; - virtual void - interpolate(const VolumeField& volField, SurfaceField& surfaceField) const = 0; + virtual void interpolate(const VolumeField& src, SurfaceField& dst) const = 0; virtual void interpolate( - const ScalarSurfaceField& faceFlux, - const VolumeField& volField, - SurfaceField& surfaceField + const ScalarSurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const = 0; // Pure virtual function for cloning @@ -122,7 +116,7 @@ class SurfaceInterpolation } void interpolate( - const ScalarSurfaceField& flux, const VolumeField& src, VolumeField& dst + const ScalarSurfaceField& flux, const VolumeField& src, SurfaceField& dst ) const { interpolationKernel_->interpolate(flux, src, dst); @@ -130,7 +124,7 @@ class SurfaceInterpolation ScalarSurfaceField interpolate(const VolumeField& src) const { - std::string nameInterpolated = "interpolated_" + volField.name; + std::string nameInterpolated = "interpolated_" + src.name; ScalarSurfaceField dst( exec_, nameInterpolated, mesh_, createCalculatedBCs>(mesh_) ); From 7110aa35b6e0d2f6f37f48305994a3395cfeb636 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 9 Feb 2025 19:51:21 +0100 Subject: [PATCH 13/35] add upwind benchmark --- benchmarks/CMakeLists.txt | 1 + .../cellCentred/interpolation/CMakeLists.txt | 5 +++ .../cellCentred/interpolation/linear.cpp | 6 ++- .../cellCentred/interpolation/upwind.cpp | 45 +++++++++++++++++++ 4 files changed, 56 insertions(+), 1 deletion(-) create mode 100644 benchmarks/finiteVolume/cellCentred/interpolation/CMakeLists.txt create mode 100644 benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp 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 index c551d6355..1989eeb03 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -10,9 +10,11 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; +using NeoFOAM::Input; -TEST_CASE("DivOperator::div", "[bench]") +TEST_CASE("linear", "[bench]") { + using TestType = NeoFOAM::scalar; auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20); NeoFOAM::Executor exec = GENERATE( @@ -24,6 +26,8 @@ TEST_CASE("DivOperator::div", "[bench]") std::string execName = std::visit([](auto e) { return e.name(); }, exec); NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + Input input = NeoFOAM::TokenList({std::string("linear")}); + auto linear = SurfaceInterpolation(exec, mesh, input); auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, surfaceBCs); diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp new file mode 100644 index 000000000..631703a82 --- /dev/null +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -0,0 +1,45 @@ +// 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 NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; +using NeoFOAM::finiteVolume::cellCentred::VolumeField; +using NeoFOAM::finiteVolume::cellCentred::SurfaceField; +using NeoFOAM::Input; + +TEST_CASE("upwind", "[bench]") +{ + using TestType = NeoFOAM::scalar; + 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>(mesh); + Input input = NeoFOAM::TokenList({std::string("upwind")}); + auto upwind = SurfaceInterpolation(exec, mesh, input); + + auto in = VolumeField(exec, "in", mesh, {}); + auto flux = SurfaceField(exec, "flux", mesh, {}); + auto out = SurfaceField(exec, "out", mesh, surfaceBCs); + + fill(flux.internalField(), NeoFOAM::one::value); + fill(in.internalField(), NeoFOAM::one::value); + + // capture the value of size as section name + DYNAMIC_SECTION("" << size) + { + + BENCHMARK(std::string(execName)) { return (upwind.interpolate(flux, in, out)); }; + } +} From 990d3b49f931c5950ffbd436eb0962179f3e60ae Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 08:58:32 +0100 Subject: [PATCH 14/35] fix formating, fix copy to host --- benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp | 1 - benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp | 1 - test/finiteVolume/cellCentred/interpolation/linear.cpp | 4 ++-- test/finiteVolume/cellCentred/interpolation/upwind.cpp | 3 ++- 4 files changed, 4 insertions(+), 5 deletions(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp index 1989eeb03..7eca684cc 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -37,7 +37,6 @@ TEST_CASE("linear", "[bench]") // 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 index 631703a82..df4bf1b8e 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -39,7 +39,6 @@ TEST_CASE("upwind", "[bench]") // capture the value of size as section name DYNAMIC_SECTION("" << size) { - BENCHMARK(std::string(execName)) { return (upwind.interpolate(flux, in, out)); }; } } diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 40bac4766..38c75248b 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -46,9 +46,9 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) linear.interpolate(in, out); out.correctBoundaryConditions(); - // FIXME: needs to copy back to host + auto outHost = out.internalField().copyToHost(); for (int i = 0; i < out.internalField().size(); i++) { - REQUIRE(out.internalField()[i] == NeoFOAM::one::value); + REQUIRE(outHost[i] == NeoFOAM::one::value); } } diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp index 8a432f1ee..ba78ea3a2 100644 --- a/test/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -49,8 +49,9 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) upwind.interpolate(flux, in, out); out.correctBoundaryConditions(); + auto outHost = out.internalField().copyToHost(); for (int i = 0; i < out.internalField().size(); i++) { - REQUIRE(out.internalField()[i] == NeoFOAM::one::value); + REQUIRE(outHost[i] == NeoFOAM::one::value); } } From c36931d51a054fd8b061be783edc429e947c85fd Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 09:04:14 +0100 Subject: [PATCH 15/35] fix copyright year --- test/finiteVolume/cellCentred/interpolation/upwind.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp index ba78ea3a2..f44c0667f 100644 --- a/test/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2024 NeoFOAM authors +// SPDX-FileCopyrightText: 2025 NeoFOAM authors #define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create // a custom main From 7f00740b6f3de0383dd2e6aa28027109fea413b3 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 13:57:41 +0100 Subject: [PATCH 16/35] add templated benchmark case, extract type for json --- .../finiteVolume/cellCentred/interpolation/linear.cpp | 10 +++++----- .../finiteVolume/cellCentred/interpolation/upwind.cpp | 2 +- scripts/catch2json.py | 8 +++++++- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp index 7eca684cc..fbaa57be1 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -1,20 +1,20 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023 NeoFOAM authors +// 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 using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; using NeoFOAM::Input; -TEST_CASE("linear", "[bench]") +TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) { - using TestType = NeoFOAM::scalar; auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20); NeoFOAM::Executor exec = GENERATE( @@ -25,14 +25,14 @@ TEST_CASE("linear", "[bench]") std::string execName = std::visit([](auto e) { return e.name(); }, exec); NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); - auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); Input input = NeoFOAM::TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, surfaceBCs); - fill(in.internalField(), NeoFOAM::one::value); + fill(in.internalField(), NeoFOAM::one::value); // capture the value of size as section name DYNAMIC_SECTION("" << size) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp index df4bf1b8e..9e77a53a5 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023 NeoFOAM authors +// SPDX-FileCopyrightText: 2025 NeoFOAM authors #define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create // a custom main 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 From a79c5375856455745e3d17235ef90af2c256e8e4 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 14:49:02 +0100 Subject: [PATCH 17/35] template also upwdind benchmark --- .../finiteVolume/cellCentred/interpolation/upwind.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp index 9e77a53a5..c2e7a026f 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -7,12 +7,14 @@ #include "NeoFOAM/NeoFOAM.hpp" #include "../../../catch_main.hpp" +#include + using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; using NeoFOAM::Input; -TEST_CASE("upwind", "[bench]") +TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) { using TestType = NeoFOAM::scalar; auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20); @@ -25,7 +27,7 @@ TEST_CASE("upwind", "[bench]") std::string execName = std::visit([](auto e) { return e.name(); }, exec); NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); - auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); Input input = NeoFOAM::TokenList({std::string("upwind")}); auto upwind = SurfaceInterpolation(exec, mesh, input); @@ -34,7 +36,7 @@ TEST_CASE("upwind", "[bench]") auto out = SurfaceField(exec, "out", mesh, surfaceBCs); fill(flux.internalField(), NeoFOAM::one::value); - fill(in.internalField(), NeoFOAM::one::value); + fill(in.internalField(), NeoFOAM::one::value); // capture the value of size as section name DYNAMIC_SECTION("" << size) From 269a3fbb4de328d2cb3374dd214fc474e5bce9b2 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 14:49:29 +0100 Subject: [PATCH 18/35] slight format change --- benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp index fbaa57be1..19d1538b3 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -6,6 +6,7 @@ #include "NeoFOAM/NeoFOAM.hpp" #include "../../../catch_main.hpp" + #include using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; From 92b8e3419bc1b484b5bcaac4315eee0e1c67cf25 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 17:40:48 +0100 Subject: [PATCH 19/35] fixup! TestType without NeoFOAM:: --- benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp index c2e7a026f..3976c2384 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -36,7 +36,7 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) auto out = SurfaceField(exec, "out", mesh, surfaceBCs); fill(flux.internalField(), NeoFOAM::one::value); - fill(in.internalField(), NeoFOAM::one::value); + fill(in.internalField(), NeoFOAM::one::value); // capture the value of size as section name DYNAMIC_SECTION("" << size) From 059cb983422754d79f00555d665754deb6334d5b Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 17:52:03 +0100 Subject: [PATCH 20/35] use Namespace NeoFoam to make test/benchmark less verbose --- .../cellCentred/interpolation/linear.cpp | 12 ++++++---- .../cellCentred/interpolation/upwind.cpp | 16 +++++++++----- .../cellCentred/interpolation/linear.cpp | 16 +++++++++----- .../cellCentred/interpolation/upwind.cpp | 22 +++++++++++-------- 4 files changed, 41 insertions(+), 25 deletions(-) diff --git a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp index 19d1538b3..af5dfd859 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/linear.cpp @@ -12,7 +12,9 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; -using NeoFOAM::Input; + +namespace NeoFOAM +{ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) { @@ -25,15 +27,15 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); + UnstructuredMesh mesh = create1DUniformMesh(exec, size); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); - Input input = NeoFOAM::TokenList({std::string("linear")}); + Input input = TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, surfaceBCs); - fill(in.internalField(), NeoFOAM::one::value); + fill(in.internalField(), one::value); // capture the value of size as section name DYNAMIC_SECTION("" << size) @@ -41,3 +43,5 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) 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 index 3976c2384..1e3b53dfe 100644 --- a/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/benchmarks/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -14,9 +14,11 @@ using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; using NeoFOAM::Input; +namespace NeoFOAM +{ + TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) { - using TestType = NeoFOAM::scalar; auto size = GENERATE(1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20); NeoFOAM::Executor exec = GENERATE( @@ -26,17 +28,17 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); + UnstructuredMesh mesh = create1DUniformMesh(exec, size); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); - Input input = NeoFOAM::TokenList({std::string("upwind")}); + Input input = TokenList({std::string("upwind")}); auto upwind = SurfaceInterpolation(exec, mesh, input); auto in = VolumeField(exec, "in", mesh, {}); - auto flux = SurfaceField(exec, "flux", mesh, {}); + auto flux = SurfaceField(exec, "flux", mesh, {}); auto out = SurfaceField(exec, "out", mesh, surfaceBCs); - fill(flux.internalField(), NeoFOAM::one::value); - fill(in.internalField(), NeoFOAM::one::value); + fill(flux.internalField(), one::value); + fill(in.internalField(), one::value); // capture the value of size as section name DYNAMIC_SECTION("" << size) @@ -44,3 +46,5 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) BENCHMARK(std::string(execName)) { return (upwind.interpolate(flux, in, out)); }; } } + +} diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 38c75248b..1ca843428 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -14,6 +14,9 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; +namespace NeoFOAM +{ + template using I = std::initializer_list; @@ -26,22 +29,22 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - auto mesh = NeoFOAM::create1DUniformMesh(exec, 10); - 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> bcs {}; for (auto patchi : I {0, 1}) { - NeoFOAM::Dictionary dict; + Dictionary dict; dict.insert("type", std::string("fixedValue")); - dict.insert("fixedValue", NeoFOAM::one::value); + dict.insert("fixedValue", one::value); bcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); } auto in = VolumeField(exec, "in", mesh, {}); auto out = SurfaceField(exec, "out", mesh, bcs); - fill(in.internalField(), NeoFOAM::one::value); + fill(in.internalField(), one::value); linear.interpolate(in, out); out.correctBoundaryConditions(); @@ -49,6 +52,7 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) auto outHost = out.internalField().copyToHost(); for (int i = 0; i < out.internalField().size(); i++) { - REQUIRE(outHost[i] == NeoFOAM::one::value); + REQUIRE(outHost[i] == one::value); } } +} diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp index f44c0667f..443ec8558 100644 --- a/test/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -13,7 +13,9 @@ using NeoFOAM::finiteVolume::cellCentred::SurfaceInterpolation; using NeoFOAM::finiteVolume::cellCentred::VolumeField; using NeoFOAM::finiteVolume::cellCentred::SurfaceField; -using NeoFOAM::Input; + +namespace NeoFOAM +{ template using I = std::initializer_list; @@ -27,24 +29,24 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - auto mesh = NeoFOAM::create1DUniformMesh(exec, 10); - Input input = NeoFOAM::TokenList({std::string("upwind")}); + auto mesh = create1DUniformMesh(exec, 10); + Input input = TokenList({std::string("upwind")}); auto upwind = SurfaceInterpolation(exec, mesh, input); std::vector> bcs {}; for (auto patchi : I {0, 1}) { - NeoFOAM::Dictionary dict; + Dictionary dict; dict.insert("type", std::string("fixedValue")); - dict.insert("fixedValue", NeoFOAM::one::value); + dict.insert("fixedValue", one::value); bcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); } auto in = VolumeField(exec, "in", mesh, {}); - auto flux = SurfaceField(exec, "flux", mesh, {}); + auto flux = SurfaceField(exec, "flux", mesh, {}); auto out = SurfaceField(exec, "out", mesh, bcs); - fill(flux.internalField(), NeoFOAM::one::value); - fill(in.internalField(), NeoFOAM::one::value); + fill(flux.internalField(), one::value); + fill(in.internalField(), one::value); upwind.interpolate(flux, in, out); out.correctBoundaryConditions(); @@ -52,6 +54,8 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) auto outHost = out.internalField().copyToHost(); for (int i = 0; i < out.internalField().size(); i++) { - REQUIRE(outHost[i] == NeoFOAM::one::value); + REQUIRE(outHost[i] == one::value); } } + +} From 391f7584c8de6f2d4646c8ed185b7297289193cb Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 15 Feb 2025 15:25:50 +0100 Subject: [PATCH 21/35] use spans() to unpack --- .../cellCentred/interpolation/linear.cpp | 12 +++++++----- .../cellCentred/interpolation/upwind.cpp | 14 ++++++++------ 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index 2097f7aae..871c92c44 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -28,11 +28,13 @@ void computeLinearInterpolation( { const auto exec = dst.exec(); auto dstS = dst.internalField().span(); - const auto srcS = src.internalField().span(); - const auto weightS = weights.internalField().span(); - const auto ownerS = dst.mesh().faceOwner().span(); - const auto neighS = dst.mesh().faceNeighbour().span(); - const auto boundS = src.boundaryField().value().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( diff --git a/src/finiteVolume/cellCentred/interpolation/upwind.cpp b/src/finiteVolume/cellCentred/interpolation/upwind.cpp index 8c655ff6e..85e7e1de2 100644 --- a/src/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/src/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -30,12 +30,14 @@ void computeUpwindInterpolation( { const auto exec = dst.exec(); auto dstS = dst.internalField().span(); - const auto srcS = src.internalField().span(); - const auto weightS = weights.internalField().span(); - const auto ownerS = dst.mesh().faceOwner().span(); - const auto neighS = dst.mesh().faceNeighbour().span(); - const auto boundS = src.boundaryField().value().span(); - const auto fluxS = flux.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( From aab61cd889e6661034791232a9a1186798a08ffd Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 15 Feb 2025 15:28:46 +0100 Subject: [PATCH 22/35] get neig/own only for internal faces --- src/finiteVolume/cellCentred/interpolation/linear.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index 871c92c44..37c7dddd0 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -41,10 +41,10 @@ void computeLinearInterpolation( exec, {0, dstS.size()}, KOKKOS_LAMBDA(const size_t facei) { - size_t own = static_cast(ownerS[facei]); - size_t nei = static_cast(neighS[facei]); if (facei < nInternalFaces) { + size_t own = static_cast(ownerS[facei]); + size_t nei = static_cast(neighS[facei]); dstS[facei] = weightS[facei] * srcS[own] + (1 - weightS[facei]) * srcS[nei]; } else From a771bfc3badd33f35fe426796d2840c8a68bc634 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Mon, 27 Jan 2025 09:28:20 +0100 Subject: [PATCH 23/35] wip --- .../cellCentred/operator/gaussGreenGrad.cpp | 41 ++++++++ .../cellCentred/operator/gradOperator.cpp | 42 ++++++++ .../cellCentred/operators/divOperator.hpp | 2 +- .../cellCentred/operators/gaussGreenDiv.cpp | 95 +++++++++++-------- 4 files changed, 137 insertions(+), 43 deletions(-) create mode 100644 benchmarks/finiteVolume/cellCentred/operator/gaussGreenGrad.cpp create mode 100644 benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp 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>(mesh); + fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); + NeoFOAM::fill(faceFlux.internalField(), 1.0); + + auto volumeBCs = fvcc::createCalculatedBCs>(mesh); + fvcc::VolumeField phi(exec, "vf", mesh, volumeBCs); + fvcc::VolumeField 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/benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp b/benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp new file mode 100644 index 000000000..78b3bc797 --- /dev/null +++ b/benchmarks/finiteVolume/cellCentred/operator/gradOperator.cpp @@ -0,0 +1,42 @@ +// 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::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); + NeoFOAM::UnstructuredMesh mesh = NeoFOAM::create1DUniformMesh(exec, size); + auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); + NeoFOAM::fill(faceFlux.internalField(), 1.0); + + auto volumeBCs = fvcc::createCalculatedBCs>(mesh); + fvcc::VolumeField phi(exec, "vf", mesh, volumeBCs); + fvcc::VolumeField 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::GradOperator(Operator::Type::Explicit, faceFlux, phi, input); + + BENCHMARK(std::string(execName)) { return (op.grad(divPhi)); }; + } +} diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index 2ed0fb16f..e3f8bc448 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 diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 890bb3779..72c057d5b 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -8,51 +8,46 @@ namespace NeoFOAM::finiteVolume::cellCentred { +/* @brief free standing function implementation of the divergence operator +** ie computes \sum_f S_f \cdot \phi_f +** +** @param faceFlux +*/ +template void computeDiv( - const SurfaceField& faceFlux, - const VolumeField& phi, - const SurfaceInterpolation& surfInterp, - Field& divPhi + const Executor& exec, + size_t nInternalFaces, + size_t nBoundaryFaces, + const std::span neighbour, + const std::span owner, + const std::span faceCells, + const std::span faceFlux, + const std::span phiF, + const std::span V, + std::span divPhi ) { - const UnstructuredMesh& mesh = phi.mesh(); - const auto exec = phi.exec(); - SurfaceField phif( - exec, "phif", mesh, createCalculatedBCs>(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(exec)) { - for (size_t i = 0; i < nInternalFaces; i++) + for (size_t i = 0; i < faceFlux.size(); i++) { - scalar flux = surfFaceFlux[i] * surfPhif[i]; - surfDivPhi[static_cast(surfOwner[i])] += flux; - surfDivPhi[static_cast(surfNeighbour[i])] -= flux; + scalar flux = faceFlux[i] * phiF[i]; + divPhi[static_cast(owner[i])] += flux; + divPhi[static_cast(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(surfFaceCells[i - nInternalFaces]); - scalar valueOwn = surfFaceFlux[i] * surfPhif[i]; - surfDivPhi[own] += valueOwn; + auto own = static_cast(faceCells[i - nInternalFaces]); + scalar valueOwn = faceFlux[i] * phiF[i]; + divPhi[own] += valueOwn; } - for (size_t celli = 0; celli < mesh.nCells(); celli++) + for (size_t celli = 0; celli < nCells; celli++) { - surfDivPhi[celli] *= 1 / surfV[celli]; + divPhi[celli] *= 1 / V[celli]; } } else @@ -61,30 +56,29 @@ void computeDiv( exec, {0, nInternalFaces}, KOKKOS_LAMBDA(const size_t i) { - scalar flux = surfFaceFlux[i] * surfPhif[i]; - Kokkos::atomic_add(&surfDivPhi[static_cast(surfOwner[i])], flux); - Kokkos::atomic_sub(&surfDivPhi[static_cast(surfNeighbour[i])], flux); + scalar flux = faceFlux[i] * phiF[i]; + Kokkos::atomic_add(&divPhi[static_cast(owner[i])], flux); + Kokkos::atomic_sub(&divPhi[static_cast(neighbour[i])], flux); } ); parallelFor( exec, - {nInternalFaces, surfPhif.size()}, + {nInternalFaces, nBoundaryFaces}, KOKKOS_LAMBDA(const size_t i) { - auto own = static_cast(surfFaceCells[i - nInternalFaces]); - scalar valueOwn = surfFaceFlux[i] * surfPhif[i]; - Kokkos::atomic_add(&surfDivPhi[own], valueOwn); + auto own = static_cast(faceCells[i - nInternalFaces]); + scalar valueOwn = faceFlux[i] * phiF[i]; + Kokkos::atomic_add(&divPhi[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) { divPhi[celli] *= 1 / V[celli]; } ); } } + void computeDiv( const SurfaceField& faceFlux, const VolumeField& phi, @@ -92,7 +86,24 @@ void computeDiv( VolumeField& divPhi ) { + const UnstructuredMesh& mesh = phi.mesh(); + const auto exec = phi.exec(); + SurfaceField phif( + exec, "phif", mesh, createCalculatedBCs>(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(); + const auto surfV = mesh.cellVolumes().span(); Field& divPhiField = divPhi.internalField(); + size_t nInternalFaces = mesh.nInternalFaces(); computeDiv(faceFlux, phi, surfInterp, divPhiField); } From e70cd488d8cdc308f48d29719e2ff4b8fad7d622 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Mon, 3 Feb 2025 19:30:08 +0100 Subject: [PATCH 24/35] wip make it compile --- .../cellCentred/operators/gaussGreenDiv.cpp | 49 +++++++++---------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 72c057d5b..44115b676 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -18,12 +18,12 @@ void computeDiv( const Executor& exec, size_t nInternalFaces, size_t nBoundaryFaces, - const std::span neighbour, - const std::span owner, - const std::span faceCells, - const std::span faceFlux, - const std::span phiF, - const std::span V, + std::span neighbour, + std::span owner, + std::span faceCells, + std::span faceFlux, + std::span phiF, + std::span V, std::span divPhi ) { @@ -79,11 +79,12 @@ void computeDiv( } +template void computeDiv( const SurfaceField& faceFlux, const VolumeField& phi, const SurfaceInterpolation& surfInterp, - VolumeField& divPhi + VolumeField& divPhi ) { const UnstructuredMesh& mesh = phi.mesh(); @@ -92,19 +93,22 @@ void computeDiv( exec, "phif", mesh, createCalculatedBCs>(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(); - const auto surfV = mesh.cellVolumes().span(); - Field& divPhiField = divPhi.internalField(); size_t nInternalFaces = mesh.nInternalFaces(); - computeDiv(faceFlux, phi, surfInterp, divPhiField); + size_t nBoundaryFaces = mesh.nBoundaryFaces(); + computeDiv( + exec, + nInternalFaces, + nBoundaryFaces, + mesh.faceNeighbour().span(), + mesh.faceOwner().span(), + mesh.boundaryMesh().faceCells().span(), + faceFlux.internalField().span(), + phif.internalField().span(), + mesh.cellVolumes().span(), + divPhi.internalField().span() + ); } GaussGreenDiv::GaussGreenDiv( @@ -117,14 +121,7 @@ void GaussGreenDiv::div( VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) { - computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); -}; - -void GaussGreenDiv::div( - Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi -) -{ - computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); }; VolumeField @@ -134,7 +131,7 @@ GaussGreenDiv::div(const SurfaceField& faceFlux, VolumeField& ph VolumeField divPhi( exec_, name, mesh_, createCalculatedBCs>(mesh_) ); - computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); return divPhi; }; From 04105586f1960848f9b17f25e129501e696ecbe5 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Mon, 3 Feb 2025 20:45:58 +0100 Subject: [PATCH 25/35] start parameterise test --- .../cellCentred/operator/divOperator.cpp | 40 ++++++++++--------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index 241acd149..2ef5e0e7c 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include "NeoFOAM/NeoFOAM.hpp" @@ -13,7 +14,8 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; using Operator = NeoFOAM::dsl::Operator; -TEST_CASE("DivOperator") + +TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector) { NeoFOAM::Executor exec = GENERATE( NeoFOAM::Executor(NeoFOAM::SerialExecutor {}), @@ -29,22 +31,22 @@ TEST_CASE("DivOperator") fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); NeoFOAM::fill(faceFlux.internalField(), 1.0); - auto volumeBCs = fvcc::createCalculatedBCs>(mesh); - fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); - NeoFOAM::fill(phi.internalField(), 1.0); - - SECTION("Construct from Token" + execName) - { - NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - fvcc::DivOperator(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); - } + auto volumeBCs = fvcc::createCalculatedBCs>(mesh); + fvcc::VolumeField vecPhi(exec, "sf", mesh, volumeBCs); + // NeoFOAM::fill(vecPhi.internalField(), NeoFOAM::Vector{1.0, 1.0, 1.0}); + + // SECTION("Construct from Token" + execName) + // { + // NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); + // fvcc::DivOperator(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); + // } } From 518bca36270610146fe404dea00ddaed620499c0 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 4 Feb 2025 20:21:05 +0100 Subject: [PATCH 26/35] wip --- .../cellCentred/operators/divOperator.hpp | 4 ++-- .../cellCentred/operators/gaussGreenDiv.cpp | 13 ++++++++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index e3f8bc448..fc8268269 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -77,7 +77,7 @@ class DivOperator : public dsl::OperatorMixin> DivOperator( dsl::Operator::Type termType, const SurfaceField& faceFlux, - VolumeField& phi, + VolumeField& phi, Input input ) : dsl::OperatorMixin>(phi.exec(), phi, termType), faceFlux_(faceFlux), @@ -86,7 +86,7 @@ class DivOperator : public dsl::OperatorMixin> DivOperator( dsl::Operator::Type termType, const SurfaceField& faceFlux, - VolumeField& phi, + VolumeField& phi, std::unique_ptr divOperatorStrategy ) : dsl::OperatorMixin>(phi.exec(), phi, termType), faceFlux_(faceFlux), diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 44115b676..2f4254263 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -84,7 +84,7 @@ void computeDiv( const SurfaceField& faceFlux, const VolumeField& phi, const SurfaceInterpolation& surfInterp, - VolumeField& divPhi + Field& divPhi ) { const UnstructuredMesh& mesh = phi.mesh(); @@ -107,7 +107,7 @@ void computeDiv( faceFlux.internalField().span(), phif.internalField().span(), mesh.cellVolumes().span(), - divPhi.internalField().span() + divPhi.span() ); } @@ -120,6 +120,13 @@ GaussGreenDiv::GaussGreenDiv( void GaussGreenDiv::div( VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) +{ + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi.internalField()); +}; + +void GaussGreenDiv::div( + Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi +) { computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); }; @@ -131,7 +138,7 @@ GaussGreenDiv::div(const SurfaceField& faceFlux, VolumeField& ph VolumeField divPhi( exec_, name, mesh_, createCalculatedBCs>(mesh_) ); - computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi.internalField()); return divPhi; }; From 0f5714aea52967e9c8f76a0ea045df0d380272fb Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Fri, 7 Feb 2025 21:41:56 +0100 Subject: [PATCH 27/35] wip make it compile --- include/NeoFOAM/dsl/explicit.hpp | 5 +- .../cellCentred/operators/divOperator.hpp | 40 +++++++++++----- .../cellCentred/operators/gaussGreenDiv.hpp | 17 +++++-- .../cellCentred/operators/gaussGreenDiv.cpp | 48 +++++++++++++++---- 4 files changed, 83 insertions(+), 27 deletions(-) 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& faceFlux, fvcc::VolumeField& phi) +Operator div(const fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return Operator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi)); + return Operator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi)); } diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index fc8268269..1c8418370 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -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& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) = 0; @@ -48,9 +51,19 @@ class DivOperatorFactory : virtual void div(Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi) = 0; + virtual void + div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi + ) = 0; + + virtual void + div(Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi) = 0; + virtual VolumeField div(const SurfaceField& faceFlux, VolumeField& phi) = 0; + virtual VolumeField + div(const SurfaceField& faceFlux, VolumeField& phi) = 0; + // Pure virtual function for cloning virtual std::unique_ptr clone() const = 0; @@ -61,14 +74,15 @@ class DivOperatorFactory : const UnstructuredMesh& mesh_; }; -class DivOperator : public dsl::OperatorMixin> +template +class DivOperator : public dsl::OperatorMixin> { public: // copy constructor DivOperator(const DivOperator& divOp) - : dsl::OperatorMixin>(divOp.exec_, divOp.field_, divOp.type_), + : dsl::OperatorMixin>(divOp.exec_, divOp.field_, divOp.type_), faceFlux_(divOp.faceFlux_), divOperatorStrategy_( divOp.divOperatorStrategy_ ? divOp.divOperatorStrategy_->clone() : nullptr @@ -80,8 +94,9 @@ class DivOperator : public dsl::OperatorMixin> VolumeField& phi, Input input ) - : dsl::OperatorMixin>(phi.exec(), phi, termType), faceFlux_(faceFlux), - divOperatorStrategy_(DivOperatorFactory::create(exec_, phi.mesh(), input)) {}; + : dsl::OperatorMixin>(phi.exec(), phi, termType), + faceFlux_(faceFlux), + divOperatorStrategy_(DivOperatorFactory::create(phi.exec(), phi.mesh(), input)) {}; DivOperator( dsl::Operator::Type termType, @@ -106,32 +121,35 @@ class DivOperator : public dsl::OperatorMixin> NF_ERROR_EXIT("DivOperatorStrategy not initialized"); } NeoFOAM::Field tmpsource(source.exec(), source.size(), 0.0); - divOperatorStrategy_->div(tmpsource, faceFlux_, field_); + divOperatorStrategy_->div(tmpsource, faceFlux_, this->getField()); source += tmpsource; } - void div(Field& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, getField()); } + void div(Field& divPhi) + { + divOperatorStrategy_->div(divPhi, faceFlux_, this->getField()); + } void div(VolumeField& 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(input)) { auto dict = std::get(input); - std::string schemeName = "div(" + faceFlux_.name + "," + field_.name + ")"; + std::string schemeName = "div(" + faceFlux_.name + "," + this->getField().name + ")"; auto tokens = dict.subDict("divSchemes").get(schemeName); - divOperatorStrategy_ = DivOperatorFactory::create(exec(), mesh, tokens); + divOperatorStrategy_ = DivOperatorFactory::create(this->exec(), mesh, tokens); } else { auto tokens = std::get(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(const Executor& exec, const UnstructuredMesh& mesh, const Input& inputs); - void + virtual void div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) override; - void + virtual void div(Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) override; - VolumeField + virtual VolumeField div(const SurfaceField& faceFlux, VolumeField& phi) override; + virtual void + div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi + ) override; + + virtual void + div(Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi + ) override; + + virtual VolumeField + div(const SurfaceField& faceFlux, VolumeField& phi) override; + std::unique_ptr clone() const override; private: diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 2f4254263..fdb139259 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -21,7 +21,7 @@ void computeDiv( std::span neighbour, std::span owner, std::span faceCells, - std::span faceFlux, + std::span faceFlux, std::span phiF, std::span V, std::span divPhi @@ -33,7 +33,7 @@ void computeDiv( { for (size_t i = 0; i < faceFlux.size(); i++) { - scalar flux = faceFlux[i] * phiF[i]; + ValueType flux = faceFlux[i] * phiF[i]; divPhi[static_cast(owner[i])] += flux; divPhi[static_cast(neighbour[i])] -= flux; } @@ -41,7 +41,7 @@ void computeDiv( for (size_t i = nInternalFaces; i < nInternalFaces + nBoundaryFaces; i++) { auto own = static_cast(faceCells[i - nInternalFaces]); - scalar valueOwn = faceFlux[i] * phiF[i]; + ValueType valueOwn = faceFlux[i] * phiF[i]; divPhi[own] += valueOwn; } @@ -56,7 +56,7 @@ void computeDiv( exec, {0, nInternalFaces}, KOKKOS_LAMBDA(const size_t i) { - scalar flux = faceFlux[i] * phiF[i]; + ValueType flux = faceFlux[i] * phiF[i]; Kokkos::atomic_add(&divPhi[static_cast(owner[i])], flux); Kokkos::atomic_sub(&divPhi[static_cast(neighbour[i])], flux); } @@ -67,7 +67,7 @@ void computeDiv( {nInternalFaces, nBoundaryFaces}, KOKKOS_LAMBDA(const size_t i) { auto own = static_cast(faceCells[i - nInternalFaces]); - scalar valueOwn = faceFlux[i] * phiF[i]; + ValueType valueOwn = faceFlux[i] * phiF[i]; Kokkos::atomic_add(&divPhi[own], valueOwn); } ); @@ -82,18 +82,20 @@ void computeDiv( template void computeDiv( const SurfaceField& faceFlux, - const VolumeField& phi, + const VolumeField& phi, const SurfaceInterpolation& surfInterp, Field& divPhi ) { const UnstructuredMesh& mesh = phi.mesh(); const auto exec = phi.exec(); - SurfaceField phif( - exec, "phif", mesh, createCalculatedBCs>(mesh) + SurfaceField phif( + exec, "phif", mesh, createCalculatedBCs>(mesh) ); - fill(phif.internalField(), 0.0); - surfInterp.interpolate(faceFlux, phi, phif); + // FIXME: needs a ValueType::zero + // fill(phif.internalField(), 0.0); + // FIXME: needs a interpolate for + // surfInterp.interpolate(faceFlux, phi, phif); size_t nInternalFaces = mesh.nInternalFaces(); size_t nBoundaryFaces = mesh.nBoundaryFaces(); @@ -142,6 +144,32 @@ GaussGreenDiv::div(const SurfaceField& faceFlux, VolumeField& ph return divPhi; }; + +void GaussGreenDiv::div( + VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi +) +{ + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi.internalField()); +}; + +void GaussGreenDiv::div( + Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi +) +{ + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); +}; + +VolumeField +GaussGreenDiv::div(const SurfaceField& faceFlux, VolumeField& phi) +{ + std::string name = "div(" + faceFlux.name + "," + phi.name + ")"; + VolumeField divPhi( + exec_, name, mesh_, createCalculatedBCs>(mesh_) + ); + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi.internalField()); + return divPhi; +}; + std::unique_ptr GaussGreenDiv::clone() const { return std::make_unique(*this); From 8859c47b1026a71070409a7c0061bccdd6d88bcd Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 8 Feb 2025 09:03:31 +0100 Subject: [PATCH 28/35] uncomment build div operator test --- .../cellCentred/operator/divOperator.cpp | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index 2ef5e0e7c..1edafb518 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -32,21 +32,21 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector NeoFOAM::fill(faceFlux.internalField(), 1.0); auto volumeBCs = fvcc::createCalculatedBCs>(mesh); - fvcc::VolumeField vecPhi(exec, "sf", mesh, volumeBCs); + fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); // NeoFOAM::fill(vecPhi.internalField(), NeoFOAM::Vector{1.0, 1.0, 1.0}); - // SECTION("Construct from Token" + execName) - // { - // NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - // fvcc::DivOperator(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); - // } + SECTION("Construct from Token" + execName) + { + NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); + fvcc::DivOperator(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); + } } From ebdf83c908c0e89529bef135bf3d512a3f54b084 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Tue, 11 Feb 2025 20:28:04 +0100 Subject: [PATCH 29/35] use interpolation and fill functions --- src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp | 6 ++---- test/finiteVolume/cellCentred/operator/divOperator.cpp | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index fdb139259..6fb67cd08 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -92,10 +92,8 @@ void computeDiv( SurfaceField phif( exec, "phif", mesh, createCalculatedBCs>(mesh) ); - // FIXME: needs a ValueType::zero - // fill(phif.internalField(), 0.0); - // FIXME: needs a interpolate for - // surfInterp.interpolate(faceFlux, phi, phif); + fill(phif.internalField(), NeoFOAM::zero::value); + surfInterp.interpolate(faceFlux, phi, phif); size_t nInternalFaces = mesh.nInternalFaces(); size_t nBoundaryFaces = mesh.nBoundaryFaces(); diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index 1edafb518..f7bd776f8 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -33,7 +33,7 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector auto volumeBCs = fvcc::createCalculatedBCs>(mesh); fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); - // NeoFOAM::fill(vecPhi.internalField(), NeoFOAM::Vector{1.0, 1.0, 1.0}); + NeoFOAM::fill(phi.internalField(), NeoFOAM::one::value); SECTION("Construct from Token" + execName) { From e97a5a0b5c4e4cbc07872ab52204b705460f16b6 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Wed, 12 Feb 2025 20:56:28 +0100 Subject: [PATCH 30/35] fix failing test --- .../cellCentred/operators/divOperator.hpp | 4 ++-- .../finiteVolume/cellCentred/operator/divOperator.cpp | 11 +++++++++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index 1c8418370..f88f53a34 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -125,12 +125,12 @@ class DivOperator : public dsl::OperatorMixin> source += tmpsource; } - void div(Field& divPhi) + void div(Field& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, this->getField()); } - void div(VolumeField& divPhi) + void div(VolumeField& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, this->getField()); } diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index f7bd776f8..a2ea7505c 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -14,6 +14,8 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; using Operator = NeoFOAM::dsl::Operator; +namespace NeoFOAM +{ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector) { @@ -25,7 +27,7 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector std::string execName = std::visit([](auto e) { return e.name(); }, exec); // TODO take 1d mesh - NeoFOAM::UnstructuredMesh mesh = NeoFOAM::createSingleCellMesh(exec); + auto mesh = create1DUniformMesh(exec, 10); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); @@ -35,6 +37,8 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); NeoFOAM::fill(phi.internalField(), NeoFOAM::one::value); + auto result = NeoFOAM::Field(exec, phi.size()); + SECTION("Construct from Token" + execName) { NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); @@ -47,6 +51,9 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector {{std::string("DivOperator"), std::string("Gauss")}, {std::string("surfaceInterpolation"), std::string("linear")}} ); - fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); + auto op = fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); + op.div(result); } } + +} From dc9315c24688928ea62251c9c4dd3a9750d674c9 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 15 Feb 2025 19:31:22 +0100 Subject: [PATCH 31/35] improve interpolation tests --- .../cellCentred/interpolation/linear.cpp | 20 ++++++++++++++----- .../cellCentred/interpolation/upwind.cpp | 9 ++++++++- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/test/finiteVolume/cellCentred/interpolation/linear.cpp b/test/finiteVolume/cellCentred/interpolation/linear.cpp index 1ca843428..449c38cdd 100644 --- a/test/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/test/finiteVolume/cellCentred/interpolation/linear.cpp @@ -32,25 +32,35 @@ TEMPLATE_TEST_CASE("linear", "", NeoFOAM::scalar, NeoFOAM::Vector) auto mesh = create1DUniformMesh(exec, 10); Input input = TokenList({std::string("linear")}); auto linear = SurfaceInterpolation(exec, mesh, input); - std::vector> bcs {}; + std::vector> vbcs {}; + std::vector> sbcs {}; for (auto patchi : I {0, 1}) { Dictionary dict; dict.insert("type", std::string("fixedValue")); dict.insert("fixedValue", one::value); - bcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); + sbcs.push_back(fvcc::SurfaceBoundary(mesh, dict, patchi)); + vbcs.push_back(fvcc::VolumeBoundary(mesh, dict, patchi)); } - auto in = VolumeField(exec, "in", mesh, {}); - auto out = SurfaceField(exec, "out", mesh, bcs); + auto in = VolumeField(exec, "in", mesh, vbcs); + auto out = SurfaceField(exec, "out", mesh, sbcs); fill(in.internalField(), one::value); + in.correctBoundaryConditions(); linear.interpolate(in, out); out.correctBoundaryConditions(); auto outHost = out.internalField().copyToHost(); - for (int i = 0; i < out.internalField().size(); i++) + auto nInternal = mesh.nInternalFaces(); + auto nBoundary = mesh.nBoundaryFaces(); + for (int i = 0; i < nInternal; i++) + { + REQUIRE(outHost[i] == one::value); + } + + for (int i = nInternal; i < nInternal + nBoundary; i++) { REQUIRE(outHost[i] == one::value); } diff --git a/test/finiteVolume/cellCentred/interpolation/upwind.cpp b/test/finiteVolume/cellCentred/interpolation/upwind.cpp index 443ec8558..69f039382 100644 --- a/test/finiteVolume/cellCentred/interpolation/upwind.cpp +++ b/test/finiteVolume/cellCentred/interpolation/upwind.cpp @@ -52,7 +52,14 @@ TEMPLATE_TEST_CASE("upwind", "", NeoFOAM::scalar, NeoFOAM::Vector) out.correctBoundaryConditions(); auto outHost = out.internalField().copyToHost(); - for (int i = 0; i < out.internalField().size(); i++) + auto nInternal = mesh.nInternalFaces(); + auto nBoundary = mesh.nBoundaryFaces(); + for (int i = 0; i < nInternal; i++) + { + REQUIRE(outHost[i] == one::value); + } + + for (int i = nInternal; i < nInternal + nBoundary; i++) { REQUIRE(outHost[i] == one::value); } From a47f1f3a567dab4be05136318d499f7f22d99b3a Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sat, 15 Feb 2025 20:07:19 +0100 Subject: [PATCH 32/35] WIP add divOperator test values --- test/finiteVolume/cellCentred/operator/divOperator.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index a2ea7505c..31016cc4d 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -36,8 +36,10 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector auto volumeBCs = fvcc::createCalculatedBCs>(mesh); fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); NeoFOAM::fill(phi.internalField(), NeoFOAM::one::value); + phi.correctBoundaryConditions(); auto result = NeoFOAM::Field(exec, phi.size()); + NeoFOAM::fill(result, NeoFOAM::zero::value); SECTION("Construct from Token" + execName) { @@ -53,6 +55,14 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector ); 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::value); + } } } From 23fac66aa002f99716f39fb60d6e468067526a69 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 16 Feb 2025 09:47:34 +0100 Subject: [PATCH 33/35] fixup! divOperator test for CPU --- .../cellCentred/operators/gaussGreenDiv.cpp | 36 ++++++++++++------- .../cellCentred/operator/divOperator.cpp | 7 +++- 2 files changed, 29 insertions(+), 14 deletions(-) diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 6fb67cd08..9ba87ffd0 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -9,9 +9,15 @@ namespace NeoFOAM::finiteVolume::cellCentred /* @brief free standing function implementation of the divergence operator -** ie computes \sum_f S_f \cdot \phi_f +** 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 void computeDiv( @@ -24,30 +30,31 @@ void computeDiv( std::span faceFlux, std::span phiF, std::span V, - std::span divPhi + std::span res ) { size_t nCells {V.size()}; // check if the executor is GPU if (std::holds_alternative(exec)) { - for (size_t i = 0; i < faceFlux.size(); i++) + for (size_t i = 0; i < nInternalFaces; i++) { ValueType flux = faceFlux[i] * phiF[i]; - divPhi[static_cast(owner[i])] += flux; - divPhi[static_cast(neighbour[i])] -= flux; + res[static_cast(owner[i])] += flux; + res[static_cast(neighbour[i])] -= flux; } for (size_t i = nInternalFaces; i < nInternalFaces + nBoundaryFaces; i++) { auto own = static_cast(faceCells[i - nInternalFaces]); ValueType valueOwn = faceFlux[i] * phiF[i]; - divPhi[own] += valueOwn; + res[own] += valueOwn; } + // TODO does it make sense to store invVol and multiply? for (size_t celli = 0; celli < nCells; celli++) { - divPhi[celli] *= 1 / V[celli]; + res[celli] *= 1 / V[celli]; } } else @@ -57,23 +64,23 @@ void computeDiv( {0, nInternalFaces}, KOKKOS_LAMBDA(const size_t i) { ValueType flux = faceFlux[i] * phiF[i]; - Kokkos::atomic_add(&divPhi[static_cast(owner[i])], flux); - Kokkos::atomic_sub(&divPhi[static_cast(neighbour[i])], flux); + Kokkos::atomic_add(&res[static_cast(owner[i])], flux); + Kokkos::atomic_sub(&res[static_cast(neighbour[i])], flux); } ); parallelFor( exec, - {nInternalFaces, nBoundaryFaces}, + {nInternalFaces, nInternalFaces + nBoundaryFaces}, KOKKOS_LAMBDA(const size_t i) { auto own = static_cast(faceCells[i - nInternalFaces]); ValueType valueOwn = faceFlux[i] * phiF[i]; - Kokkos::atomic_add(&divPhi[own], valueOwn); + Kokkos::atomic_add(&res[own], valueOwn); } ); parallelFor( - exec, {0, nCells}, KOKKOS_LAMBDA(const size_t celli) { divPhi[celli] *= 1 / V[celli]; } + exec, {0, nCells}, KOKKOS_LAMBDA(const size_t celli) { res[celli] *= 1 / V[celli]; } ); } } @@ -92,9 +99,12 @@ void computeDiv( SurfaceField phif( exec, "phif", mesh, createCalculatedBCs>(mesh) ); - fill(phif.internalField(), NeoFOAM::zero::value); + // fill(phif.internalField(), NeoFOAM::zero::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( diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index 31016cc4d..bbe0a5ad4 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -26,16 +26,21 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector ); std::string execName = std::visit([](auto e) { return e.name(); }, exec); - // TODO take 1d mesh auto mesh = create1DUniformMesh(exec, 10); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + // compute correspondin uniform faceFlux fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); NeoFOAM::fill(faceFlux.internalField(), 1.0); + auto boundFaceFlux = faceFlux.internalField().span(); + // face on the left side has different orientation + boundFaceFlux[9] = -1.0; + auto volumeBCs = fvcc::createCalculatedBCs>(mesh); fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); NeoFOAM::fill(phi.internalField(), NeoFOAM::one::value); + NeoFOAM::fill(phi.boundaryField().value(), NeoFOAM::one::value); phi.correctBoundaryConditions(); auto result = NeoFOAM::Field(exec, phi.size()); From 7c0d396d436173d20a92ff43766002cf52255ff7 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 16 Feb 2025 12:53:31 +0100 Subject: [PATCH 34/35] rename test file, set correct flux value on boundary via kokkos function --- .../operator/{divOperator.cpp => gaussGreenDiv.cpp} | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) rename test/finiteVolume/cellCentred/operator/{divOperator.cpp => gaussGreenDiv.cpp} (91%) diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp similarity index 91% rename from test/finiteVolume/cellCentred/operator/divOperator.cpp rename to test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp index bbe0a5ad4..0893de394 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp @@ -29,13 +29,17 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector auto mesh = create1DUniformMesh(exec, 10); auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); - // compute correspondin uniform faceFlux + // compute corresponding uniform faceFlux + // TODO this should be handled outside of the unit test fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); NeoFOAM::fill(faceFlux.internalField(), 1.0); auto boundFaceFlux = faceFlux.internalField().span(); // face on the left side has different orientation - boundFaceFlux[9] = -1.0; - + parallelFor( + exec, + {mesh.nCells() - 1, mesh.nCells()}, + KOKKOS_LAMBDA(const size_t i) { boundFaceFlux[i] = -1.0; } + ); auto volumeBCs = fvcc::createCalculatedBCs>(mesh); fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); From 112b5929e78e4a8951610f311ca2f330e872a6e6 Mon Sep 17 00:00:00 2001 From: Gregor Olenik Date: Sun, 16 Feb 2025 13:06:13 +0100 Subject: [PATCH 35/35] clean-up --- .../cellCentred/operators/gaussGreenGrad.cpp | 44 +++++++++++-------- .../cellCentred/operator/CMakeLists.txt | 2 +- .../cellCentred/operator/gaussGreenDiv.cpp | 18 ++++---- 3 files changed, 36 insertions(+), 28 deletions(-) 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& phi, - const SurfaceInterpolation& surfInterp, - VolumeField& gradPhi + const VolumeField& in, const SurfaceInterpolation& surfInterp, VolumeField& out ) { - const UnstructuredMesh& mesh = gradPhi.mesh(); - const auto exec = gradPhi.exec(); + const UnstructuredMesh& mesh = out.mesh(); + const auto exec = out.exec(); SurfaceField phif( exec, "phif", mesh, createCalculatedBCs>(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(exec)) { for (size_t i = 0; i < nInternalFaces; i++) { - Vector flux = sSf[i] * surfPhif[i]; + Vector flux = faceAreaS[i] * surfPhif[i]; surfGradPhi[static_cast(surfOwner[i])] += flux; surfGradPhi[static_cast(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(surfOwner[i])], flux); Kokkos::atomic_sub(&surfGradPhi[static_cast(surfNeighbour[i])], flux); } @@ -70,7 +78,7 @@ void computeGrad( {nInternalFaces, surfPhif.size()}, KOKKOS_LAMBDA(const size_t i) { size_t own = static_cast(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/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 index 0893de394..a9ffe8101 100644 --- a/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp +++ b/test/finiteVolume/cellCentred/operator/gaussGreenDiv.cpp @@ -27,12 +27,12 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector std::string execName = std::visit([](auto e) { return e.name(); }, exec); auto mesh = create1DUniformMesh(exec, 10); - auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); + auto surfaceBCs = fvcc::createCalculatedBCs>(mesh); // compute corresponding uniform faceFlux // TODO this should be handled outside of the unit test - fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); - NeoFOAM::fill(faceFlux.internalField(), 1.0); + fvcc::SurfaceField faceFlux(exec, "sf", mesh, surfaceBCs); + fill(faceFlux.internalField(), 1.0); auto boundFaceFlux = faceFlux.internalField().span(); // face on the left side has different orientation parallelFor( @@ -43,22 +43,22 @@ TEMPLATE_TEST_CASE("DivOperator", "[template]", NeoFOAM::scalar, NeoFOAM::Vector auto volumeBCs = fvcc::createCalculatedBCs>(mesh); fvcc::VolumeField phi(exec, "sf", mesh, volumeBCs); - NeoFOAM::fill(phi.internalField(), NeoFOAM::one::value); - NeoFOAM::fill(phi.boundaryField().value(), NeoFOAM::one::value); + fill(phi.internalField(), one::value); + fill(phi.boundaryField().value(), one::value); phi.correctBoundaryConditions(); - auto result = NeoFOAM::Field(exec, phi.size()); - NeoFOAM::fill(result, NeoFOAM::zero::value); + auto result = Field(exec, phi.size()); + fill(result, zero::value); SECTION("Construct from Token" + execName) { - NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); + Input input = TokenList({std::string("Gauss"), std::string("linear")}); fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); } SECTION("Construct from Dictionary" + execName) { - NeoFOAM::Input input = NeoFOAM::Dictionary( + Input input = Dictionary( {{std::string("DivOperator"), std::string("Gauss")}, {std::string("surfaceInterpolation"), std::string("linear")}} );