From 467ef17b83816a6334fad650d6673e850a5652ef Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 2 Feb 2025 10:01:40 +0100 Subject: [PATCH 01/21] moved getGkoExecutor into source file --- include/NeoFOAM/linearAlgebra/utilities.hpp | 9 +-------- src/CMakeLists.txt | 1 + src/linearAlgebra/utilities.cpp | 21 +++++++++++++++++++++ 3 files changed, 23 insertions(+), 8 deletions(-) create mode 100644 src/linearAlgebra/utilities.cpp diff --git a/include/NeoFOAM/linearAlgebra/utilities.hpp b/include/NeoFOAM/linearAlgebra/utilities.hpp index 92bb0d7d3..8ce3ffed0 100644 --- a/include/NeoFOAM/linearAlgebra/utilities.hpp +++ b/include/NeoFOAM/linearAlgebra/utilities.hpp @@ -14,14 +14,7 @@ namespace NeoFOAM::la { -std::shared_ptr getGkoExecutor(Executor exec) -{ - return std::visit( - [](auto concreteExec) - { return gko::ext::kokkos::create_executor(concreteExec.underlyingExec()); }, - exec - ); -} +std::shared_ptr getGkoExecutor(Executor exec); namespace detail { diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 739d97e44..cd31ab41f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,6 +27,7 @@ target_sources( "executor/CPUExecutor.cpp" "executor/GPUExecutor.cpp" "executor/serialExecutor.cpp" + "linearAlgebra/utilities.cpp" "mesh/unstructured/boundaryMesh.cpp" "mesh/unstructured/unstructuredMesh.cpp" "finiteVolume/cellCentred/stencil/stencilDataBase.cpp" diff --git a/src/linearAlgebra/utilities.cpp b/src/linearAlgebra/utilities.cpp new file mode 100644 index 000000000..31edfb3d3 --- /dev/null +++ b/src/linearAlgebra/utilities.cpp @@ -0,0 +1,21 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2024-2025 NeoFOAM authors + +#pragma once + +#include "NeoFOAM/linearAlgebra/utilities.hpp" + + +namespace NeoFOAM::la +{ + +std::shared_ptr getGkoExecutor(Executor exec) +{ + return std::visit( + [](auto concreteExec) + { return gko::ext::kokkos::create_executor(concreteExec.underlyingExec()); }, + exec + ); +} + +} From 28996f6cb9617a99e1bcbb6c75e7488575ee066f Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 2 Feb 2025 15:04:53 +0100 Subject: [PATCH 02/21] moved utilites into source file --- CMakePresets.json | 3 +- .../operators/fvccSparsityPattern.cpp | 2 +- .../cellCentred/operators/gaussGreenDiv.cpp | 18 ++-- src/linearAlgebra/utilities.cpp | 2 - test/linearAlgebra/linearAlgebra.cpp | 2 +- test/timeIntegration/CMakeLists.txt | 1 + test/timeIntegration/forwardEuler.cpp | 91 ------------------- .../implicitTimeIntegration.cpp | 85 +++++++++++++++++ 8 files changed, 102 insertions(+), 102 deletions(-) delete mode 100644 test/timeIntegration/forwardEuler.cpp create mode 100644 test/timeIntegration/implicitTimeIntegration.cpp diff --git a/CMakePresets.json b/CMakePresets.json index 69bf55aa6..43ba62ce3 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -11,7 +11,8 @@ "hidden": true, "binaryDir": "${sourceDir}/build/${presetName}", "cacheVariables": { - "Kokkos_ENABLE_SERIAL": "ON" + "Kokkos_ENABLE_SERIAL": "ON", + "GINKGO_BUILD_PAPI_SDE": "OFF" }, "generator": "Ninja" }, diff --git a/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp b/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp index 42942be41..3cb963b72 100644 --- a/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp +++ b/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp @@ -118,7 +118,7 @@ void SparsityPattern::update() NeoFOAM::la::CSRMatrix csrMatrix(values, colIdx, rowPtrs); NeoFOAM::Field rhs(exec, nCells, 0.0); - ls_ = NeoFOAM::la::LinearSystem(csrMatrix, rhs); + ls_ = NeoFOAM::la::LinearSystem(csrMatrix, rhs, "fvcc"); } const NeoFOAM::la::LinearSystem& diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 0bf38e40c..b8418e4f9 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -142,17 +142,23 @@ void GaussGreenDiv::div( // lower triangular part - // add neighbour contribution + // add neighbour contribution upper std::size_t rowNeiStart = rowPtrs[nei]; - values[rowNeiStart + neiOffs[facei]] += -flux * weight; - Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], -flux * weight); + // values[rowNeiStart + neiOffs[facei]] += -flux * weight; + // Kokkos::atomic_add(&values[rowNeiStart + diagOffs[nei]], -flux * weight); + NeoFOAM::scalar value = flux * (1 - weight); + values[rowNeiStart + neiOffs[facei]] += flux * (1 - weight); + Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], flux * (1 - weight)); // upper triangular part - // add owner contribution + // add owner contribution lower std::size_t rowOwnStart = rowPtrs[own]; - values[rowOwnStart + ownOffs[facei]] += flux * (1 - weight); - Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], flux * (1 - weight)); + // values[rowOwnStart + ownOffs[facei]] += flux * (1 - weight); + // Kokkos::atomic_add(&values[rowOwnStart + diagOffs[own]], flux * (1 - weight)); + value = -flux * weight; + values[rowOwnStart + ownOffs[facei]] += value; + Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], value); } ); }; diff --git a/src/linearAlgebra/utilities.cpp b/src/linearAlgebra/utilities.cpp index 31edfb3d3..033013bea 100644 --- a/src/linearAlgebra/utilities.cpp +++ b/src/linearAlgebra/utilities.cpp @@ -1,8 +1,6 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2024-2025 NeoFOAM authors -#pragma once - #include "NeoFOAM/linearAlgebra/utilities.hpp" diff --git a/test/linearAlgebra/linearAlgebra.cpp b/test/linearAlgebra/linearAlgebra.cpp index 665f24781..d3323c1bd 100644 --- a/test/linearAlgebra/linearAlgebra.cpp +++ b/test/linearAlgebra/linearAlgebra.cpp @@ -53,7 +53,7 @@ TEST_CASE("MatrixAssembly - Ginkgo") NeoFOAM::la::CSRMatrix csrMatrix(values, colIdx, rowPtrs); NeoFOAM::Field rhs(exec, 3, 2.0); - NeoFOAM::la::LinearSystem linearSystem(csrMatrix, rhs); + NeoFOAM::la::LinearSystem linearSystem(csrMatrix, rhs, "custom"); NeoFOAM::Field x(exec, {0.0, 0.0, 0.0}); NeoFOAM::Dictionary solverDict; diff --git a/test/timeIntegration/CMakeLists.txt b/test/timeIntegration/CMakeLists.txt index a98fe0876..4639c08dd 100644 --- a/test/timeIntegration/CMakeLists.txt +++ b/test/timeIntegration/CMakeLists.txt @@ -2,6 +2,7 @@ # SPDX-FileCopyrightText: 2024 NeoFOAM authors neofoam_unit_test(timeIntegration) +neofoam_unit_test(implicitTimeIntegration) if(NOT WIN32) neofoam_unit_test(rungeKutta) endif() diff --git a/test/timeIntegration/forwardEuler.cpp b/test/timeIntegration/forwardEuler.cpp deleted file mode 100644 index 3300f1238..000000000 --- a/test/timeIntegration/forwardEuler.cpp +++ /dev/null @@ -1,91 +0,0 @@ -// 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 "NeoFOAM/NeoFOAM.hpp" - - -namespace dsl = NeoFOAM::DSL; - -class Divergence -{ - -public: - - std::string display() const { return "Divergence"; } - - void explicitOperation(NeoFOAM::Field& source, NeoFOAM::scalar scale) - { - auto sourceField = source.span(); - NeoFOAM::parallelFor( - source.exec(), - {0, source.size()}, - KOKKOS_LAMBDA(const size_t i) { sourceField[i] += 1.0 * scale; } - ); - } - - dsl::EqnTerm::Type getType() const { return termType_; } - - fvcc::VolumeField* volumeField() const { return nullptr; } - - const NeoFOAM::Executor& exec() const { return exec_; } - - std::size_t nCells() const { return nCells_; } - - dsl::EqnTerm::Type termType_; - NeoFOAM::Executor exec_; - std::size_t nCells_; -}; - -class TimeTerm -{ - -public: - - std::string name() const { return "TimeTerm"; } - - void explicitOperation(NeoFOAM::Field& source, NeoFOAM::scalar scale) - { - auto sourceField = source.span(); - NeoFOAM::parallelFor( - source.exec(), - {0, source.size()}, - KOKKOS_LAMBDA(const size_t i) { sourceField[i] += 1 * scale; } - ); - } - - dsl::EqnTerm::Type getType() const { return termType_; } - - fvcc::VolumeField* volumeField() const { return nullptr; } - - const NeoFOAM::Executor& exec() const { return exec_; } - - std::size_t nCells() const { return nCells_; } - - dsl::EqnTerm::Type termType_; - const NeoFOAM::Executor exec_; - std::size_t nCells_; -}; - - -TEST_CASE("TimeIntegration") -{ - auto exec = NeoFOAM::SerialExecutor(); - namespace fvcc = NeoFOAM::finiteVolume::cellCentred; - - NeoFOAM::Dictionary dict; - dict.insert("type", std::string("forwardEuler")); - - dsl::EqnTerm divTerm = Divergence(dsl::EqnTerm::Type::Explicit, exec, 1); - - dsl::EqnTerm ddtTerm = TimeTerm(dsl::EqnTerm::Type::Temporal, exec, 1); - - dsl::EqnSystem eqnSys = ddtTerm + divTerm; - - fvcc::TimeIntegration timeIntergrator(eqnSys, dict); -} diff --git a/test/timeIntegration/implicitTimeIntegration.cpp b/test/timeIntegration/implicitTimeIntegration.cpp new file mode 100644 index 000000000..a9de82ef9 --- /dev/null +++ b/test/timeIntegration/implicitTimeIntegration.cpp @@ -0,0 +1,85 @@ +// 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 "../dsl/common.hpp" + +#include "NeoFOAM/NeoFOAM.hpp" + +// only needed for msvc +template class NeoFOAM::timeIntegration::ForwardEuler; + +struct CreateField +{ + std::string name; + const NeoFOAM::UnstructuredMesh& mesh; + NeoFOAM::scalar value = 0; + std::int64_t timeIndex = 0; + std::int64_t iterationIndex = 0; + std::int64_t subCycleIndex = 0; + + NeoFOAM::Document operator()(NeoFOAM::Database& db) + { + std::vector> bcs {}; + NeoFOAM::Field internalField(mesh.exec(), mesh.nCells(), value); + fvcc::VolumeField vf( + mesh.exec(), name, mesh, internalField, bcs, db, "", "" + ); + return NeoFOAM::Document( + {{"name", vf.name}, + {"timeIndex", timeIndex}, + {"iterationIndex", iterationIndex}, + {"subCycleIndex", subCycleIndex}, + {"field", vf}}, + fvcc::validateFieldDoc + ); + } +}; + +TEST_CASE("TimeIntegration") +{ + 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::Database db; + auto mesh = NeoFOAM::createSingleCellMesh(exec); + fvcc::FieldCollection& fieldCollection = fvcc::FieldCollection::instance(db, "fieldCollection"); + + NeoFOAM::Dictionary fvSchemes; + NeoFOAM::Dictionary ddtSchemes; + ddtSchemes.insert("type", std::string("forwardEuler")); + fvSchemes.insert("ddtSchemes", ddtSchemes); + NeoFOAM::Dictionary fvSolution; + + fvcc::VolumeField& vf = + fieldCollection.registerField>( + CreateField {.name = "vf", .mesh = mesh, .value = 2.0, .timeIndex = 1} + ); + + SECTION("Create expression and perform explicitOperation on " + execName) + { + auto dummy = Dummy(vf); + Operator ddtOperator = NeoFOAM::dsl::temporal::ddt(vf); + + // ddt(U) = f + auto eqn = ddtOperator + dummy; + double dt {2.0}; + double time {1.0}; + + + // int(ddt(U)) + f = 0 + // (U^1-U^0)/dt = -f + // U^1 = - f * dt + U^0, where dt = 2, f=1, U^0=2.0 -> U^1=-2.0 + NeoFOAM::dsl::solve(eqn, vf, time, dt, fvSchemes, fvSolution); + REQUIRE(getField(vf.internalField()) == -2.0); + } +} From fab9ebda1be4a8717fd78323f7306629eb6b97c1 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 2 Feb 2025 15:05:14 +0100 Subject: [PATCH 03/21] added linearSytem for fvcc --- .../operators/fvccLinearSystem.hpp | 138 ++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp new file mode 100644 index 000000000..c3af8cc47 --- /dev/null +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp @@ -0,0 +1,138 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2025 NeoFOAM authors + +#pragma once + +#include "NeoFOAM/linearAlgebra.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" + +namespace NeoFOAM::finiteVolume::cellCentred +{ + +template +class LinearSystem +{ +public: + + LinearSystem(VolumeField& psi) + : psi_(psi), ls_(SparsityPattern::readOrCreate(psi.mesh()).linearSystem()), + sparsityPattern_(SparsityPattern::readOrCreate(psi.mesh())) { + + }; + + LinearSystem( + VolumeField& psi, + const la::LinearSystem& ls, + std::shared_ptr sparsityPattern + ) + : psi_(psi), ls_(ls), sparsityPattern_(sparsityPattern) {}; + + LinearSystem(const LinearSystem& ls) + : psi_(ls.psi_), ls_(ls.ls_), sparsityPattern_(ls.sparsityPattern_) {}; + + ~LinearSystem() = default; + + [[nodiscard]] la::LinearSystem& linearSystem() { return ls_; } + [[nodiscard]] SparsityPattern& sparsityPattern() + { + if (!sparsityPattern_) + { + NF_THROW("fvcc:LinearSystem:sparsityPattern: sparsityPattern is null"); + } + return *sparsityPattern_; + } + + [[nodiscard]] const la::LinearSystem& linearSystem() const { return ls_; } + [[nodiscard]] const SparsityPattern& sparsityPattern() const + { + if (!sparsityPattern_) + { + NF_THROW("fvcc:LinearSystem:sparsityPattern: sparsityPattern is null"); + } + return *sparsityPattern_; + } + + const Executor& exec() const { return ls_.exec(); } + + void diag(Field& field) + { + if (field.size() != sparsityPattern_->diagOffset().size()) + { + NF_THROW("fvcc:LinearSystem:diag: field size mismatch"); + } + const auto diagOffset = sparsityPattern_->diagOffset().span(); + const auto rowPtrs = ls_.matrix().rowPtrs(); + std::span fieldSpan = field.span(); + std::span values = ls_.matrix().values(); + NeoFOAM::parallelFor( + exec(), + {0, diagOffset.size()}, + KOKKOS_LAMBDA(const std::size_t celli) { + auto diagOffsetCelli = diagOffset[celli]; + fieldSpan[celli] = values[rowPtrs[celli] + diagOffsetCelli]; + } + ); + } + + Field diagIndex() + { + Field diagIndex(exec(), sparsityPattern_->diagOffset().size()); + const auto diagOffset = sparsityPattern_->diagOffset().span(); + auto diagIndexSpan = diagIndex.span(); + const auto rowPtrs = ls_.matrix().rowPtrs(); + NeoFOAM::parallelFor( + exec(), + {0, diagIndex.size()}, + KOKKOS_LAMBDA(const std::size_t celli) { + diagIndexSpan[celli] = rowPtrs[celli] + diagOffset[celli]; + } + ); + return diagIndex; + } + +private: + + VolumeField& psi_; + la::LinearSystem ls_; + std::shared_ptr sparsityPattern_; +}; + +template +VolumeField +operator&(const LinearSystem ls, const VolumeField& psi) +{ + VolumeField resultField( + psi.exec(), + "ls_" + psi.name, + psi.mesh(), + psi.internalField(), + psi.boundaryField(), + psi.boundaryConditions() + ); + + auto [result, b, x] = + spans(resultField.internalField(), ls.linearSystem().rhs(), psi.internalField()); + + const std::span values = ls.linearSystem().matrix().values(); + const std::span colIdxs = ls.linearSystem().matrix().colIdxs(); + const std::span rowPtrs = ls.linearSystem().matrix().rowPtrs(); + + NeoFOAM::parallelFor( + resultField.exec(), + {0, result.size()}, + KOKKOS_LAMBDA(const std::size_t rowi) { + IndexType rowStart = rowPtrs[rowi]; + IndexType rowEnd = rowPtrs[rowi + 1]; + ValueType sum = 0.0; + for (IndexType coli = rowStart; coli < rowEnd; coli++) + { + sum += values[coli] * x[colIdxs[coli]]; + } + result[rowi] = sum - b[rowi]; + } + ); + + return resultField; +} + +} From 7d4544366ca89f45ce1b6e92fd8f7a7305975640 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Thu, 6 Feb 2025 20:37:48 +0100 Subject: [PATCH 04/21] fixup! --- .../finiteVolume/cellCentred/operators/fvccLinearSystem.hpp | 2 +- src/linearAlgebra/utilities.cpp | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp index c3af8cc47..29c8959b6 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp @@ -15,7 +15,7 @@ class LinearSystem public: LinearSystem(VolumeField& psi) - : psi_(psi), ls_(SparsityPattern::readOrCreate(psi.mesh()).linearSystem()), + : psi_(psi), ls_(SparsityPattern::readOrCreate(psi.mesh())->linearSystem()), sparsityPattern_(SparsityPattern::readOrCreate(psi.mesh())) { }; diff --git a/src/linearAlgebra/utilities.cpp b/src/linearAlgebra/utilities.cpp index 31edfb3d3..033013bea 100644 --- a/src/linearAlgebra/utilities.cpp +++ b/src/linearAlgebra/utilities.cpp @@ -1,8 +1,6 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2024-2025 NeoFOAM authors -#pragma once - #include "NeoFOAM/linearAlgebra/utilities.hpp" From aace259c273f82d142765ba57d6721a56e58003a Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sat, 8 Feb 2025 19:12:51 +0100 Subject: [PATCH 05/21] solve expression with ginkgo --- include/NeoFOAM/dsl/solver.hpp | 54 ++++++++++++++++++- .../cellCentred/fields/geometricField.hpp | 3 ++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index b0cb0f7d4..48cbe571c 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -14,10 +14,57 @@ #include "NeoFOAM/dsl/expression.hpp" #include "NeoFOAM/timeIntegration/timeIntegration.hpp" +#include "NeoFOAM/linearAlgebra/ginkgo.hpp" + namespace NeoFOAM::dsl { +template +NeoFOAM::la::LinearSystem +ginkgoMatrix(Expression& exp, FieldType& solution) +{ + using ValueType = typename FieldType::ElementType; + auto vol = solution.mesh().cellVolumes().span(); + // la::CSRMatrix matrix(values, colIdxs, rowPtrs); + // return la::LinearSystem(matrix, rhs, sparsityPattern); + auto expSource = exp.explicitOperation(solution.mesh().nCells()); + auto expSourceSpan = expSource.span(); + + auto ls = exp.implicitOperation(); + Field rhs(solution.exec(), ls.rhs().data(), ls.rhs().size()); + auto rhsSpan = rhs.span(); + // we subtract the explicit source term from the rhs + NeoFOAM::parallelFor( + solution.exec(), + {0, rhs.size()}, + KOKKOS_LAMBDA(const size_t i) { rhsSpan[i] -= expSourceSpan[i] * vol[i]; } + ); + + Field mValues( + solution.exec(), ls.matrix().values().data(), ls.matrix().values().size() + ); + Field mColIdxs(solution.exec(), ls.matrix().colIdxs().size()); + auto mColIdxsSpan = ls.matrix().colIdxs(); + NeoFOAM::parallelFor( + mColIdxs, KOKKOS_LAMBDA(const size_t i) { return int(mColIdxsSpan[i]); } + ); + + Field mRowPtrs(solution.exec(), ls.matrix().rowPtrs().size()); + auto mRowPtrsSpan = ls.matrix().rowPtrs(); + NeoFOAM::parallelFor( + mRowPtrs, KOKKOS_LAMBDA(const size_t i) { return int(mRowPtrsSpan[i]); } + ); + + auto values = ls.matrix().values(); + + + la::CSRMatrix matrix(mValues, mColIdxs, mRowPtrs); + + NeoFOAM::la::LinearSystem convertedLs(matrix, rhs, ls.sparsityPattern()); + return convertedLs; +} + /* @brief solve an expression * * @param exp - Expression which is to be solved/updated. @@ -50,8 +97,13 @@ void solve( } else { - NF_ERROR_EXIT("Not implemented."); // solve sparse matrix system + using ValueType = typename FieldType::ElementType; + auto ls = ginkgoMatrix(exp, solution); + + + NeoFOAM::la::ginkgo::BiCGStab solver(solution.exec(), fvSolution); + solver.solve(ls, solution.internalField()); } } diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp index 89e6c2b3d..655c9f63f 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp @@ -27,6 +27,9 @@ class GeometricFieldMixin { public: + + typedef ValueType ElementType; + /** * @brief Constructor for GeometricFieldMixin. * From 0da0f1d3bcc1a796fee197c1dc67a9fb7b13a80c Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sat, 8 Feb 2025 19:13:24 +0100 Subject: [PATCH 06/21] add sourceTerm --- include/NeoFOAM/dsl/explicit.hpp | 6 ++++++ .../finiteVolume/cellCentred/operators/sourceTerm.hpp | 4 +++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/include/NeoFOAM/dsl/explicit.hpp b/include/NeoFOAM/dsl/explicit.hpp index e2fe0297e..a56e7a25a 100644 --- a/include/NeoFOAM/dsl/explicit.hpp +++ b/include/NeoFOAM/dsl/explicit.hpp @@ -12,6 +12,7 @@ // TODO we should get rid of this include since it includes details // from a general implementation #include "NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp" namespace fvcc = NeoFOAM::finiteVolume::cellCentred; @@ -24,5 +25,10 @@ div(const fvcc::SurfaceField& faceFlux, fvcc::VolumeField& coeff, fvcc::VolumeField& phi) +{ + return Operator(fvcc::SourceTerm(dsl::Operator::Type::Explicit, coeff, phi)); +} + } // namespace NeoFOAM diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp index 1aa4d9378..cead24c27 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp @@ -30,6 +30,7 @@ class SourceTerm : public dsl::OperatorMixin> void explicitOperation(Field& source) { auto operatorScaling = getCoefficient(); + const auto vol = coefficients_.mesh().cellVolumes().span(); auto [sourceSpan, fieldSpan, coeff] = spans(source, field_.internalField(), coefficients_.internalField()); NeoFOAM::parallelFor( @@ -49,6 +50,7 @@ class SourceTerm : public dsl::OperatorMixin> void implicitOperation(la::LinearSystem& ls) { const auto operatorScaling = getCoefficient(); + const auto vol = coefficients_.mesh().cellVolumes().span(); const auto [diagOffs, coeff] = spans(sparsityPattern_->diagOffset(), coefficients_.internalField()); const auto rowPtrs = ls.matrix().rowPtrs(); @@ -59,7 +61,7 @@ class SourceTerm : public dsl::OperatorMixin> {0, coeff.size()}, KOKKOS_LAMBDA(const size_t celli) { std::size_t idx = rowPtrs[celli] + diagOffs[celli]; - values[idx] += operatorScaling[celli] * coeff[celli]; + values[idx] += operatorScaling[celli] * coeff[celli] * vol[celli]; } ); } From e2d69b313102c7c3438a9b2b43de5953cbf52a71 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sat, 8 Feb 2025 19:14:07 +0100 Subject: [PATCH 07/21] WIP ddtOperator --- include/NeoFOAM/dsl/implicit.hpp | 10 ++++++++++ .../cellCentred/operators/ddtOperator.hpp | 12 +++++++++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/include/NeoFOAM/dsl/implicit.hpp b/include/NeoFOAM/dsl/implicit.hpp index 0a19d9df3..0e1a76cde 100644 --- a/include/NeoFOAM/dsl/implicit.hpp +++ b/include/NeoFOAM/dsl/implicit.hpp @@ -19,4 +19,14 @@ namespace NeoFOAM::dsl::imp Operator ddt(fvcc::VolumeField& phi) { return dsl::temporal::ddt(phi); } +Operator Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) +{ + return Operator(fvcc::SourceTerm(dsl::Operator::Type::Implicit, coeff, phi)); +} + +Operator div(fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) +{ + return Operator(fvcc::DivOperator(dsl::Operator::Type::Implicit, faceFlux, phi)); +} + } // namespace NeoFOAM diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index f4ad512bc..a835402a0 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -22,7 +22,9 @@ class DdtOperator : public dsl::OperatorMixin> DdtOperator(dsl::Operator::Type termType, VolumeField& field) : dsl::OperatorMixin>(field.exec(), field, termType), - sparsityPattern_(SparsityPattern::readOrCreate(field.mesh())) {}; + sparsityPattern_(SparsityPattern::readOrCreate(field.mesh())) { + + }; void explicitOperation(Field& source) { @@ -45,17 +47,21 @@ class DdtOperator : public dsl::OperatorMixin> void implicitOperation(la::LinearSystem& ls) { + const scalar rDeltat = 1 / 1.0; const auto operatorScaling = getCoefficient(); const auto [diagOffs, oldField] = spans(sparsityPattern_->diagOffset(), oldTime(field_).internalField()); const auto rowPtrs = ls.matrix().rowPtrs(); const auto colIdxs = ls.matrix().colIdxs(); std::span values = ls.matrix().values(); + std::span rhs = ls.rhs().span(); NeoFOAM::parallelFor( ls.exec(), - {0, coeff.size()}, + {0, oldField.size()}, KOKKOS_LAMBDA(const size_t celli) { std::size_t idx = rowPtrs[celli] + diagOffs[celli]; + values[idx] += operatorScaling[celli] * rDeltat; + rhs[celli] += operatorScaling[celli] * rDeltat * oldField[celli]; } ); } @@ -70,7 +76,7 @@ class DdtOperator : public dsl::OperatorMixin> private: - const VolumeField& coefficients_; + // const VolumeField& coefficients_; const std::shared_ptr sparsityPattern_; }; From 05dbdcab0dc453c1b31e430341f5f85e71d71c4d Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sat, 8 Feb 2025 19:14:34 +0100 Subject: [PATCH 08/21] implicit Operator for GaussGreenDiv --- include/NeoFOAM/finiteVolume/cellCentred.hpp | 2 ++ .../cellCentred/operators/divOperator.hpp | 20 ++++++++++++++ .../cellCentred/operators/gaussGreenDiv.hpp | 2 ++ .../cellCentred/operators/gaussGreenDiv.cpp | 27 ++++++++++--------- 4 files changed, 39 insertions(+), 12 deletions(-) diff --git a/include/NeoFOAM/finiteVolume/cellCentred.hpp b/include/NeoFOAM/finiteVolume/cellCentred.hpp index 665af19db..70ad29da5 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred.hpp @@ -14,6 +14,8 @@ #include "cellCentred/operators/divOperator.hpp" #include "cellCentred/operators/gaussGreenDiv.hpp" +#include "cellCentred/operators/ddtOperator.hpp" +#include "cellCentred/operators/sourceTerm.hpp" #include "cellCentred/interpolation/linear.hpp" #include "cellCentred/interpolation/upwind.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index 80419b81d..5b489cbe4 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -42,6 +42,8 @@ class DivOperatorFactory : virtual ~DivOperatorFactory() {} // Virtual destructor + virtual la::LinearSystem createEmptyLinearSystem() const = 0; + virtual void div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) = 0; @@ -116,6 +118,24 @@ class DivOperator : public dsl::OperatorMixin> source += tmpsource; } + la::LinearSystem createEmptyLinearSystem() const + { + if (divOperatorStrategy_ == nullptr) + { + NF_ERROR_EXIT("DivOperatorStrategy not initialized"); + } + return divOperatorStrategy_->createEmptyLinearSystem(); + } + + void implicitOperation(la::LinearSystem& ls) + { + if (divOperatorStrategy_ == nullptr) + { + NF_ERROR_EXIT("DivOperatorStrategy not initialized"); + } + divOperatorStrategy_->div(ls, faceFlux_, field_); + } + void div(Field& divPhi) { divOperatorStrategy_->div(divPhi, faceFlux_, getField()); } void div(la::LinearSystem& ls) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp index e1a67033f..eca993164 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp @@ -25,6 +25,8 @@ class GaussGreenDiv : public DivOperatorFactory::Register GaussGreenDiv(const Executor& exec, const UnstructuredMesh& mesh, const Input& inputs); + la::LinearSystem createEmptyLinearSystem() const override; + void div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) override; diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index b8418e4f9..b7eb58462 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -103,6 +103,12 @@ GaussGreenDiv::GaussGreenDiv( surfaceInterpolation_(exec, mesh, inputs), sparsityPattern_(SparsityPattern::readOrCreate(mesh)) {}; + +la::LinearSystem GaussGreenDiv::createEmptyLinearSystem() const +{ + return sparsityPattern_->linearSystem(); +}; + void GaussGreenDiv::div( VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ) @@ -137,28 +143,25 @@ void GaussGreenDiv::div( KOKKOS_LAMBDA(const size_t facei) { scalar flux = sFaceFlux[facei]; scalar weight = 0.5; + scalar value = 0; std::size_t own = static_cast(owner[facei]); std::size_t nei = static_cast(neighbour[facei]); - // lower triangular part - // add neighbour contribution upper std::size_t rowNeiStart = rowPtrs[nei]; - // values[rowNeiStart + neiOffs[facei]] += -flux * weight; - // Kokkos::atomic_add(&values[rowNeiStart + diagOffs[nei]], -flux * weight); - NeoFOAM::scalar value = flux * (1 - weight); - values[rowNeiStart + neiOffs[facei]] += flux * (1 - weight); - Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], flux * (1 - weight)); + std::size_t rowOwnStart = rowPtrs[own]; + + value = -weight * flux; + // scalar valueNei = (1 - weight) * flux; + values[rowNeiStart + neiOffs[facei]] += value; + Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], value); // upper triangular part // add owner contribution lower - std::size_t rowOwnStart = rowPtrs[own]; - // values[rowOwnStart + ownOffs[facei]] += flux * (1 - weight); - // Kokkos::atomic_add(&values[rowOwnStart + diagOffs[own]], flux * (1 - weight)); - value = -flux * weight; + value = flux * (1 - weight); values[rowOwnStart + ownOffs[facei]] += value; - Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], value); + Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], value); } ); }; From b0951ee25a22a80b1c683d70bac63a5a34ae49a9 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 9 Feb 2025 00:04:07 +0100 Subject: [PATCH 09/21] replace Operator with Spatial and TemporalOperator --- include/NeoFOAM/dsl.hpp | 3 +- include/NeoFOAM/dsl/ddt.hpp | 14 +- include/NeoFOAM/dsl/explicit.hpp | 11 +- include/NeoFOAM/dsl/expression.hpp | 119 ++++----- include/NeoFOAM/dsl/implicit.hpp | 18 +- include/NeoFOAM/dsl/solver.hpp | 3 +- .../dsl/{operator.hpp => spatialOperator.hpp} | 60 ++--- include/NeoFOAM/dsl/temporalOperator.hpp | 231 ++++++++++++++++++ .../cellCentred/operators/ddtOperator.hpp | 13 +- .../cellCentred/operators/divOperator.hpp | 10 +- .../cellCentred/operators/sourceTerm.hpp | 6 +- src/CMakeLists.txt | 3 +- src/dsl/operator.cpp | 60 ----- src/dsl/spatialOperator.cpp | 67 +++++ src/dsl/temporalOperator.cpp | 70 ++++++ test/dsl/CMakeLists.txt | 3 +- test/dsl/common.hpp | 84 ++++++- test/dsl/expression.cpp | 8 +- .../dsl/{operator.cpp => spatialOperator.cpp} | 29 ++- test/dsl/temporalOperator.cpp | 137 +++++++++++ .../cellCentred/operator/divOperator.cpp | 6 +- .../implicitTimeIntegration.cpp | 2 +- test/timeIntegration/rungeKutta.cpp | 10 +- test/timeIntegration/timeIntegration.cpp | 2 +- 24 files changed, 747 insertions(+), 222 deletions(-) rename include/NeoFOAM/dsl/{operator.hpp => spatialOperator.hpp} (83%) create mode 100644 include/NeoFOAM/dsl/temporalOperator.hpp delete mode 100644 src/dsl/operator.cpp create mode 100644 src/dsl/spatialOperator.cpp create mode 100644 src/dsl/temporalOperator.cpp rename test/dsl/{operator.cpp => spatialOperator.cpp} (80%) create mode 100644 test/dsl/temporalOperator.cpp diff --git a/include/NeoFOAM/dsl.hpp b/include/NeoFOAM/dsl.hpp index 50a2941c0..6a588a36e 100644 --- a/include/NeoFOAM/dsl.hpp +++ b/include/NeoFOAM/dsl.hpp @@ -6,5 +6,6 @@ #include "dsl/explicit.hpp" #include "dsl/expression.hpp" #include "dsl/implicit.hpp" -#include "dsl/operator.hpp" +#include "dsl/spatialOperator.hpp" +#include "dsl/temporalOperator.hpp" #include "dsl/solver.hpp" diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp index 992596fe3..583f985c5 100644 --- a/include/NeoFOAM/dsl/ddt.hpp +++ b/include/NeoFOAM/dsl/ddt.hpp @@ -5,7 +5,7 @@ #pragma once #include "NeoFOAM/fields/field.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" namespace NeoFOAM::dsl::temporal { @@ -16,17 +16,23 @@ class Ddt : public OperatorMixin public: - Ddt(FieldType& field) : OperatorMixin(field.exec(), field, Operator::Type::Temporal) + Ddt(FieldType& field) + : OperatorMixin(field.exec(), field, SpatialOperator::Type::Implicit) {} std::string getName() const { return "TimeOperator"; } - void explicitOperation([[maybe_unused]] Field& source, [[maybe_unused]] scalar scale) + void + explicitOperation([[maybe_unused]] Field& source, NeoFOAM::scalar t, NeoFOAM::scalar dt) { NF_ERROR_EXIT("Not implemented"); } - void implicitOperation([[maybe_unused]] Field& phi) + void implicitOperation( + la::LinearSystem& ls, + NeoFOAM::scalar t, + NeoFOAM::scalar dt + ) { NF_ERROR_EXIT("Not implemented"); } diff --git a/include/NeoFOAM/dsl/explicit.hpp b/include/NeoFOAM/dsl/explicit.hpp index a56e7a25a..eb3c423ff 100644 --- a/include/NeoFOAM/dsl/explicit.hpp +++ b/include/NeoFOAM/dsl/explicit.hpp @@ -5,7 +5,7 @@ #pragma once #include "NeoFOAM/fields/field.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" #include "NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp" @@ -19,15 +19,16 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; namespace NeoFOAM::dsl::exp { -Operator +SpatialOperator div(const fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return Operator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi)); + return SpatialOperator(fvcc::DivOperator(dsl::SpatialOperator::Type::Explicit, faceFlux, phi)); } -Operator Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) +SpatialOperator +Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) { - return Operator(fvcc::SourceTerm(dsl::Operator::Type::Explicit, coeff, phi)); + return SpatialOperator(fvcc::SourceTerm(dsl::SpatialOperator::Type::Explicit, coeff, phi)); } diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp index f67484808..d2ff1ad2e 100644 --- a/include/NeoFOAM/dsl/expression.hpp +++ b/include/NeoFOAM/dsl/expression.hpp @@ -9,7 +9,8 @@ #include "NeoFOAM/core/primitives/scalar.hpp" #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/linearAlgebra/linearSystem.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" +#include "NeoFOAM/dsl/temporalOperator.hpp" #include "NeoFOAM/core/error.hpp" namespace la = NeoFOAM::la; @@ -22,13 +23,11 @@ class Expression { public: - Expression(const Executor& exec) - : exec_(exec), temporalOperators_(), implicitOperators_(), explicitOperators_() - {} + Expression(const Executor& exec) : exec_(exec), temporalOperators_(), spatialOperators_() {} Expression(const Expression& exp) : exec_(exp.exec_), temporalOperators_(exp.temporalOperators_), - implicitOperators_(exp.implicitOperators_), explicitOperators_(exp.explicitOperators_) + spatialOperators_(exp.spatialOperators_) {} void build(const NeoFOAM::Dictionary& input) @@ -37,11 +36,7 @@ class Expression { op.build(input); } - for (auto& op : implicitOperators_) - { - op.build(input); - } - for (auto& op : explicitOperators_) + for (auto& op : spatialOperators_) { op.build(input); } @@ -57,43 +52,45 @@ class Expression /* @brief perform all explicit operation and accumulate the result */ Field explicitOperation(Field& source) { - for (auto& oper : explicitOperators_) + for (auto& oper : spatialOperators_) { - oper.explicitOperation(source); + if (oper.getType() == SpatialOperator::Type::Explicit) + { + oper.explicitOperation(source); + } } + // for (auto& oper : temporalOperators_) + // { + // if (oper.getType() == SpatialOperator::Type::Explicit) + // { + // oper.explicitOperation(source); + // } + // } return source; } /* @brief perform all implicit operation and accumulate the result */ la::LinearSystem implicitOperation() { - if (implicitOperators_.empty()) + // TODO: implement + // if (implicitOperators_.empty()) + // { + // NF_ERROR_EXIT("No implicit operators in the expression"); + // } + auto ls = spatialOperators_[0].createEmptyLinearSystem(); + for (auto& oper : spatialOperators_) { - NF_ERROR_EXIT("No implicit operators in the expression"); - } - auto ls = implicitOperators_[0].createEmptyLinearSystem(); - for (auto& oper : implicitOperators_) - { - oper.implicitOperation(ls); + if (oper.getType() == SpatialOperator::Type::Implicit) + { + oper.implicitOperation(ls); + } } return ls; } - void addOperator(const Operator& oper) - { - switch (oper.getType()) - { - case Operator::Type::Temporal: - temporalOperators_.push_back(oper); - break; - case Operator::Type::Implicit: - implicitOperators_.push_back(oper); - break; - case Operator::Type::Explicit: - explicitOperators_.push_back(oper); - break; - } - } + void addOperator(const SpatialOperator& oper) { spatialOperators_.push_back(oper); } + + void addOperator(const TemporalOperator& oper) { temporalOperators_.push_back(oper); } void addExpression(const Expression& equation) { @@ -101,35 +98,24 @@ class Expression { temporalOperators_.push_back(oper); } - for (auto& oper : equation.implicitOperators_) + for (auto& oper : equation.spatialOperators_) { - implicitOperators_.push_back(oper); - } - for (auto& oper : equation.explicitOperators_) - { - explicitOperators_.push_back(oper); + spatialOperators_.push_back(oper); } } /* @brief getter for the total number of terms in the equation */ - size_t size() const - { - return temporalOperators_.size() + implicitOperators_.size() + explicitOperators_.size(); - } + size_t size() const { return temporalOperators_.size() + spatialOperators_.size(); } // getters - const std::vector& temporalOperators() const { return temporalOperators_; } - - const std::vector& implicitOperators() const { return implicitOperators_; } + const std::vector& temporalOperators() const { return temporalOperators_; } - const std::vector& explicitOperators() const { return explicitOperators_; } + const std::vector& spatialOperators() const { return spatialOperators_; } - std::vector& temporalOperators() { return temporalOperators_; } + std::vector& temporalOperators() { return temporalOperators_; } - std::vector& implicitOperators() { return implicitOperators_; } - - std::vector& explicitOperators() { return explicitOperators_; } + std::vector& spatialOperators() { return spatialOperators_; } const Executor& exec() const { return exec_; } @@ -137,11 +123,9 @@ class Expression const Executor exec_; - std::vector temporalOperators_; - - std::vector implicitOperators_; + std::vector temporalOperators_; - std::vector explicitOperators_; + std::vector spatialOperators_; }; [[nodiscard]] inline Expression operator+(Expression lhs, const Expression& rhs) @@ -150,13 +134,22 @@ class Expression return lhs; } -[[nodiscard]] inline Expression operator+(Expression lhs, const Operator& rhs) +[[nodiscard]] inline Expression operator+(Expression lhs, const SpatialOperator& rhs) { lhs.addOperator(rhs); return lhs; } -[[nodiscard]] inline Expression operator+(const Operator& lhs, const Operator& rhs) +template +[[nodiscard]] inline Expression operator+(leftOperator lhs, rightOperator rhs) +{ + Expression expr(lhs.exec()); + expr.addOperator(lhs); + expr.addOperator(rhs); + return expr; +} + +[[nodiscard]] inline Expression operator+(const SpatialOperator& lhs, const SpatialOperator& rhs) { Expression expr(lhs.exec()); expr.addOperator(lhs); @@ -171,11 +164,7 @@ class Expression { expr.addOperator(scale * oper); } - for (const auto& oper : es.implicitOperators()) - { - expr.addOperator(scale * oper); - } - for (const auto& oper : es.explicitOperators()) + for (const auto& oper : es.spatialOperators()) { expr.addOperator(scale * oper); } @@ -188,13 +177,13 @@ class Expression return lhs; } -[[nodiscard]] inline Expression operator-(Expression lhs, const Operator& rhs) +[[nodiscard]] inline Expression operator-(Expression lhs, const SpatialOperator& rhs) { lhs.addOperator(-1.0 * rhs); return lhs; } -[[nodiscard]] inline Expression operator-(const Operator& lhs, const Operator& rhs) +[[nodiscard]] inline Expression operator-(const SpatialOperator& lhs, const SpatialOperator& rhs) { Expression expr(lhs.exec()); expr.addOperator(lhs); diff --git a/include/NeoFOAM/dsl/implicit.hpp b/include/NeoFOAM/dsl/implicit.hpp index 0e1a76cde..06bbcf8bc 100644 --- a/include/NeoFOAM/dsl/implicit.hpp +++ b/include/NeoFOAM/dsl/implicit.hpp @@ -5,7 +5,8 @@ #pragma once #include "NeoFOAM/fields/field.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" +#include "NeoFOAM/dsl/temporalOperator.hpp" #include "NeoFOAM/dsl/ddt.hpp" #include "NeoFOAM/finiteVolume/cellCentred.hpp" #include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" @@ -17,16 +18,21 @@ namespace NeoFOAM::dsl::imp { -Operator ddt(fvcc::VolumeField& phi) { return dsl::temporal::ddt(phi); } +TemporalOperator ddt(fvcc::VolumeField& phi) +{ + return fvcc::DdtOperator(dsl::SpatialOperator::Type::Implicit, phi); +} -Operator Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) +SpatialOperator +Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) { - return Operator(fvcc::SourceTerm(dsl::Operator::Type::Implicit, coeff, phi)); + return SpatialOperator(fvcc::SourceTerm(dsl::SpatialOperator::Type::Implicit, coeff, phi)); } -Operator div(fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) +SpatialOperator +div(fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return Operator(fvcc::DivOperator(dsl::Operator::Type::Implicit, faceFlux, phi)); + return SpatialOperator(fvcc::DivOperator(dsl::SpatialOperator::Type::Implicit, faceFlux, phi)); } } // namespace NeoFOAM diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index 48cbe571c..38eaab720 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -84,7 +84,8 @@ void solve( [[maybe_unused]] const Dictionary& fvSolution ) { - if (exp.temporalOperators().size() == 0 && exp.implicitOperators().size() == 0) + // todo fix + if (exp.temporalOperators().size() == 0 && exp.spatialOperators().size() == 0) { NF_ERROR_EXIT("No temporal or implicit terms to solve."); } diff --git a/include/NeoFOAM/dsl/operator.hpp b/include/NeoFOAM/dsl/spatialOperator.hpp similarity index 83% rename from include/NeoFOAM/dsl/operator.hpp rename to include/NeoFOAM/dsl/spatialOperator.hpp index a597f6e79..6a84a6750 100644 --- a/include/NeoFOAM/dsl/operator.hpp +++ b/include/NeoFOAM/dsl/spatialOperator.hpp @@ -16,13 +16,6 @@ namespace la = NeoFOAM::la; namespace NeoFOAM::dsl { -template -concept HasTemporalOperator = requires(T t) { - { - t.temporalOperation(std::declval&>(), std::declval()) - } -> std::same_as; // Adjust return type and arguments as needed -}; - template concept HasExplicitOperator = requires(T t) { { @@ -37,6 +30,9 @@ concept HasImplicitOperator = requires(T t) { } -> std::same_as; // Adjust return type and arguments as needed }; +template +concept IsSpatialOperator = HasExplicitOperator || HasImplicitOperator; + /* @class Operator * @brief A class to represent an operator in NeoFOAMs dsl * @@ -49,7 +45,7 @@ concept HasImplicitOperator = requires(T t) { * * @ingroup dsl */ -class Operator +class SpatialOperator { public: @@ -60,24 +56,24 @@ class Operator Explicit }; - template - Operator(T cls) : model_(std::make_unique>(std::move(cls))) + template + SpatialOperator(T cls) : model_(std::make_unique>(std::move(cls))) {} - Operator(const Operator& eqnOperator); + SpatialOperator(const SpatialOperator& eqnOperator); - Operator(Operator&& eqnOperator); + SpatialOperator(SpatialOperator&& eqnOperator); void explicitOperation(Field& source); - void temporalOperation(Field& field); - void implicitOperation(la::LinearSystem& ls); la::LinearSystem createEmptyLinearSystem() const; - /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - Operator::Type getType() const; + /* returns the fundamental type of an operator, ie explicit, implicit */ + SpatialOperator::Type getType() const; + + std::string getName() const; Coeff& getCoefficient(); @@ -101,8 +97,6 @@ class Operator virtual void explicitOperation(Field& source) = 0; - virtual void temporalOperation(Field& field) = 0; - virtual void implicitOperation(la::LinearSystem& ls) = 0; virtual la::LinearSystem createEmptyLinearSystem() const = 0; @@ -113,8 +107,8 @@ class Operator /* returns the name of the operator */ virtual std::string getName() const = 0; - /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - virtual Operator::Type getType() const = 0; + /* returns the fundamental type of an operator, ie explicit, implicit */ + virtual SpatialOperator::Type getType() const = 0; /* @brief get the associated coefficient for this term */ virtual Coeff& getCoefficient() = 0; @@ -147,14 +141,6 @@ class Operator } } - virtual void temporalOperation(Field& field) override - { - if constexpr (HasTemporalOperator) - { - concreteOp_.temporalOperation(field); - } - } - virtual void implicitOperation(la::LinearSystem& ls) override { if constexpr (HasImplicitOperator) @@ -188,7 +174,7 @@ class Operator virtual void build(const Input& input) override { concreteOp_.build(input); } /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - Operator::Type getType() const override { return concreteOp_.getType(); } + SpatialOperator::Type getType() const override { return concreteOp_.getType(); } /* @brief Get the executor */ const Executor& exec() const override { return concreteOp_.exec(); } @@ -212,19 +198,19 @@ class Operator }; -Operator operator*(scalar scalarCoeff, Operator rhs); +SpatialOperator operator*(scalar scalarCoeff, SpatialOperator rhs); -Operator operator*(const Field& coeffField, Operator rhs); +SpatialOperator operator*(const Field& coeffField, SpatialOperator rhs); -Operator operator*(const Coeff& coeff, Operator rhs); +SpatialOperator operator*(const Coeff& coeff, SpatialOperator rhs); template requires std::invocable -Operator operator*([[maybe_unused]] CoeffFunction coeffFunc, const Operator& lhs) +SpatialOperator operator*([[maybe_unused]] CoeffFunction coeffFunc, const SpatialOperator& lhs) { // TODO implement NF_ERROR_EXIT("Not implemented"); - Operator result = lhs; + SpatialOperator result = lhs; // if (!result.getCoefficient().useSpan) // { // result.setField(std::make_shared>(result.exec(), result.nCells(), 1.0)); @@ -245,10 +231,10 @@ class OperatorMixin public: - OperatorMixin(const Executor exec, FieldType& field, Operator::Type type) + OperatorMixin(const Executor exec, FieldType& field, SpatialOperator::Type type) : exec_(exec), coeffs_(), field_(field), type_(type) {}; - Operator::Type getType() const { return type_; } + SpatialOperator::Type getType() const { return type_; } virtual ~OperatorMixin() = default; @@ -273,6 +259,6 @@ class OperatorMixin FieldType& field_; - Operator::Type type_; + SpatialOperator::Type type_; }; } // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/temporalOperator.hpp b/include/NeoFOAM/dsl/temporalOperator.hpp new file mode 100644 index 000000000..47c3b604a --- /dev/null +++ b/include/NeoFOAM/dsl/temporalOperator.hpp @@ -0,0 +1,231 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#pragma once + +#include +#include + +#include "NeoFOAM/core/primitives/scalar.hpp" +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/linearAlgebra/linearSystem.hpp" +#include "NeoFOAM/core/input.hpp" +#include "NeoFOAM/dsl/coeff.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" + +namespace la = NeoFOAM::la; + +namespace NeoFOAM::dsl +{ + + +template +concept HasTemporalExplicitOperator = requires(T t) { + { + t.explicitOperation( + std::declval&>(), + std::declval(), + std::declval() + ) + } -> std::same_as; // Adjust return type and arguments as needed +}; + +template +concept HasTemporalImplicitOperator = requires(T t) { + { + t.implicitOperation( + std::declval&>(), + std::declval(), + std::declval() + ) + } -> std::same_as; // Adjust return type and arguments as needed +}; + +template +concept HasTemporalOperator = HasTemporalExplicitOperator || HasTemporalImplicitOperator; + +/* @class TemporalOperator + * @brief A class to represent an TemporalOperator in NeoFOAMs dsl + * + * The design here is based on the type erasure design pattern + * see https://www.youtube.com/watch?v=4eeESJQk-mw + * + * Motivation for using type erasure is that concrete implementation + * of TemporalOperator e.g Divergence, Laplacian, etc can be stored in a vector of + * TemporalOperator + * + * @ingroup dsl + */ +class TemporalOperator +{ +public: + + + template + TemporalOperator(T cls) : model_(std::make_unique>(std::move(cls))) + {} + + TemporalOperator(const TemporalOperator& eqnOperator); + + TemporalOperator(TemporalOperator&& eqnOperator); + + void explicitOperation(Field& source, scalar t, scalar dt); + + void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt); + + la::LinearSystem createEmptyLinearSystem() const; + + /* returns the fundamental type of an operator, ie explicit, implicit */ + SpatialOperator::Type getType() const; + + std::string getName() const; + + Coeff& getCoefficient(); + + Coeff getCoefficient() const; + + /* @brief Given an input this function reads required coeffs */ + void build(const Input& input); + + /* @brief Get the executor */ + const Executor& exec() const; + + +private: + + /* @brief Base class defining the concept of a term. This effectively + * defines what functions need to be implemented by a concrete Operator implementation + * */ + struct TemporalOperatorConcept + { + virtual ~TemporalOperatorConcept() = default; + + virtual void explicitOperation(Field& source, scalar t, scalar dt) = 0; + + virtual void + implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) = 0; + + virtual la::LinearSystem createEmptyLinearSystem() const = 0; + + /* @brief Given an input this function reads required coeffs */ + virtual void build(const Input& input) = 0; + + /* returns the name of the operator */ + virtual std::string getName() const = 0; + + /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ + virtual SpatialOperator::Type getType() const = 0; + + /* @brief get the associated coefficient for this term */ + virtual Coeff& getCoefficient() = 0; + + /* @brief get the associated coefficient for this term */ + virtual Coeff getCoefficient() const = 0; + + /* @brief Get the executor */ + virtual const Executor& exec() const = 0; + + // The Prototype Design Pattern + virtual std::unique_ptr clone() const = 0; + }; + + // Templated derived class to implement the type-specific behavior + template + struct TemporalOperatorModel : TemporalOperatorConcept + { + /* @brief build with concrete TemporalOperator */ + TemporalOperatorModel(ConcreteTemporalOperatorType concreteOp) + : concreteOp_(std::move(concreteOp)) + {} + + /* returns the name of the operator */ + std::string getName() const override { return concreteOp_.getName(); } + + virtual void explicitOperation(Field& source, scalar t, scalar dt) override + { + if constexpr (HasTemporalExplicitOperator) + { + concreteOp_.explicitOperation(source, t, dt); + } + } + + virtual void + implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) override + { + if constexpr (HasTemporalImplicitOperator) + { + concreteOp_.implicitOperation(ls, t, dt); + } + } + + virtual la::LinearSystem createEmptyLinearSystem() const override + { + // if constexpr (HasTemporalImplicitOperator) + // { + return concreteOp_.createEmptyLinearSystem(); + // } + throw std::runtime_error("Implicit operation not implemented"); + // only need to avoid compiler warning about missing return statement + // this code path should never be reached as we call implicitOperation on an explicit + // operator + NeoFOAM::Field values(exec(), 1, 0.0); + NeoFOAM::Field colIdx(exec(), 1, 0); + NeoFOAM::Field rowPtrs(exec(), 2, 0); + NeoFOAM::la::CSRMatrix csrMatrix( + values, colIdx, rowPtrs + ); + + NeoFOAM::Field rhs(exec(), 1, 0.0); + return la::LinearSystem(csrMatrix, rhs, "custom"); + } + + /* @brief Given an input this function reads required coeffs */ + virtual void build(const Input& input) override { concreteOp_.build(input); } + + /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ + SpatialOperator::Type getType() const override { return concreteOp_.getType(); } + + /* @brief Get the executor */ + const Executor& exec() const override { return concreteOp_.exec(); } + + /* @brief get the associated coefficient for this term */ + virtual Coeff& getCoefficient() override { return concreteOp_.getCoefficient(); } + + /* @brief get the associated coefficient for this term */ + virtual Coeff getCoefficient() const override { return concreteOp_.getCoefficient(); } + + // The Prototype Design Pattern + std::unique_ptr clone() const override + { + return std::make_unique(*this); + } + + ConcreteTemporalOperatorType concreteOp_; + }; + + std::unique_ptr model_; +}; + + +TemporalOperator operator*(scalar scalarCoeff, TemporalOperator rhs); + +TemporalOperator operator*(const Field& coeffField, TemporalOperator rhs); + +TemporalOperator operator*(const Coeff& coeff, TemporalOperator rhs); + +template + requires std::invocable +TemporalOperator operator*([[maybe_unused]] CoeffFunction coeffFunc, const TemporalOperator& lhs) +{ + // TODO implement + NF_ERROR_EXIT("Not implemented"); + TemporalOperator result = lhs; + // if (!result.getCoefficient().useSpan) + // { + // result.setField(std::make_shared>(result.exec(), result.nCells(), 1.0)); + // } + // map(result.exec(), result.getCoefficient().values, scaleFunc); + return result; +} + + +} // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index a835402a0..5c8449867 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -6,7 +6,7 @@ #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/core/executor/executor.hpp" #include "NeoFOAM/core/input.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" @@ -20,14 +20,15 @@ class DdtOperator : public dsl::OperatorMixin> public: - DdtOperator(dsl::Operator::Type termType, VolumeField& field) + DdtOperator(dsl::SpatialOperator::Type termType, VolumeField& field) : dsl::OperatorMixin>(field.exec(), field, termType), sparsityPattern_(SparsityPattern::readOrCreate(field.mesh())) { }; - void explicitOperation(Field& source) + void explicitOperation(Field& source, scalar t, scalar dt) { + const scalar rDeltat = 1 / dt; auto operatorScaling = getCoefficient(); auto [sourceSpan, field, oldField] = spans(source, field_.internalField(), oldTime(field_).internalField()); @@ -35,7 +36,7 @@ class DdtOperator : public dsl::OperatorMixin> source.exec(), source.range(), KOKKOS_LAMBDA(const size_t celli) { - sourceSpan[celli] += field[celli] - oldField[celli]; + sourceSpan[celli] += rDeltat * (field[celli] - oldField[celli]); } ); } @@ -45,9 +46,9 @@ class DdtOperator : public dsl::OperatorMixin> return sparsityPattern_->linearSystem(); } - void implicitOperation(la::LinearSystem& ls) + void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) { - const scalar rDeltat = 1 / 1.0; + const scalar rDeltat = 1 / dt; const auto operatorScaling = getCoefficient(); const auto [diagOffs, oldField] = spans(sparsityPattern_->diagOffset(), oldTime(field_).internalField()); diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index 5b489cbe4..da5f17354 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -7,7 +7,7 @@ #include "NeoFOAM/linearAlgebra/linearSystem.hpp" #include "NeoFOAM/core/executor/executor.hpp" #include "NeoFOAM/core/input.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp" @@ -83,7 +83,7 @@ class DivOperator : public dsl::OperatorMixin> ) {}; DivOperator( - dsl::Operator::Type termType, + dsl::SpatialOperator::Type termType, const SurfaceField& faceFlux, VolumeField& phi, Input input @@ -92,7 +92,7 @@ class DivOperator : public dsl::OperatorMixin> divOperatorStrategy_(DivOperatorFactory::create(exec_, phi.mesh(), input)) {}; DivOperator( - dsl::Operator::Type termType, + dsl::SpatialOperator::Type termType, const SurfaceField& faceFlux, VolumeField& phi, std::unique_ptr divOperatorStrategy @@ -101,7 +101,9 @@ class DivOperator : public dsl::OperatorMixin> divOperatorStrategy_(std::move(divOperatorStrategy)) {}; DivOperator( - dsl::Operator::Type termType, const SurfaceField& faceFlux, VolumeField& phi + dsl::SpatialOperator::Type termType, + const SurfaceField& faceFlux, + VolumeField& phi ) : dsl::OperatorMixin>(phi.exec(), phi, termType), faceFlux_(faceFlux), divOperatorStrategy_(nullptr) {}; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp index cead24c27..d524d68ed 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp @@ -6,7 +6,7 @@ #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/core/executor/executor.hpp" #include "NeoFOAM/core/input.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" @@ -21,7 +21,9 @@ class SourceTerm : public dsl::OperatorMixin> public: SourceTerm( - dsl::Operator::Type termType, VolumeField& coefficients, VolumeField& field + dsl::SpatialOperator::Type termType, + VolumeField& coefficients, + VolumeField& field ) : dsl::OperatorMixin>(field.exec(), field, termType), coefficients_(coefficients), diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cd31ab41f..97d3c2027 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,7 +23,8 @@ target_sources( "core/demangle.cpp" "core/tokenList.cpp" "dsl/coeff.cpp" - "dsl/operator.cpp" + "dsl/spatialOperator.cpp" + "dsl/temporalOperator.cpp" "executor/CPUExecutor.cpp" "executor/GPUExecutor.cpp" "executor/serialExecutor.cpp" diff --git a/src/dsl/operator.cpp b/src/dsl/operator.cpp deleted file mode 100644 index b85c018d9..000000000 --- a/src/dsl/operator.cpp +++ /dev/null @@ -1,60 +0,0 @@ -// SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors - -#include "NeoFOAM/dsl/operator.hpp" - -namespace NeoFOAM::dsl -{ - - -Operator::Operator(const Operator& eqnOperator) : model_ {eqnOperator.model_->clone()} {} - -Operator::Operator(Operator&& eqnOperator) : model_ {std::move(eqnOperator.model_)} {} - -void Operator::explicitOperation(Field& source) { model_->explicitOperation(source); } - -void Operator::temporalOperation(Field& field) { model_->temporalOperation(field); } - -void Operator::implicitOperation(la::LinearSystem& ls) -{ - model_->implicitOperation(ls); -} - -la::LinearSystem Operator::createEmptyLinearSystem() const -{ - return model_->createEmptyLinearSystem(); -} - -Operator::Type Operator::getType() const { return model_->getType(); } - -Coeff& Operator::getCoefficient() { return model_->getCoefficient(); } - -Coeff Operator::getCoefficient() const { return model_->getCoefficient(); } - -void Operator::build(const Input& input) { model_->build(input); } - -const Executor& Operator::exec() const { return model_->exec(); } - -Operator operator*(scalar scalarCoeff, Operator rhs) -{ - Operator result = rhs; - result.getCoefficient() *= scalarCoeff; - return result; -} - -Operator operator*(const Field& coeffField, Operator rhs) -{ - Operator result = rhs; - result.getCoefficient() *= Coeff(coeffField); - return result; -} - -Operator operator*(const Coeff& coeff, Operator rhs) -{ - Operator result = rhs; - result.getCoefficient() *= coeff; - return result; -} - - -} // namespace NeoFOAM::dsl diff --git a/src/dsl/spatialOperator.cpp b/src/dsl/spatialOperator.cpp new file mode 100644 index 000000000..ae1a97693 --- /dev/null +++ b/src/dsl/spatialOperator.cpp @@ -0,0 +1,67 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors + +#include "NeoFOAM/dsl/spatialOperator.hpp" + +namespace NeoFOAM::dsl +{ + + +SpatialOperator::SpatialOperator(const SpatialOperator& eqnOperator) + : model_ {eqnOperator.model_->clone()} +{} + +SpatialOperator::SpatialOperator(SpatialOperator&& eqnOperator) + : model_ {std::move(eqnOperator.model_)} +{} + +void SpatialOperator::explicitOperation(Field& source) +{ + model_->explicitOperation(source); +} + +void SpatialOperator::implicitOperation(la::LinearSystem& ls) +{ + model_->implicitOperation(ls); +} + +la::LinearSystem SpatialOperator::createEmptyLinearSystem() const +{ + return model_->createEmptyLinearSystem(); +} + +SpatialOperator::Type SpatialOperator::getType() const { return model_->getType(); } + +std::string SpatialOperator::getName() const { return model_->getName(); } + +Coeff& SpatialOperator::getCoefficient() { return model_->getCoefficient(); } + +Coeff SpatialOperator::getCoefficient() const { return model_->getCoefficient(); } + +void SpatialOperator::build(const Input& input) { model_->build(input); } + +const Executor& SpatialOperator::exec() const { return model_->exec(); } + +SpatialOperator operator*(scalar scalarCoeff, SpatialOperator rhs) +{ + SpatialOperator result = rhs; + result.getCoefficient() *= scalarCoeff; + return result; +} + +SpatialOperator operator*(const Field& coeffField, SpatialOperator rhs) +{ + SpatialOperator result = rhs; + result.getCoefficient() *= Coeff(coeffField); + return result; +} + +SpatialOperator operator*(const Coeff& coeff, SpatialOperator rhs) +{ + SpatialOperator result = rhs; + result.getCoefficient() *= coeff; + return result; +} + + +} // namespace NeoFOAM::dsl diff --git a/src/dsl/temporalOperator.cpp b/src/dsl/temporalOperator.cpp new file mode 100644 index 000000000..43b120d24 --- /dev/null +++ b/src/dsl/temporalOperator.cpp @@ -0,0 +1,70 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors + +#include "NeoFOAM/dsl/temporalOperator.hpp" + +namespace NeoFOAM::dsl +{ + + +TemporalOperator::TemporalOperator(const TemporalOperator& eqnOperator) + : model_ {eqnOperator.model_->clone()} +{} + +TemporalOperator::TemporalOperator(TemporalOperator&& eqnOperator) + : model_ {std::move(eqnOperator.model_)} +{} + +void TemporalOperator::explicitOperation(Field& source, scalar t, scalar dt) +{ + model_->explicitOperation(source, t, dt); +} + + +void TemporalOperator::implicitOperation( + la::LinearSystem& ls, scalar t, scalar dt +) +{ + model_->implicitOperation(ls, t, dt); +} + +la::LinearSystem TemporalOperator::createEmptyLinearSystem() const +{ + return model_->createEmptyLinearSystem(); +} + +SpatialOperator::Type TemporalOperator::getType() const { return model_->getType(); } + +std::string TemporalOperator::getName() const { return model_->getName(); } + +Coeff& TemporalOperator::getCoefficient() { return model_->getCoefficient(); } + +Coeff TemporalOperator::getCoefficient() const { return model_->getCoefficient(); } + +void TemporalOperator::build(const Input& input) { model_->build(input); } + +const Executor& TemporalOperator::exec() const { return model_->exec(); } + +TemporalOperator operator*(scalar scalarCoeff, TemporalOperator rhs) +{ + TemporalOperator result = rhs; + result.getCoefficient() *= scalarCoeff; + return result; +} + +TemporalOperator operator*(const Field& coeffField, TemporalOperator rhs) +{ + TemporalOperator result = rhs; + result.getCoefficient() *= Coeff(coeffField); + return result; +} + +TemporalOperator operator*(const Coeff& coeff, TemporalOperator rhs) +{ + TemporalOperator result = rhs; + result.getCoefficient() *= coeff; + return result; +} + + +} // namespace NeoFOAM::dsl diff --git a/test/dsl/CMakeLists.txt b/test/dsl/CMakeLists.txt index 48ace3536..75abe98c4 100644 --- a/test/dsl/CMakeLists.txt +++ b/test/dsl/CMakeLists.txt @@ -3,4 +3,5 @@ neofoam_unit_test(coeff) neofoam_unit_test(expression) -neofoam_unit_test(operator) +neofoam_unit_test(spatialOperator) +neofoam_unit_test(temporalOperator) diff --git a/test/dsl/common.hpp b/test/dsl/common.hpp index 7fc3e0ba1..0f4b8c30f 100644 --- a/test/dsl/common.hpp +++ b/test/dsl/common.hpp @@ -13,22 +13,24 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; using Field = NeoFOAM::Field; using Coeff = NeoFOAM::dsl::Coeff; -using Operator = NeoFOAM::dsl::Operator; +using SpatialOperator = NeoFOAM::dsl::SpatialOperator; using Executor = NeoFOAM::Executor; using VolumeField = fvcc::VolumeField; using OperatorMixin = NeoFOAM::dsl::OperatorMixin; using BoundaryFields = NeoFOAM::BoundaryFields; -/* A dummy implementation of a Operator - * following the Operator interface */ +/* A dummy implementation of a SpatialOperator + * following the SpatialOperator interface */ class Dummy : public OperatorMixin { public: - Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} + Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) + {} - Dummy(VolumeField& field, Operator::Type type) : OperatorMixin(field.exec(), field, type) {} + Dummy(VolumeField& field, SpatialOperator::Type type) : OperatorMixin(field.exec(), field, type) + {} void explicitOperation(Field& source) { @@ -83,6 +85,78 @@ class Dummy : public OperatorMixin std::string getName() const { return "Dummy"; } }; +/* A dummy implementation of a SpatialOperator + * following the SpatialOperator interface */ +class TemporalDummy : public OperatorMixin +{ + +public: + + TemporalDummy(VolumeField& field) + : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) + {} + + TemporalDummy(VolumeField& field, SpatialOperator::Type type) + : OperatorMixin(field.exec(), field, type) + {} + + void explicitOperation(Field& source, NeoFOAM::scalar t, NeoFOAM::scalar dt) + { + auto sourceSpan = source.span(); + auto fieldSpan = field_.internalField().span(); + auto coeff = getCoefficient(); + NeoFOAM::parallelFor( + source.exec(), + source.range(), + KOKKOS_LAMBDA(const size_t i) { sourceSpan[i] += coeff[i] * fieldSpan[i]; } + ); + } + + void implicitOperation( + la::LinearSystem& ls, + NeoFOAM::scalar t, + NeoFOAM::scalar dt + ) + { + auto values = ls.matrix().values(); + auto rhs = ls.rhs().span(); + auto fieldSpan = field_.internalField().span(); + auto coeff = getCoefficient(); + + // update diag + NeoFOAM::parallelFor( + exec(), + {0, values.size()}, + KOKKOS_LAMBDA(const size_t i) { values[i] += coeff[i] * fieldSpan[i]; } + ); + + // update rhs + NeoFOAM::parallelFor( + exec(), + ls.rhs().range(), + KOKKOS_LAMBDA(const size_t i) { rhs[i] += coeff[i] * fieldSpan[i]; } + ); + } + + la::LinearSystem createEmptyLinearSystem() const + { + NeoFOAM::Field values(exec(), 1, 0.0); + NeoFOAM::Field colIdx(exec(), 1, 0.0); + NeoFOAM::Field rowPtrs(exec(), {0, 1}); + NeoFOAM::la::CSRMatrix csrMatrix( + values, colIdx, rowPtrs + ); + + NeoFOAM::Field rhs(exec(), 1, 0.0); + NeoFOAM::la::LinearSystem linearSystem( + csrMatrix, rhs, "diagonal" + ); + return linearSystem; + } + + std::string getName() const { return "TemporalDummy"; } +}; + NeoFOAM::scalar getField(const NeoFOAM::Field& source) { auto sourceField = source.copyToHost(); diff --git a/test/dsl/expression.cpp b/test/dsl/expression.cpp index d7a4a913f..63cfeead7 100644 --- a/test/dsl/expression.cpp +++ b/test/dsl/expression.cpp @@ -53,12 +53,12 @@ TEST_CASE("Expression") SECTION("Create equation and perform implicitOperation on " + execName) { - auto a = Dummy(vf, Operator::Type::Implicit); - auto b = Dummy(vf, Operator::Type::Implicit); + auto a = Dummy(vf, SpatialOperator::Type::Implicit); + auto b = Dummy(vf, SpatialOperator::Type::Implicit); auto eqnA = a + b; - auto eqnB = - fB * Dummy(vf, Operator::Type::Implicit) + 2 * Dummy(vf, Operator::Type::Implicit); + auto eqnB = fB * Dummy(vf, SpatialOperator::Type::Implicit) + + 2 * Dummy(vf, SpatialOperator::Type::Implicit); auto eqnC = Expression(2 * a - b); auto eqnD = 3 * (2 * a - b); auto eqnE = (2 * a - b) + (2 * a - b); diff --git a/test/dsl/operator.cpp b/test/dsl/spatialOperator.cpp similarity index 80% rename from test/dsl/operator.cpp rename to test/dsl/spatialOperator.cpp index 1d2f1ed31..8f17f6bf6 100644 --- a/test/dsl/operator.cpp +++ b/test/dsl/spatialOperator.cpp @@ -4,7 +4,9 @@ // a custom main #include "common.hpp" -TEST_CASE("Operator") +namespace dsl = NeoFOAM::dsl; + +TEST_CASE("SpatialOperator") { NeoFOAM::Executor exec = GENERATE( NeoFOAM::Executor(NeoFOAM::SerialExecutor {}), @@ -16,17 +18,17 @@ TEST_CASE("Operator") auto mesh = NeoFOAM::createSingleCellMesh(exec); - SECTION("Operator creation on " + execName) + SECTION("SpatialOperator creation on " + execName) { Field fA(exec, 1, 2.0); BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); std::vector> bcs {}; auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - auto b = Dummy(vf); + dsl::SpatialOperator b = Dummy(vf); REQUIRE(b.getName() == "Dummy"); - REQUIRE(b.getType() == Operator::Type::Explicit); + REQUIRE(b.getType() == SpatialOperator::Type::Explicit); } SECTION("Supports Coefficients Explicit " + execName) @@ -38,9 +40,9 @@ TEST_CASE("Operator") BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - auto c = 2 * Dummy(vf); - auto d = fB * Dummy(vf); - auto e = Coeff(-3, fB) * Dummy(vf); + dsl::SpatialOperator c = 2.0 * dsl::SpatialOperator(Dummy(vf)); + dsl::SpatialOperator d = fB * dsl::SpatialOperator(Dummy(vf)); + dsl::SpatialOperator e = Coeff(-3, fB) * dsl::SpatialOperator(Dummy(vf)); [[maybe_unused]] auto coeffC = c.getCoefficient(); [[maybe_unused]] auto coeffD = d.getCoefficient(); @@ -71,10 +73,10 @@ TEST_CASE("Operator") std::vector> bcs {}; auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - auto b = Dummy(vf, Operator::Type::Implicit); + dsl::SpatialOperator b = Dummy(vf, SpatialOperator::Type::Implicit); REQUIRE(b.getName() == "Dummy"); - REQUIRE(b.getType() == Operator::Type::Implicit); + REQUIRE(b.getType() == SpatialOperator::Type::Implicit); auto ls = b.createEmptyLinearSystem(); REQUIRE(ls.matrix().nValues() == 1); @@ -91,9 +93,12 @@ TEST_CASE("Operator") BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - auto c = 2 * Dummy(vf, Operator::Type::Implicit); - auto d = fB * Dummy(vf, Operator::Type::Implicit); - auto e = Coeff(-3, fB) * Dummy(vf, Operator::Type::Implicit); + dsl::SpatialOperator c = + 2 * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); + dsl::SpatialOperator d = + fB * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); + dsl::SpatialOperator e = + Coeff(-3, fB) * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); [[maybe_unused]] auto coeffC = c.getCoefficient(); [[maybe_unused]] auto coeffD = d.getCoefficient(); diff --git a/test/dsl/temporalOperator.cpp b/test/dsl/temporalOperator.cpp new file mode 100644 index 000000000..a6d3b865e --- /dev/null +++ b/test/dsl/temporalOperator.cpp @@ -0,0 +1,137 @@ +// 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 "common.hpp" + +namespace dsl = NeoFOAM::dsl; + +TEST_CASE("TemporalOperator") +{ + 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::createSingleCellMesh(exec); + + SECTION("Operator creation on " + execName) + { + Field fA(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + + std::vector> bcs {}; + auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); + dsl::TemporalOperator b = TemporalDummy(vf); + + REQUIRE(b.getName() == "TemporalDummy"); + REQUIRE(b.getType() == dsl::SpatialOperator::Type::Explicit); + } + + SECTION("Supports Coefficients Explicit " + execName) + { + std::vector> bcs {}; + NeoFOAM::scalar t = 0.0; + NeoFOAM::scalar dt = 0.1; + + Field fA(exec, 1, 2.0); + Field fB(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); + + dsl::TemporalOperator c = 2 * dsl::TemporalOperator(TemporalDummy(vf)); + dsl::TemporalOperator d = fB * dsl::TemporalOperator(TemporalDummy(vf)); + dsl::TemporalOperator e = Coeff(-3, fB) * dsl::TemporalOperator(TemporalDummy(vf)); + + [[maybe_unused]] auto coeffC = c.getCoefficient(); + [[maybe_unused]] auto coeffD = d.getCoefficient(); + [[maybe_unused]] auto coeffE = e.getCoefficient(); + + Field source(exec, 1, 2.0); + c.explicitOperation(source, t, dt); + + // 2 += 2 * 2 + auto hostSourceC = source.copyToHost(); + REQUIRE(hostSourceC.span()[0] == 6.0); + + // 6 += 2 * 2 + d.explicitOperation(source, t, dt); + auto hostSourceD = source.copyToHost(); + REQUIRE(hostSourceD.span()[0] == 10.0); + + // 10 += - 6 * 2 + e.explicitOperation(source, t, dt); + auto hostSourceE = source.copyToHost(); + REQUIRE(hostSourceE.span()[0] == -2.0); + } + + SECTION("Implicit Operations " + execName) + { + Field fA(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + + std::vector> bcs {}; + auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); + dsl::TemporalOperator b = TemporalDummy(vf, SpatialOperator::Type::Implicit); + + REQUIRE(b.getName() == "TemporalDummy"); + REQUIRE(b.getType() == SpatialOperator::Type::Implicit); + + auto ls = b.createEmptyLinearSystem(); + REQUIRE(ls.matrix().nValues() == 1); + REQUIRE(ls.matrix().nColIdxs() == 1); + REQUIRE(ls.matrix().nRows() == 1); + } + + SECTION("Supports Coefficients Implicit " + execName) + { + std::vector> bcs {}; + NeoFOAM::scalar t = 0.0; + NeoFOAM::scalar dt = 0.1; + + Field fA(exec, 1, 2.0); + Field fB(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); + + auto c = 2 * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); + auto d = fB * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); + auto e = Coeff(-3, fB) + * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); + + [[maybe_unused]] auto coeffC = c.getCoefficient(); + [[maybe_unused]] auto coeffD = d.getCoefficient(); + [[maybe_unused]] auto coeffE = e.getCoefficient(); + + // Field source(exec, 1, 2.0); + auto ls = c.createEmptyLinearSystem(); + c.implicitOperation(ls, t, dt); + + // c = 2 * 2 + auto hostRhsC = ls.rhs().copyToHost(); + REQUIRE(hostRhsC.span()[0] == 4.0); + auto hostLsC = ls.copyToHost(); + REQUIRE(hostLsC.matrix().values()[0] == 4.0); + + + // d= 2 * 2 + ls = d.createEmptyLinearSystem(); + d.implicitOperation(ls, t, dt); + auto hostRhsD = ls.rhs().copyToHost(); + REQUIRE(hostRhsD.span()[0] == 4.0); + auto hostLsD = ls.copyToHost(); + REQUIRE(hostLsD.matrix().values()[0] == 4.0); + + + // e = - -3 * 2 * 2 = -12 + ls = e.createEmptyLinearSystem(); + e.implicitOperation(ls, t, dt); + auto hostRhsE = ls.rhs().copyToHost(); + REQUIRE(hostRhsE.span()[0] == -12.0); + auto hostLsE = ls.copyToHost(); + REQUIRE(hostLsE.matrix().values()[0] == -12.0); + } +} diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index 241acd149..bfda71e11 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -11,7 +11,7 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; -using Operator = NeoFOAM::dsl::Operator; +using SpatialOperator = NeoFOAM::dsl::SpatialOperator; TEST_CASE("DivOperator") { @@ -36,7 +36,7 @@ TEST_CASE("DivOperator") SECTION("Construct from Token" + execName) { NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); + fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); } SECTION("Construct from Dictionary" + execName) @@ -45,6 +45,6 @@ TEST_CASE("DivOperator") {{std::string("DivOperator"), std::string("Gauss")}, {std::string("surfaceInterpolation"), std::string("linear")}} ); - fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); + fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); } } diff --git a/test/timeIntegration/implicitTimeIntegration.cpp b/test/timeIntegration/implicitTimeIntegration.cpp index a9de82ef9..500fed4b7 100644 --- a/test/timeIntegration/implicitTimeIntegration.cpp +++ b/test/timeIntegration/implicitTimeIntegration.cpp @@ -68,7 +68,7 @@ TEST_CASE("TimeIntegration") SECTION("Create expression and perform explicitOperation on " + execName) { auto dummy = Dummy(vf); - Operator ddtOperator = NeoFOAM::dsl::temporal::ddt(vf); + NeoFOAM::dsl::TemporalOperator ddtOperator = NeoFOAM::dsl::imp::ddt(vf); // ddt(U) = f auto eqn = ddtOperator + dummy; diff --git a/test/timeIntegration/rungeKutta.cpp b/test/timeIntegration/rungeKutta.cpp index 5267153ef..3db667d43 100644 --- a/test/timeIntegration/rungeKutta.cpp +++ b/test/timeIntegration/rungeKutta.cpp @@ -17,7 +17,8 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; using Field = NeoFOAM::Field; using Coeff = NeoFOAM::dsl::Coeff; -using Operator = NeoFOAM::dsl::Operator; +using SpatialOperator = NeoFOAM::dsl::SpatialOperator; +using TemporalOperator = NeoFOAM::dsl::TemporalOperator; using Executor = NeoFOAM::Executor; using VolumeField = fvcc::VolumeField; using OperatorMixin = NeoFOAM::dsl::OperatorMixin; @@ -32,7 +33,9 @@ class YSquared : public OperatorMixin public: - YSquared(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} + YSquared(VolumeField& field) + : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) + {} void explicitOperation(Field& source) const { @@ -119,7 +122,8 @@ TEST_CASE("TimeIntegration - Runge Kutta") vfOld.internalField() = initialValue; // Set expression - Operator ddtOp = Ddt(vfOld); + TemporalOperator ddtOp = NeoFOAM::dsl::imp::ddt(vfOld); + auto divOp = YSquared(vfOld); auto eqn = ddtOp + divOp; diff --git a/test/timeIntegration/timeIntegration.cpp b/test/timeIntegration/timeIntegration.cpp index a9de82ef9..500fed4b7 100644 --- a/test/timeIntegration/timeIntegration.cpp +++ b/test/timeIntegration/timeIntegration.cpp @@ -68,7 +68,7 @@ TEST_CASE("TimeIntegration") SECTION("Create expression and perform explicitOperation on " + execName) { auto dummy = Dummy(vf); - Operator ddtOperator = NeoFOAM::dsl::temporal::ddt(vf); + NeoFOAM::dsl::TemporalOperator ddtOperator = NeoFOAM::dsl::imp::ddt(vf); // ddt(U) = f auto eqn = ddtOperator + dummy; From 0603f4be2fce7f4a76ec5599cae073b6abc3bd78 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 9 Feb 2025 09:52:12 +0100 Subject: [PATCH 10/21] Fix: divOperator test still used Operator --- benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp b/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp index 3d0d2e522..3290af481 100644 --- a/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp @@ -34,7 +34,7 @@ TEST_CASE("DivOperator::div", "[bench]") DYNAMIC_SECTION("" << size) { NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - auto op = fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); + auto op = fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); BENCHMARK(std::string(execName)) { return (op.div(divPhi)); }; } From e26c3b830a8913be34a5cd3eee41b22a5394bc56 Mon Sep 17 00:00:00 2001 From: Henning Scheufler <16381681+HenningScheufler@users.noreply.github.com> Date: Sun, 9 Feb 2025 21:14:33 +0100 Subject: [PATCH 11/21] Apply suggestions from code review Co-authored-by: Gregor Olenik --- include/NeoFOAM/dsl/ddt.hpp | 8 ++++---- include/NeoFOAM/dsl/solver.hpp | 4 +++- .../cellCentred/operators/fvccLinearSystem.hpp | 5 +---- src/dsl/spatialOperator.cpp | 2 +- src/dsl/temporalOperator.cpp | 2 +- 5 files changed, 10 insertions(+), 11 deletions(-) diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp index 583f985c5..0270ad884 100644 --- a/include/NeoFOAM/dsl/ddt.hpp +++ b/include/NeoFOAM/dsl/ddt.hpp @@ -23,15 +23,15 @@ class Ddt : public OperatorMixin std::string getName() const { return "TimeOperator"; } void - explicitOperation([[maybe_unused]] Field& source, NeoFOAM::scalar t, NeoFOAM::scalar dt) + explicitOperation([[maybe_unused]] Field& source, [[maybe_unused]] scalar t, [[maybe_unused]] scalar dt) { NF_ERROR_EXIT("Not implemented"); } void implicitOperation( - la::LinearSystem& ls, - NeoFOAM::scalar t, - NeoFOAM::scalar dt + la::LinearSystem& ls, + scalar t, + scalar dt ) { NF_ERROR_EXIT("Not implemented"); diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index 38eaab720..c1716eb50 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -14,7 +14,9 @@ #include "NeoFOAM/dsl/expression.hpp" #include "NeoFOAM/timeIntegration/timeIntegration.hpp" +#if NF_WITH_GINKGO #include "NeoFOAM/linearAlgebra/ginkgo.hpp" +#endif namespace NeoFOAM::dsl @@ -84,7 +86,7 @@ void solve( [[maybe_unused]] const Dictionary& fvSolution ) { - // todo fix + // FIXME: if (exp.temporalOperators().size() == 0 && exp.spatialOperators().size() == 0) { NF_ERROR_EXIT("No temporal or implicit terms to solve."); diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp index 29c8959b6..f06206493 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp @@ -56,10 +56,7 @@ class LinearSystem void diag(Field& field) { - if (field.size() != sparsityPattern_->diagOffset().size()) - { - NF_THROW("fvcc:LinearSystem:diag: field size mismatch"); - } + NF_ASSERT_EQUAL(field.size() , sparsityPattern_->diagOffset().size()); const auto diagOffset = sparsityPattern_->diagOffset().span(); const auto rowPtrs = ls_.matrix().rowPtrs(); std::span fieldSpan = field.span(); diff --git a/src/dsl/spatialOperator.cpp b/src/dsl/spatialOperator.cpp index ae1a97693..488b1477e 100644 --- a/src/dsl/spatialOperator.cpp +++ b/src/dsl/spatialOperator.cpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors #include "NeoFOAM/dsl/spatialOperator.hpp" diff --git a/src/dsl/temporalOperator.cpp b/src/dsl/temporalOperator.cpp index 43b120d24..049311a1b 100644 --- a/src/dsl/temporalOperator.cpp +++ b/src/dsl/temporalOperator.cpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors #include "NeoFOAM/dsl/temporalOperator.hpp" From 4185e23ca400d9f756e1c9a5f683c873c0fe513d Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 9 Feb 2025 21:26:04 +0100 Subject: [PATCH 12/21] Introduced Operator::Type --- .../cellCentred/operator/divOperator.cpp | 2 +- include/NeoFOAM/dsl/ddt.hpp | 16 +++++++--------- include/NeoFOAM/dsl/explicit.hpp | 4 ++-- include/NeoFOAM/dsl/expression.hpp | 6 +++--- include/NeoFOAM/dsl/implicit.hpp | 6 +++--- include/NeoFOAM/dsl/operator.hpp | 19 +++++++++++++++++++ include/NeoFOAM/dsl/spatialOperator.hpp | 13 +++++++------ include/NeoFOAM/dsl/temporalOperator.hpp | 8 ++++---- .../cellCentred/operators/ddtOperator.hpp | 2 +- .../cellCentred/operators/divOperator.hpp | 8 +++----- .../cellCentred/operators/sourceTerm.hpp | 4 +--- src/dsl/spatialOperator.cpp | 2 +- src/dsl/temporalOperator.cpp | 2 +- test/dsl/common.hpp | 13 +++++-------- test/dsl/expression.cpp | 8 ++++---- test/dsl/spatialOperator.cpp | 14 ++++++-------- test/dsl/temporalOperator.cpp | 13 ++++++------- .../cellCentred/operator/divOperator.cpp | 6 +++--- test/timeIntegration/rungeKutta.cpp | 4 +--- test/timeIntegration/timeIntegration.cpp | 2 +- 20 files changed, 79 insertions(+), 73 deletions(-) create mode 100644 include/NeoFOAM/dsl/operator.hpp diff --git a/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp b/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp index 3290af481..3d0d2e522 100644 --- a/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/benchmarks/finiteVolume/cellCentred/operator/divOperator.cpp @@ -34,7 +34,7 @@ TEST_CASE("DivOperator::div", "[bench]") DYNAMIC_SECTION("" << size) { NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - auto op = fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); + auto op = fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); BENCHMARK(std::string(execName)) { return (op.div(divPhi)); }; } diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp index 0270ad884..b98452dee 100644 --- a/include/NeoFOAM/dsl/ddt.hpp +++ b/include/NeoFOAM/dsl/ddt.hpp @@ -16,23 +16,21 @@ class Ddt : public OperatorMixin public: - Ddt(FieldType& field) - : OperatorMixin(field.exec(), field, SpatialOperator::Type::Implicit) + Ddt(FieldType& field) : OperatorMixin(field.exec(), field, Operator::Type::Implicit) {} std::string getName() const { return "TimeOperator"; } - void - explicitOperation([[maybe_unused]] Field& source, [[maybe_unused]] scalar t, [[maybe_unused]] scalar dt) + void explicitOperation( + [[maybe_unused]] Field& source, + [[maybe_unused]] scalar t, + [[maybe_unused]] scalar dt + ) { NF_ERROR_EXIT("Not implemented"); } - void implicitOperation( - la::LinearSystem& ls, - scalar t, - scalar dt - ) + void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) { NF_ERROR_EXIT("Not implemented"); } diff --git a/include/NeoFOAM/dsl/explicit.hpp b/include/NeoFOAM/dsl/explicit.hpp index eb3c423ff..ac33ea23f 100644 --- a/include/NeoFOAM/dsl/explicit.hpp +++ b/include/NeoFOAM/dsl/explicit.hpp @@ -22,13 +22,13 @@ namespace NeoFOAM::dsl::exp SpatialOperator div(const fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return SpatialOperator(fvcc::DivOperator(dsl::SpatialOperator::Type::Explicit, faceFlux, phi)); + return SpatialOperator(fvcc::DivOperator(dsl::Operator::Type::Explicit, faceFlux, phi)); } SpatialOperator Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) { - return SpatialOperator(fvcc::SourceTerm(dsl::SpatialOperator::Type::Explicit, coeff, phi)); + return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Explicit, coeff, phi)); } diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp index d2ff1ad2e..df7a97217 100644 --- a/include/NeoFOAM/dsl/expression.hpp +++ b/include/NeoFOAM/dsl/expression.hpp @@ -54,14 +54,14 @@ class Expression { for (auto& oper : spatialOperators_) { - if (oper.getType() == SpatialOperator::Type::Explicit) + if (oper.getType() == Operator::Type::Explicit) { oper.explicitOperation(source); } } // for (auto& oper : temporalOperators_) // { - // if (oper.getType() == SpatialOperator::Type::Explicit) + // if (oper.getType() == Operator::Type::Explicit) // { // oper.explicitOperation(source); // } @@ -80,7 +80,7 @@ class Expression auto ls = spatialOperators_[0].createEmptyLinearSystem(); for (auto& oper : spatialOperators_) { - if (oper.getType() == SpatialOperator::Type::Implicit) + if (oper.getType() == Operator::Type::Implicit) { oper.implicitOperation(ls); } diff --git a/include/NeoFOAM/dsl/implicit.hpp b/include/NeoFOAM/dsl/implicit.hpp index 06bbcf8bc..1a80ae7e1 100644 --- a/include/NeoFOAM/dsl/implicit.hpp +++ b/include/NeoFOAM/dsl/implicit.hpp @@ -20,19 +20,19 @@ namespace NeoFOAM::dsl::imp TemporalOperator ddt(fvcc::VolumeField& phi) { - return fvcc::DdtOperator(dsl::SpatialOperator::Type::Implicit, phi); + return fvcc::DdtOperator(dsl::Operator::Type::Implicit, phi); } SpatialOperator Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) { - return SpatialOperator(fvcc::SourceTerm(dsl::SpatialOperator::Type::Implicit, coeff, phi)); + return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Implicit, coeff, phi)); } SpatialOperator div(fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return SpatialOperator(fvcc::DivOperator(dsl::SpatialOperator::Type::Implicit, faceFlux, phi)); + return SpatialOperator(fvcc::DivOperator(dsl::Operator::Type::Implicit, faceFlux, phi)); } } // namespace NeoFOAM diff --git a/include/NeoFOAM/dsl/operator.hpp b/include/NeoFOAM/dsl/operator.hpp new file mode 100644 index 000000000..12d8e5af4 --- /dev/null +++ b/include/NeoFOAM/dsl/operator.hpp @@ -0,0 +1,19 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#pragma once + +namespace NeoFOAM::dsl +{ + +class Operator +{ +public: + + enum class Type + { + Implicit, + Explicit + }; +}; + +} // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/spatialOperator.hpp b/include/NeoFOAM/dsl/spatialOperator.hpp index 6a84a6750..a1080c6c3 100644 --- a/include/NeoFOAM/dsl/spatialOperator.hpp +++ b/include/NeoFOAM/dsl/spatialOperator.hpp @@ -10,6 +10,7 @@ #include "NeoFOAM/linearAlgebra/linearSystem.hpp" #include "NeoFOAM/core/input.hpp" #include "NeoFOAM/dsl/coeff.hpp" +#include "NeoFOAM/dsl/operator.hpp" namespace la = NeoFOAM::la; @@ -71,7 +72,7 @@ class SpatialOperator la::LinearSystem createEmptyLinearSystem() const; /* returns the fundamental type of an operator, ie explicit, implicit */ - SpatialOperator::Type getType() const; + Operator::Type getType() const; std::string getName() const; @@ -108,7 +109,7 @@ class SpatialOperator virtual std::string getName() const = 0; /* returns the fundamental type of an operator, ie explicit, implicit */ - virtual SpatialOperator::Type getType() const = 0; + virtual Operator::Type getType() const = 0; /* @brief get the associated coefficient for this term */ virtual Coeff& getCoefficient() = 0; @@ -174,7 +175,7 @@ class SpatialOperator virtual void build(const Input& input) override { concreteOp_.build(input); } /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - SpatialOperator::Type getType() const override { return concreteOp_.getType(); } + Operator::Type getType() const override { return concreteOp_.getType(); } /* @brief Get the executor */ const Executor& exec() const override { return concreteOp_.exec(); } @@ -231,10 +232,10 @@ class OperatorMixin public: - OperatorMixin(const Executor exec, FieldType& field, SpatialOperator::Type type) + OperatorMixin(const Executor exec, FieldType& field, Operator::Type type) : exec_(exec), coeffs_(), field_(field), type_(type) {}; - SpatialOperator::Type getType() const { return type_; } + Operator::Type getType() const { return type_; } virtual ~OperatorMixin() = default; @@ -259,6 +260,6 @@ class OperatorMixin FieldType& field_; - SpatialOperator::Type type_; + Operator::Type type_; }; } // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/temporalOperator.hpp b/include/NeoFOAM/dsl/temporalOperator.hpp index 47c3b604a..4a739daeb 100644 --- a/include/NeoFOAM/dsl/temporalOperator.hpp +++ b/include/NeoFOAM/dsl/temporalOperator.hpp @@ -10,7 +10,7 @@ #include "NeoFOAM/linearAlgebra/linearSystem.hpp" #include "NeoFOAM/core/input.hpp" #include "NeoFOAM/dsl/coeff.hpp" -#include "NeoFOAM/dsl/spatialOperator.hpp" +#include "NeoFOAM/dsl/operator.hpp" namespace la = NeoFOAM::la; @@ -75,7 +75,7 @@ class TemporalOperator la::LinearSystem createEmptyLinearSystem() const; /* returns the fundamental type of an operator, ie explicit, implicit */ - SpatialOperator::Type getType() const; + Operator::Type getType() const; std::string getName() const; @@ -113,7 +113,7 @@ class TemporalOperator virtual std::string getName() const = 0; /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - virtual SpatialOperator::Type getType() const = 0; + virtual Operator::Type getType() const = 0; /* @brief get the associated coefficient for this term */ virtual Coeff& getCoefficient() = 0; @@ -182,7 +182,7 @@ class TemporalOperator virtual void build(const Input& input) override { concreteOp_.build(input); } /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ - SpatialOperator::Type getType() const override { return concreteOp_.getType(); } + Operator::Type getType() const override { return concreteOp_.getType(); } /* @brief Get the executor */ const Executor& exec() const override { return concreteOp_.exec(); } diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index 5c8449867..0b6eb5e39 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -20,7 +20,7 @@ class DdtOperator : public dsl::OperatorMixin> public: - DdtOperator(dsl::SpatialOperator::Type termType, VolumeField& field) + DdtOperator(dsl::Operator::Type termType, VolumeField& field) : dsl::OperatorMixin>(field.exec(), field, termType), sparsityPattern_(SparsityPattern::readOrCreate(field.mesh())) { diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index da5f17354..f6baa9984 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp @@ -83,7 +83,7 @@ class DivOperator : public dsl::OperatorMixin> ) {}; DivOperator( - dsl::SpatialOperator::Type termType, + dsl::Operator::Type termType, const SurfaceField& faceFlux, VolumeField& phi, Input input @@ -92,7 +92,7 @@ class DivOperator : public dsl::OperatorMixin> divOperatorStrategy_(DivOperatorFactory::create(exec_, phi.mesh(), input)) {}; DivOperator( - dsl::SpatialOperator::Type termType, + dsl::Operator::Type termType, const SurfaceField& faceFlux, VolumeField& phi, std::unique_ptr divOperatorStrategy @@ -101,9 +101,7 @@ class DivOperator : public dsl::OperatorMixin> divOperatorStrategy_(std::move(divOperatorStrategy)) {}; DivOperator( - dsl::SpatialOperator::Type termType, - const SurfaceField& faceFlux, - VolumeField& phi + dsl::Operator::Type termType, const SurfaceField& faceFlux, VolumeField& phi ) : dsl::OperatorMixin>(phi.exec(), phi, termType), faceFlux_(faceFlux), divOperatorStrategy_(nullptr) {}; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp index d524d68ed..d1818b8f5 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp @@ -21,9 +21,7 @@ class SourceTerm : public dsl::OperatorMixin> public: SourceTerm( - dsl::SpatialOperator::Type termType, - VolumeField& coefficients, - VolumeField& field + dsl::Operator::Type termType, VolumeField& coefficients, VolumeField& field ) : dsl::OperatorMixin>(field.exec(), field, termType), coefficients_(coefficients), diff --git a/src/dsl/spatialOperator.cpp b/src/dsl/spatialOperator.cpp index 488b1477e..69fc1ed82 100644 --- a/src/dsl/spatialOperator.cpp +++ b/src/dsl/spatialOperator.cpp @@ -30,7 +30,7 @@ la::LinearSystem SpatialOperator::createEmptyLinearSystem() co return model_->createEmptyLinearSystem(); } -SpatialOperator::Type SpatialOperator::getType() const { return model_->getType(); } +Operator::Type SpatialOperator::getType() const { return model_->getType(); } std::string SpatialOperator::getName() const { return model_->getName(); } diff --git a/src/dsl/temporalOperator.cpp b/src/dsl/temporalOperator.cpp index 049311a1b..0bba489c2 100644 --- a/src/dsl/temporalOperator.cpp +++ b/src/dsl/temporalOperator.cpp @@ -33,7 +33,7 @@ la::LinearSystem TemporalOperator::createEmptyLinearSystem() c return model_->createEmptyLinearSystem(); } -SpatialOperator::Type TemporalOperator::getType() const { return model_->getType(); } +Operator::Type TemporalOperator::getType() const { return model_->getType(); } std::string TemporalOperator::getName() const { return model_->getName(); } diff --git a/test/dsl/common.hpp b/test/dsl/common.hpp index 0f4b8c30f..5f14ffdaf 100644 --- a/test/dsl/common.hpp +++ b/test/dsl/common.hpp @@ -13,7 +13,7 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; using Field = NeoFOAM::Field; using Coeff = NeoFOAM::dsl::Coeff; -using SpatialOperator = NeoFOAM::dsl::SpatialOperator; +using Operator = NeoFOAM::dsl::Operator; using Executor = NeoFOAM::Executor; using VolumeField = fvcc::VolumeField; using OperatorMixin = NeoFOAM::dsl::OperatorMixin; @@ -26,11 +26,9 @@ class Dummy : public OperatorMixin public: - Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) - {} + Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} - Dummy(VolumeField& field, SpatialOperator::Type type) : OperatorMixin(field.exec(), field, type) - {} + Dummy(VolumeField& field, Operator::Type type) : OperatorMixin(field.exec(), field, type) {} void explicitOperation(Field& source) { @@ -92,11 +90,10 @@ class TemporalDummy : public OperatorMixin public: - TemporalDummy(VolumeField& field) - : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) + TemporalDummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} - TemporalDummy(VolumeField& field, SpatialOperator::Type type) + TemporalDummy(VolumeField& field, Operator::Type type) : OperatorMixin(field.exec(), field, type) {} diff --git a/test/dsl/expression.cpp b/test/dsl/expression.cpp index 63cfeead7..d7a4a913f 100644 --- a/test/dsl/expression.cpp +++ b/test/dsl/expression.cpp @@ -53,12 +53,12 @@ TEST_CASE("Expression") SECTION("Create equation and perform implicitOperation on " + execName) { - auto a = Dummy(vf, SpatialOperator::Type::Implicit); - auto b = Dummy(vf, SpatialOperator::Type::Implicit); + auto a = Dummy(vf, Operator::Type::Implicit); + auto b = Dummy(vf, Operator::Type::Implicit); auto eqnA = a + b; - auto eqnB = fB * Dummy(vf, SpatialOperator::Type::Implicit) - + 2 * Dummy(vf, SpatialOperator::Type::Implicit); + auto eqnB = + fB * Dummy(vf, Operator::Type::Implicit) + 2 * Dummy(vf, Operator::Type::Implicit); auto eqnC = Expression(2 * a - b); auto eqnD = 3 * (2 * a - b); auto eqnE = (2 * a - b) + (2 * a - b); diff --git a/test/dsl/spatialOperator.cpp b/test/dsl/spatialOperator.cpp index 8f17f6bf6..ec15b597d 100644 --- a/test/dsl/spatialOperator.cpp +++ b/test/dsl/spatialOperator.cpp @@ -28,7 +28,7 @@ TEST_CASE("SpatialOperator") dsl::SpatialOperator b = Dummy(vf); REQUIRE(b.getName() == "Dummy"); - REQUIRE(b.getType() == SpatialOperator::Type::Explicit); + REQUIRE(b.getType() == Operator::Type::Explicit); } SECTION("Supports Coefficients Explicit " + execName) @@ -73,10 +73,10 @@ TEST_CASE("SpatialOperator") std::vector> bcs {}; auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - dsl::SpatialOperator b = Dummy(vf, SpatialOperator::Type::Implicit); + dsl::SpatialOperator b = Dummy(vf, Operator::Type::Implicit); REQUIRE(b.getName() == "Dummy"); - REQUIRE(b.getType() == SpatialOperator::Type::Implicit); + REQUIRE(b.getType() == Operator::Type::Implicit); auto ls = b.createEmptyLinearSystem(); REQUIRE(ls.matrix().nValues() == 1); @@ -93,12 +93,10 @@ TEST_CASE("SpatialOperator") BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - dsl::SpatialOperator c = - 2 * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); - dsl::SpatialOperator d = - fB * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); + dsl::SpatialOperator c = 2 * dsl::SpatialOperator(Dummy(vf, Operator::Type::Implicit)); + dsl::SpatialOperator d = fB * dsl::SpatialOperator(Dummy(vf, Operator::Type::Implicit)); dsl::SpatialOperator e = - Coeff(-3, fB) * dsl::SpatialOperator(Dummy(vf, SpatialOperator::Type::Implicit)); + Coeff(-3, fB) * dsl::SpatialOperator(Dummy(vf, Operator::Type::Implicit)); [[maybe_unused]] auto coeffC = c.getCoefficient(); [[maybe_unused]] auto coeffD = d.getCoefficient(); diff --git a/test/dsl/temporalOperator.cpp b/test/dsl/temporalOperator.cpp index a6d3b865e..7eb444829 100644 --- a/test/dsl/temporalOperator.cpp +++ b/test/dsl/temporalOperator.cpp @@ -28,7 +28,7 @@ TEST_CASE("TemporalOperator") dsl::TemporalOperator b = TemporalDummy(vf); REQUIRE(b.getName() == "TemporalDummy"); - REQUIRE(b.getType() == dsl::SpatialOperator::Type::Explicit); + REQUIRE(b.getType() == dsl::Operator::Type::Explicit); } SECTION("Supports Coefficients Explicit " + execName) @@ -75,10 +75,10 @@ TEST_CASE("TemporalOperator") std::vector> bcs {}; auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - dsl::TemporalOperator b = TemporalDummy(vf, SpatialOperator::Type::Implicit); + dsl::TemporalOperator b = TemporalDummy(vf, Operator::Type::Implicit); REQUIRE(b.getName() == "TemporalDummy"); - REQUIRE(b.getType() == SpatialOperator::Type::Implicit); + REQUIRE(b.getType() == Operator::Type::Implicit); auto ls = b.createEmptyLinearSystem(); REQUIRE(ls.matrix().nValues() == 1); @@ -97,10 +97,9 @@ TEST_CASE("TemporalOperator") BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); auto vf = VolumeField(exec, "vf", mesh, fA, bf, bcs); - auto c = 2 * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); - auto d = fB * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); - auto e = Coeff(-3, fB) - * dsl::TemporalOperator(TemporalDummy(vf, SpatialOperator::Type::Implicit)); + auto c = 2 * dsl::TemporalOperator(TemporalDummy(vf, Operator::Type::Implicit)); + auto d = fB * dsl::TemporalOperator(TemporalDummy(vf, Operator::Type::Implicit)); + auto e = Coeff(-3, fB) * dsl::TemporalOperator(TemporalDummy(vf, Operator::Type::Implicit)); [[maybe_unused]] auto coeffC = c.getCoefficient(); [[maybe_unused]] auto coeffD = d.getCoefficient(); diff --git a/test/finiteVolume/cellCentred/operator/divOperator.cpp b/test/finiteVolume/cellCentred/operator/divOperator.cpp index bfda71e11..241acd149 100644 --- a/test/finiteVolume/cellCentred/operator/divOperator.cpp +++ b/test/finiteVolume/cellCentred/operator/divOperator.cpp @@ -11,7 +11,7 @@ namespace fvcc = NeoFOAM::finiteVolume::cellCentred; -using SpatialOperator = NeoFOAM::dsl::SpatialOperator; +using Operator = NeoFOAM::dsl::Operator; TEST_CASE("DivOperator") { @@ -36,7 +36,7 @@ TEST_CASE("DivOperator") SECTION("Construct from Token" + execName) { NeoFOAM::Input input = NeoFOAM::TokenList({std::string("Gauss"), std::string("linear")}); - fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); + fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); } SECTION("Construct from Dictionary" + execName) @@ -45,6 +45,6 @@ TEST_CASE("DivOperator") {{std::string("DivOperator"), std::string("Gauss")}, {std::string("surfaceInterpolation"), std::string("linear")}} ); - fvcc::DivOperator(SpatialOperator::Type::Explicit, faceFlux, phi, input); + fvcc::DivOperator(Operator::Type::Explicit, faceFlux, phi, input); } } diff --git a/test/timeIntegration/rungeKutta.cpp b/test/timeIntegration/rungeKutta.cpp index 3db667d43..dd53e0cc6 100644 --- a/test/timeIntegration/rungeKutta.cpp +++ b/test/timeIntegration/rungeKutta.cpp @@ -33,9 +33,7 @@ class YSquared : public OperatorMixin public: - YSquared(VolumeField& field) - : OperatorMixin(field.exec(), field, SpatialOperator::Type::Explicit) - {} + YSquared(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} void explicitOperation(Field& source) const { diff --git a/test/timeIntegration/timeIntegration.cpp b/test/timeIntegration/timeIntegration.cpp index 500fed4b7..7ee7c3599 100644 --- a/test/timeIntegration/timeIntegration.cpp +++ b/test/timeIntegration/timeIntegration.cpp @@ -71,7 +71,7 @@ TEST_CASE("TimeIntegration") NeoFOAM::dsl::TemporalOperator ddtOperator = NeoFOAM::dsl::imp::ddt(vf); // ddt(U) = f - auto eqn = ddtOperator + dummy; + NeoFOAM::dsl::Expression eqn = ddtOperator + dummy; double dt {2.0}; double time {1.0}; From cf9539f6d5bee946e1e235e88900b0d5cd6fcd55 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Sun, 9 Feb 2025 21:32:35 +0100 Subject: [PATCH 13/21] removed enum Type from spatialOperator --- include/NeoFOAM/dsl/spatialOperator.hpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/include/NeoFOAM/dsl/spatialOperator.hpp b/include/NeoFOAM/dsl/spatialOperator.hpp index a1080c6c3..1867c04ee 100644 --- a/include/NeoFOAM/dsl/spatialOperator.hpp +++ b/include/NeoFOAM/dsl/spatialOperator.hpp @@ -50,13 +50,6 @@ class SpatialOperator { public: - enum class Type - { - Temporal, - Implicit, - Explicit - }; - template SpatialOperator(T cls) : model_(std::make_unique>(std::move(cls))) {} From 75da695b36da3ec6caae56761e51e4894ad32e9d Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 16:48:14 +0100 Subject: [PATCH 14/21] removed unnecessary return statements --- include/NeoFOAM/core/parallelAlgorithms.hpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/include/NeoFOAM/core/parallelAlgorithms.hpp b/include/NeoFOAM/core/parallelAlgorithms.hpp index 087b09f99..396728a8b 100644 --- a/include/NeoFOAM/core/parallelAlgorithms.hpp +++ b/include/NeoFOAM/core/parallelAlgorithms.hpp @@ -124,7 +124,7 @@ void parallelReduce( const NeoFOAM::Executor& exec, std::pair range, Kernel kernel, T& value ) { - return std::visit([&](const auto& e) { return parallelReduce(e, range, kernel, value); }, exec); + std::visit([&](const auto& e) { parallelReduce(e, range, kernel, value); }, exec); } @@ -159,9 +159,7 @@ void parallelReduce( template void parallelReduce(Field& field, Kernel kernel, T& value) { - return std::visit( - [&](const auto& e) { return parallelReduce(e, field, kernel, value); }, field.exec() - ); + std::visit([&](const auto& e) { parallelReduce(e, field, kernel, value); }, field.exec()); } template @@ -203,9 +201,7 @@ void parallelScan( ReturnType& returnValue ) { - return std::visit( - [&](const auto& e) { return parallelScan(e, range, kernel, returnValue); }, exec - ); + std::visit([&](const auto& e) { parallelScan(e, range, kernel, returnValue); }, exec); } From 501ecc486357f17c0d099e1615cd8dbbb3b4c28f Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 16:53:34 +0100 Subject: [PATCH 15/21] change interface of the ddtOperator expression API --- include/NeoFOAM/dsl/expression.hpp | 48 +++++++---- .../cellCentred/operators/ddtOperator.hpp | 10 ++- .../NeoFOAM/timeIntegration/backwardEuler.hpp | 84 +++++++++++++++++++ .../NeoFOAM/timeIntegration/forwardEuler.hpp | 4 +- .../NeoFOAM/timeIntegration/rungeKutta.hpp | 5 +- .../timeIntegration/timeIntegration.hpp | 19 +++-- src/timeIntegration/rungeKutta.cpp | 2 +- src/timeIntegration/timeIntegration.cpp | 3 + 8 files changed, 142 insertions(+), 33 deletions(-) create mode 100644 include/NeoFOAM/timeIntegration/backwardEuler.hpp diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp index df7a97217..0c88fe9ee 100644 --- a/include/NeoFOAM/dsl/expression.hpp +++ b/include/NeoFOAM/dsl/expression.hpp @@ -52,42 +52,54 @@ class Expression /* @brief perform all explicit operation and accumulate the result */ Field explicitOperation(Field& source) { - for (auto& oper : spatialOperators_) + for (auto& op : spatialOperators_) + { + if (op.getType() == Operator::Type::Explicit) + { + op.explicitOperation(source); + } + } + return source; + } + + Field explicitOperation(Field& source, scalar t, scalar dt) + { + for (auto& op : temporalOperators_) { - if (oper.getType() == Operator::Type::Explicit) + if (op.getType() == Operator::Type::Explicit) { - oper.explicitOperation(source); + op.explicitOperation(source, t, dt); } } - // for (auto& oper : temporalOperators_) - // { - // if (oper.getType() == Operator::Type::Explicit) - // { - // oper.explicitOperation(source); - // } - // } return source; } /* @brief perform all implicit operation and accumulate the result */ la::LinearSystem implicitOperation() { - // TODO: implement - // if (implicitOperators_.empty()) - // { - // NF_ERROR_EXIT("No implicit operators in the expression"); - // } auto ls = spatialOperators_[0].createEmptyLinearSystem(); - for (auto& oper : spatialOperators_) + for (auto& op : spatialOperators_) { - if (oper.getType() == Operator::Type::Implicit) + if (op.getType() == Operator::Type::Implicit) { - oper.implicitOperation(ls); + op.implicitOperation(ls); } } return ls; } + void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) + { + for (auto& op : temporalOperators_) + { + if (op.getType() == Operator::Type::Implicit) + { + op.implicitOperation(ls, t, dt); + } + } + } + + void addOperator(const SpatialOperator& oper) { spatialOperators_.push_back(oper); } void addOperator(const TemporalOperator& oper) { temporalOperators_.push_back(oper); } diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index 0b6eb5e39..e6fa05fab 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -29,6 +29,7 @@ class DdtOperator : public dsl::OperatorMixin> void explicitOperation(Field& source, scalar t, scalar dt) { const scalar rDeltat = 1 / dt; + const auto vol = getField().mesh().cellVolumes().span(); auto operatorScaling = getCoefficient(); auto [sourceSpan, field, oldField] = spans(source, field_.internalField(), oldTime(field_).internalField()); @@ -36,7 +37,7 @@ class DdtOperator : public dsl::OperatorMixin> source.exec(), source.range(), KOKKOS_LAMBDA(const size_t celli) { - sourceSpan[celli] += rDeltat * (field[celli] - oldField[celli]); + sourceSpan[celli] += rDeltat * (field[celli] - oldField[celli]) * vol[celli]; } ); } @@ -49,6 +50,7 @@ class DdtOperator : public dsl::OperatorMixin> void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) { const scalar rDeltat = 1 / dt; + const auto vol = getField().mesh().cellVolumes().span(); const auto operatorScaling = getCoefficient(); const auto [diagOffs, oldField] = spans(sparsityPattern_->diagOffset(), oldTime(field_).internalField()); @@ -61,8 +63,8 @@ class DdtOperator : public dsl::OperatorMixin> {0, oldField.size()}, KOKKOS_LAMBDA(const size_t celli) { std::size_t idx = rowPtrs[celli] + diagOffs[celli]; - values[idx] += operatorScaling[celli] * rDeltat; - rhs[celli] += operatorScaling[celli] * rDeltat * oldField[celli]; + values[idx] += operatorScaling[celli] * rDeltat * vol[celli]; + rhs[celli] += operatorScaling[celli] * rDeltat * oldField[celli] * vol[celli]; } ); } @@ -73,7 +75,7 @@ class DdtOperator : public dsl::OperatorMixin> // do nothing } - std::string getName() const { return "DivOperator"; } + std::string getName() const { return "DdtOperator"; } private: diff --git a/include/NeoFOAM/timeIntegration/backwardEuler.hpp b/include/NeoFOAM/timeIntegration/backwardEuler.hpp new file mode 100644 index 000000000..b01eebd32 --- /dev/null +++ b/include/NeoFOAM/timeIntegration/backwardEuler.hpp @@ -0,0 +1,84 @@ +// SPDX-License-Identifier: MIT +// +// SPDX-FileCopyrightText: 2023 NeoFOAM authors + +#pragma once + +#include "NeoFOAM/core/database/fieldCollection.hpp" +#include "NeoFOAM/core/database/oldTimeCollection.hpp" +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/timeIntegration/timeIntegration.hpp" +#include "NeoFOAM/dsl/solver.hpp" + +#if NF_WITH_GINKGO +#include "NeoFOAM/linearAlgebra/ginkgo.hpp" +#endif + +namespace NeoFOAM::timeIntegration +{ + +template +class BackwardEuler : + public TimeIntegratorBase::template Register< + BackwardEuler> +{ + +public: + + using Expression = NeoFOAM::dsl::Expression; + using Base = + TimeIntegratorBase::template Register>; + + BackwardEuler(const Dictionary& schemeDict, const Dictionary& solutionDict) + : Base(schemeDict, solutionDict) + {} + + static std::string name() { return "backwardEuler"; } + + static std::string doc() { return "first order time integration method"; } + + static std::string schema() { return "none"; } + + void solve( + Expression& eqn, SolutionFieldType& solutionField, [[maybe_unused]] scalar t, scalar dt + ) override + { + auto source = eqn.explicitOperation(solutionField.size()); + SolutionFieldType& oldSolutionField = + NeoFOAM::finiteVolume::cellCentred::oldTime(solutionField); + + // solutionField.internalField() = oldSolutionField.internalField() - source * dt; + // solutionField.correctBoundaryConditions(); + // solve sparse matrix system + using ValueType = typename SolutionFieldType::ElementType; + auto ls = eqn.implicitOperation(); + auto values = ls.matrix().values(); + eqn.implicitOperation(ls, t, dt); + auto ginkgoLs = NeoFOAM::dsl::ginkgoMatrix(ls, solutionField); + + NeoFOAM::Dictionary fvSolutionDict {}; + fvSolutionDict.insert("maxIters", 100); + fvSolutionDict.insert("relTol", float(1e-8)); + + auto internalFieldSpan = solutionField.internalField().span(); + + + NeoFOAM::la::ginkgo::BiCGStab solver(solutionField.exec(), fvSolutionDict); + solver.solve(ginkgoLs, solutionField.internalField()); + + // check if executor is GPU + if (std::holds_alternative(eqn.exec())) + { + Kokkos::fence(); + } + oldSolutionField.internalField() = solutionField.internalField(); + }; + + std::unique_ptr> clone() const override + { + return std::make_unique(*this); + } +}; + + +} // namespace NeoFOAM diff --git a/include/NeoFOAM/timeIntegration/forwardEuler.hpp b/include/NeoFOAM/timeIntegration/forwardEuler.hpp index aa172a87f..d6a41dac6 100644 --- a/include/NeoFOAM/timeIntegration/forwardEuler.hpp +++ b/include/NeoFOAM/timeIntegration/forwardEuler.hpp @@ -23,7 +23,9 @@ class ForwardEuler : using Base = TimeIntegratorBase::template Register>; - ForwardEuler(const Dictionary& dict) : Base(dict) {} + ForwardEuler(const Dictionary& schemeDict, const Dictionary& solutionDict) + : Base(schemeDict, solutionDict) + {} static std::string name() { return "forwardEuler"; } diff --git a/include/NeoFOAM/timeIntegration/rungeKutta.hpp b/include/NeoFOAM/timeIntegration/rungeKutta.hpp index e8724d4d2..e14c60703 100644 --- a/include/NeoFOAM/timeIntegration/rungeKutta.hpp +++ b/include/NeoFOAM/timeIntegration/rungeKutta.hpp @@ -58,7 +58,6 @@ class RungeKutta : using Expression = NeoFOAM::dsl::Expression; using Base = TimeIntegratorBase::template Register>; - using Base::dict_; /** * @brief Default constructor. @@ -74,7 +73,9 @@ class RungeKutta : * @brief Constructor that initializes the RungeKutta solver with a dictionary configuration. * @param dict The dictionary containing configuration parameters. */ - RungeKutta(const Dictionary& dict) : Base(dict) {} + RungeKutta(const Dictionary& schemeDict, const Dictionary& solutionDict) + : Base(schemeDict, solutionDict) + {} /** * @brief Copy constructor. diff --git a/include/NeoFOAM/timeIntegration/timeIntegration.hpp b/include/NeoFOAM/timeIntegration/timeIntegration.hpp index 233912dc1..6170b6e85 100644 --- a/include/NeoFOAM/timeIntegration/timeIntegration.hpp +++ b/include/NeoFOAM/timeIntegration/timeIntegration.hpp @@ -18,7 +18,9 @@ namespace NeoFOAM::timeIntegration */ template class TimeIntegratorBase : - public RuntimeSelectionFactory, Parameters> + public RuntimeSelectionFactory< + TimeIntegratorBase, + Parameters> { public: @@ -27,7 +29,9 @@ class TimeIntegratorBase : static std::string name() { return "timeIntegrationFactory"; } - TimeIntegratorBase(const Dictionary& dict) : dict_(dict) {} + TimeIntegratorBase(const Dictionary& schemeDict, const Dictionary& solutionDict) + : schemeDict_(schemeDict), solutionDict_(solutionDict) + {} virtual ~TimeIntegratorBase() {} @@ -40,7 +44,8 @@ class TimeIntegratorBase : protected: - const Dictionary& dict_; + const Dictionary& schemeDict_; + const Dictionary& solutionDict_; }; /** @@ -63,10 +68,10 @@ class TimeIntegration TimeIntegration(TimeIntegration&& timeIntegrator) : timeIntegratorStrategy_(std::move(timeIntegrator.timeIntegratorStrategy_)) {}; - TimeIntegration(const Dictionary& dict) - : timeIntegratorStrategy_( - TimeIntegratorBase::create(dict.get("type"), dict) - ) {}; + TimeIntegration(const Dictionary& schemeDict, const Dictionary& solutionDict) + : timeIntegratorStrategy_(TimeIntegratorBase::create( + schemeDict.get("type"), schemeDict, solutionDict + )) {}; void solve(Expression& eqn, SolutionFieldType& sol, scalar t, scalar dt) { diff --git a/src/timeIntegration/rungeKutta.cpp b/src/timeIntegration/rungeKutta.cpp index 4286d5d41..408f82479 100644 --- a/src/timeIntegration/rungeKutta.cpp +++ b/src/timeIntegration/rungeKutta.cpp @@ -131,7 +131,7 @@ void RungeKutta::initODEMemory(const scalar t) ERKStepSetTableNum( ark, NeoFOAM::sundials::stringToERKTable( - this->dict_.template get("Runge-Kutta-Method") + this->schemeDict_.template get("Runge-Kutta-Method") ) ); ARKodeSetUserData(ark, pdeExpr_.get()); diff --git a/src/timeIntegration/timeIntegration.cpp b/src/timeIntegration/timeIntegration.cpp index a55541487..e63afc23a 100644 --- a/src/timeIntegration/timeIntegration.cpp +++ b/src/timeIntegration/timeIntegration.cpp @@ -4,6 +4,7 @@ #include "NeoFOAM/timeIntegration/timeIntegration.hpp" #include "NeoFOAM/timeIntegration/forwardEuler.hpp" +#include "NeoFOAM/timeIntegration/backwardEuler.hpp" namespace fvcc = NeoFOAM::finiteVolume::cellCentred; @@ -12,4 +13,6 @@ namespace NeoFOAM::timeIntegration template class ForwardEuler>; +template class BackwardEuler>; + } // namespace NeoFOAM::dsl From 76eb70d49699417b9eaaf8cb7d2431677bbd1bb7 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 16:54:11 +0100 Subject: [PATCH 16/21] hardcoded upwind --- src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index b7eb58462..2e2c8ef6a 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -142,7 +142,8 @@ void GaussGreenDiv::div( {0, nInternalFaces}, KOKKOS_LAMBDA(const size_t facei) { scalar flux = sFaceFlux[facei]; - scalar weight = 0.5; + // scalar weight = 0.5; + scalar weight = flux >= 0 ? 1 : 0; scalar value = 0; std::size_t own = static_cast(owner[facei]); std::size_t nei = static_cast(neighbour[facei]); From 5bb17e64c21c7b97e63b385d642f91747ed7fd29 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 16:54:56 +0100 Subject: [PATCH 17/21] changed segmentedField --- include/NeoFOAM/fields/segmentedField.hpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/include/NeoFOAM/fields/segmentedField.hpp b/include/NeoFOAM/fields/segmentedField.hpp index 3738aa370..02deb0c12 100644 --- a/include/NeoFOAM/fields/segmentedField.hpp +++ b/include/NeoFOAM/fields/segmentedField.hpp @@ -14,7 +14,7 @@ namespace NeoFOAM * * The offsets are computed by a prefix sum of the input values. So, with given * input of {1, 2, 3, 4, 5} the offsets are {0, 1, 3, 6, 10, 15}. - * Note that the length of offSpan must be length of valSpan + 1 + * Note that the length of offSpan must be length of intervals + 1 * and are all elements of offSpan are required to be zero * * @param[in] in The values to compute the offsets from. @@ -24,14 +24,16 @@ template IndexType segmentsFromIntervals(const Field& intervals, Field& offsets) { IndexType finalValue = 0; - auto inSpan = intervals.span(); - auto offsSpan = offsets.span(); - NF_ASSERT_EQUAL(inSpan.size() + 1, offsSpan.size()); + const auto inSpan = intervals.span(); + // skip the first element of the offsets + // assumed to be zero + auto offsSpan = offsets.span().subspan(1); + NF_ASSERT_EQUAL(inSpan.size(), offsSpan.size()); NeoFOAM::parallelScan( intervals.exec(), - {1, offsSpan.size()}, - KOKKOS_LAMBDA(const std::size_t i, NeoFOAM::localIdx& update, const bool final) { - update += inSpan[i - 1]; + {0, offsSpan.size()}, + KOKKOS_LAMBDA(const std::size_t i, IndexType& update, const bool final) { + update += inSpan[i]; if (final) { offsSpan[i] = update; From 72be729d8546a4c4693348757da39de27cc019bc Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 16:55:19 +0100 Subject: [PATCH 18/21] added ginkgoMatrix conversion --- include/NeoFOAM/dsl/solver.hpp | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index c1716eb50..19f255301 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -22,6 +22,35 @@ namespace NeoFOAM::dsl { +template +NeoFOAM::la::LinearSystem ginkgoMatrix( + NeoFOAM::la::LinearSystem& ls, FieldType& solution +) +{ + using ValueType = typename FieldType::ElementType; + Field rhs(solution.exec(), ls.rhs().data(), ls.rhs().size()); + + Field mValues( + solution.exec(), ls.matrix().values().data(), ls.matrix().values().size() + ); + Field mColIdxs(solution.exec(), ls.matrix().colIdxs().size()); + auto mColIdxsSpan = ls.matrix().colIdxs(); + NeoFOAM::parallelFor( + mColIdxs, KOKKOS_LAMBDA(const size_t i) { return int(mColIdxsSpan[i]); } + ); + + Field mRowPtrs(solution.exec(), ls.matrix().rowPtrs().size()); + auto mRowPtrsSpan = ls.matrix().rowPtrs(); + NeoFOAM::parallelFor( + mRowPtrs, KOKKOS_LAMBDA(const size_t i) { return int(mRowPtrsSpan[i]); } + ); + + la::CSRMatrix matrix(mValues, mColIdxs, mRowPtrs); + + NeoFOAM::la::LinearSystem convertedLs(matrix, rhs, ls.sparsityPattern()); + return convertedLs; +} + template NeoFOAM::la::LinearSystem ginkgoMatrix(Expression& exp, FieldType& solution) @@ -95,7 +124,9 @@ void solve( if (exp.temporalOperators().size() > 0) { // integrate equations in time - timeIntegration::TimeIntegration timeIntegrator(fvSchemes.subDict("ddtSchemes")); + timeIntegration::TimeIntegration timeIntegrator( + fvSchemes.subDict("ddtSchemes"), fvSolution + ); timeIntegrator.solve(exp, solution, t, dt); } else From 5596901dfe49a7cc0eb31503eec34ad16e98846b Mon Sep 17 00:00:00 2001 From: Henning Scheufler <16381681+HenningScheufler@users.noreply.github.com> Date: Fri, 14 Feb 2025 16:59:11 +0100 Subject: [PATCH 19/21] Apply suggestions from code review Co-authored-by: Gregor Olenik --- include/NeoFOAM/dsl/ddt.hpp | 2 +- include/NeoFOAM/dsl/solver.hpp | 5 +---- include/NeoFOAM/dsl/spatialOperator.hpp | 2 +- include/NeoFOAM/dsl/temporalOperator.hpp | 2 +- .../finiteVolume/cellCentred/operators/ddtOperator.hpp | 1 - 5 files changed, 4 insertions(+), 8 deletions(-) diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp index b98452dee..ee8ca6b5c 100644 --- a/include/NeoFOAM/dsl/ddt.hpp +++ b/include/NeoFOAM/dsl/ddt.hpp @@ -30,7 +30,7 @@ class Ddt : public OperatorMixin NF_ERROR_EXIT("Not implemented"); } - void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) + void implicitOperation( [[maybe_unused]] la::LinearSystem& ls,[[maybe_unused]] scalar t, [[maybe_unused]] scalar dt) { NF_ERROR_EXIT("Not implemented"); } diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index 19f255301..c2a392372 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -57,8 +57,6 @@ ginkgoMatrix(Expression& exp, FieldType& solution) { using ValueType = typename FieldType::ElementType; auto vol = solution.mesh().cellVolumes().span(); - // la::CSRMatrix matrix(values, colIdxs, rowPtrs); - // return la::LinearSystem(matrix, rhs, sparsityPattern); auto expSource = exp.explicitOperation(solution.mesh().nCells()); auto expSourceSpan = expSource.span(); @@ -92,8 +90,7 @@ ginkgoMatrix(Expression& exp, FieldType& solution) la::CSRMatrix matrix(mValues, mColIdxs, mRowPtrs); - NeoFOAM::la::LinearSystem convertedLs(matrix, rhs, ls.sparsityPattern()); - return convertedLs; + return {matrix, rhs, ls.sparsityPattern()}; } /* @brief solve an expression diff --git a/include/NeoFOAM/dsl/spatialOperator.hpp b/include/NeoFOAM/dsl/spatialOperator.hpp index 1867c04ee..8b0c13519 100644 --- a/include/NeoFOAM/dsl/spatialOperator.hpp +++ b/include/NeoFOAM/dsl/spatialOperator.hpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors #pragma once #include diff --git a/include/NeoFOAM/dsl/temporalOperator.hpp b/include/NeoFOAM/dsl/temporalOperator.hpp index 4a739daeb..08c4560ff 100644 --- a/include/NeoFOAM/dsl/temporalOperator.hpp +++ b/include/NeoFOAM/dsl/temporalOperator.hpp @@ -1,5 +1,5 @@ // SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +// SPDX-FileCopyrightText: 2023-2025 NeoFOAM authors #pragma once #include diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index e6fa05fab..61ec36ef9 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -79,7 +79,6 @@ class DdtOperator : public dsl::OperatorMixin> private: - // const VolumeField& coefficients_; const std::shared_ptr sparsityPattern_; }; From f2691b15a6aec3d462d0c520a434ed90e18ed1bd Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 17:15:11 +0100 Subject: [PATCH 20/21] renamed files of sparsityPattern and linearSystem --- .../finiteVolume/cellCentred/operators/ddtOperator.hpp | 2 +- .../finiteVolume/cellCentred/operators/gaussGreenDiv.hpp | 2 +- .../operators/{fvccLinearSystem.hpp => linearSystem.hpp} | 9 ++++++--- .../finiteVolume/cellCentred/operators/sourceTerm.hpp | 2 +- .../{fvccSparsityPattern.hpp => sparsityPattern.hpp} | 0 src/CMakeLists.txt | 2 +- .../{fvccSparsityPattern.cpp => sparsityPattern.cpp} | 2 +- 7 files changed, 11 insertions(+), 8 deletions(-) rename include/NeoFOAM/finiteVolume/cellCentred/operators/{fvccLinearSystem.hpp => linearSystem.hpp} (93%) rename include/NeoFOAM/finiteVolume/cellCentred/operators/{fvccSparsityPattern.hpp => sparsityPattern.hpp} (100%) rename src/finiteVolume/cellCentred/operators/{fvccSparsityPattern.cpp => sparsityPattern.cpp} (98%) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index 61ec36ef9..d59ea8048 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -8,7 +8,7 @@ #include "NeoFOAM/core/input.hpp" #include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/mesh/unstructured.hpp" -#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp index eca993164..ce604e60d 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp @@ -8,7 +8,7 @@ #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp" #include "NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp" -#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp similarity index 93% rename from include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp rename to include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp index f06206493..66f291db2 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccLinearSystem.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp @@ -3,8 +3,11 @@ #pragma once -#include "NeoFOAM/linearAlgebra.hpp" -#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/linearAlgebra/CSRMatrix.hpp" +#include "NeoFOAM/linearAlgebra/linearSystem.hpp" +#include "NeoFOAM/linearAlgebra/ginkgo.hpp" + +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { @@ -56,7 +59,7 @@ class LinearSystem void diag(Field& field) { - NF_ASSERT_EQUAL(field.size() , sparsityPattern_->diagOffset().size()); + NF_ASSERT_EQUAL(field.size(), sparsityPattern_->diagOffset().size()); const auto diagOffset = sparsityPattern_->diagOffset().span(); const auto rowPtrs = ls_.matrix().rowPtrs(); std::span fieldSpan = field.span(); diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp index d1818b8f5..0bd1eb520 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp @@ -8,7 +8,7 @@ #include "NeoFOAM/core/input.hpp" #include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/mesh/unstructured.hpp" -#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp similarity index 100% rename from include/NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp rename to include/NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 97d3c2027..761477ef3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -38,7 +38,7 @@ target_sources( "finiteVolume/cellCentred/boundary/boundary.cpp" "finiteVolume/cellCentred/operators/gaussGreenGrad.cpp" "finiteVolume/cellCentred/operators/gaussGreenDiv.cpp" - "finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp" + "finiteVolume/cellCentred/operators/sparsityPattern.cpp" "finiteVolume/cellCentred/interpolation/linear.cpp" "finiteVolume/cellCentred/interpolation/upwind.cpp" "timeIntegration/timeIntegration.cpp" diff --git a/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp b/src/finiteVolume/cellCentred/operators/sparsityPattern.cpp similarity index 98% rename from src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp rename to src/finiteVolume/cellCentred/operators/sparsityPattern.cpp index 3cb963b72..7e4dc1842 100644 --- a/src/finiteVolume/cellCentred/operators/fvccSparsityPattern.cpp +++ b/src/finiteVolume/cellCentred/operators/sparsityPattern.cpp @@ -1,7 +1,7 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors -#include "NeoFOAM/finiteVolume/cellCentred/operators/fvccSparsityPattern.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" #include "NeoFOAM/fields/segmentedField.hpp" namespace NeoFOAM::finiteVolume::cellCentred From b316af0ebb7ffd4d21953393f41a0281ca54a1f4 Mon Sep 17 00:00:00 2001 From: Henning Scheufler Date: Fri, 14 Feb 2025 17:19:34 +0100 Subject: [PATCH 21/21] update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 402935f56..d60b93137 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,6 @@ # Version 0.2.0 (unreleased) ## Features +- add temporal operators inside expressions and the ability to solve linear system [#259 ](https://github.com/exasim-project/NeoFOAM/pull/259) - Add basic Ginkgo solver interface [#250](https://github.com/exasim-project/NeoFOAM/pull/250) - support for implicit operators in the DSL [#246](https://github.com/exasim-project/NeoFOAM/pull/246) - add segmentedField to represent vector of vectors [#202](https://github.com/exasim-project/NeoFOAM/pull/202)