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) 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); } 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..ee8ca6b5c 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,21 @@ class Ddt : public OperatorMixin public: - Ddt(FieldType& field) : OperatorMixin(field.exec(), field, Operator::Type::Temporal) + 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 scale) + void explicitOperation( + [[maybe_unused]] Field& source, + [[maybe_unused]] scalar t, + [[maybe_unused]] scalar dt + ) { NF_ERROR_EXIT("Not implemented"); } - void implicitOperation([[maybe_unused]] Field& phi) + 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/explicit.hpp b/include/NeoFOAM/dsl/explicit.hpp index e2fe0297e..ac33ea23f 100644 --- a/include/NeoFOAM/dsl/explicit.hpp +++ b/include/NeoFOAM/dsl/explicit.hpp @@ -5,23 +5,30 @@ #pragma once #include "NeoFOAM/fields/field.hpp" -#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/dsl/spatialOperator.hpp" #include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" #include "NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp" // TODO we should get rid of this include since it includes details // from a general implementation #include "NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp" namespace fvcc = NeoFOAM::finiteVolume::cellCentred; namespace NeoFOAM::dsl::exp { -Operator +SpatialOperator div(const fvcc::SurfaceField& faceFlux, fvcc::VolumeField& phi) { - return Operator(fvcc::DivOperator(dsl::Operator::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::Operator::Type::Explicit, coeff, phi)); } diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp index f67484808..0c88fe9ee 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,79 +52,82 @@ class Expression /* @brief perform all explicit operation and accumulate the result */ Field explicitOperation(Field& source) { - for (auto& oper : explicitOperators_) + for (auto& op : spatialOperators_) { - oper.explicitOperation(source); + if (op.getType() == Operator::Type::Explicit) + { + op.explicitOperation(source); + } } return source; } - /* @brief perform all implicit operation and accumulate the result */ - la::LinearSystem implicitOperation() + Field explicitOperation(Field& source, scalar t, scalar dt) { - if (implicitOperators_.empty()) + for (auto& op : temporalOperators_) { - NF_ERROR_EXIT("No implicit operators in the expression"); + if (op.getType() == Operator::Type::Explicit) + { + op.explicitOperation(source, t, dt); + } } - auto ls = implicitOperators_[0].createEmptyLinearSystem(); - for (auto& oper : implicitOperators_) + return source; + } + + /* @brief perform all implicit operation and accumulate the result */ + la::LinearSystem implicitOperation() + { + auto ls = spatialOperators_[0].createEmptyLinearSystem(); + for (auto& op : spatialOperators_) { - oper.implicitOperation(ls); + if (op.getType() == Operator::Type::Implicit) + { + op.implicitOperation(ls); + } } return ls; } - void addOperator(const Operator& oper) + void implicitOperation(la::LinearSystem& ls, scalar t, scalar dt) { - switch (oper.getType()) + for (auto& op : temporalOperators_) { - 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; + 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); } + void addExpression(const Expression& equation) { for (auto& oper : equation.temporalOperators_) { 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 +135,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 +146,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 +176,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 +189,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 0a19d9df3..1a80ae7e1 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,6 +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::Operator::Type::Implicit, phi); +} + +SpatialOperator +Source(fvcc::VolumeField& coeff, fvcc::VolumeField& phi) +{ + return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Implicit, coeff, phi)); +} + +SpatialOperator +div(fvcc::SurfaceField& faceFlux, fvcc::VolumeField& 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 index a597f6e79..12d8e5af4 100644 --- a/include/NeoFOAM/dsl/operator.hpp +++ b/include/NeoFOAM/dsl/operator.hpp @@ -2,277 +2,18 @@ // 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" - -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) { - { - t.explicitOperation(std::declval&>()) - } -> std::same_as; // Adjust return type and arguments as needed -}; - -template -concept HasImplicitOperator = requires(T t) { - { - t.implicitOperation(std::declval&>()) - } -> std::same_as; // Adjust return type and arguments as needed -}; - -/* @class Operator - * @brief A class to represent an operator 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 Operators e.g Divergence, Laplacian, etc can be stored in a vector of - * Operators - * - * @ingroup dsl - */ class Operator { public: enum class Type { - Temporal, Implicit, Explicit }; - - template - Operator(T cls) : model_(std::make_unique>(std::move(cls))) - {} - - Operator(const Operator& eqnOperator); - - Operator(Operator&& 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; - - 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 OperatorConcept - { - virtual ~OperatorConcept() = default; - - 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; - - /* @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 Operator::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 OperatorModel : OperatorConcept - { - /* @brief build with concrete operator */ - OperatorModel(ConcreteOperatorType concreteOp) : concreteOp_(std::move(concreteOp)) {} - - /* returns the name of the operator */ - std::string getName() const override { return concreteOp_.getName(); } - - virtual void explicitOperation(Field& source) override - { - if constexpr (HasExplicitOperator) - { - concreteOp_.explicitOperation(source); - } - } - - virtual void temporalOperation(Field& field) override - { - if constexpr (HasTemporalOperator) - { - concreteOp_.temporalOperation(field); - } - } - - virtual void implicitOperation(la::LinearSystem& ls) override - { - if constexpr (HasImplicitOperator) - { - concreteOp_.implicitOperation(ls); - } - } - - virtual la::LinearSystem createEmptyLinearSystem() const override - { - if constexpr (HasImplicitOperator) - { - 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 */ - Operator::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); - } - - ConcreteOperatorType concreteOp_; - }; - - std::unique_ptr model_; }; - -Operator operator*(scalar scalarCoeff, Operator rhs); - -Operator operator*(const Field& coeffField, Operator rhs); - -Operator operator*(const Coeff& coeff, Operator rhs); - -template - requires std::invocable -Operator operator*([[maybe_unused]] CoeffFunction coeffFunc, const Operator& lhs) -{ - // TODO implement - NF_ERROR_EXIT("Not implemented"); - Operator 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; -} - -/* @class OperatorMixin - * @brief A mixin class to simplify implementations of concrete operators - * in NeoFOAMs dsl - * - * @ingroup dsl - */ -template -class OperatorMixin -{ - -public: - - OperatorMixin(const Executor exec, FieldType& field, Operator::Type type) - : exec_(exec), coeffs_(), field_(field), type_(type) {}; - - Operator::Type getType() const { return type_; } - - virtual ~OperatorMixin() = default; - - virtual const Executor& exec() const final { return exec_; } - - Coeff& getCoefficient() { return coeffs_; } - - const Coeff& getCoefficient() const { return coeffs_; } - - FieldType& getField() { return field_; } - - const FieldType& getField() const { return field_; } - - /* @brief Given an input this function reads required coeffs */ - void build([[maybe_unused]] const Input& input) {} - -protected: - - const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.) - - Coeff coeffs_; - - FieldType& field_; - - Operator::Type type_; -}; } // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/solver.hpp b/include/NeoFOAM/dsl/solver.hpp index b0cb0f7d4..c2a392372 100644 --- a/include/NeoFOAM/dsl/solver.hpp +++ b/include/NeoFOAM/dsl/solver.hpp @@ -14,10 +14,85 @@ #include "NeoFOAM/dsl/expression.hpp" #include "NeoFOAM/timeIntegration/timeIntegration.hpp" +#if NF_WITH_GINKGO +#include "NeoFOAM/linearAlgebra/ginkgo.hpp" +#endif + 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) +{ + using ValueType = typename FieldType::ElementType; + auto vol = solution.mesh().cellVolumes().span(); + 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); + + return {matrix, rhs, ls.sparsityPattern()}; +} + /* @brief solve an expression * * @param exp - Expression which is to be solved/updated. @@ -37,7 +112,8 @@ void solve( [[maybe_unused]] const Dictionary& fvSolution ) { - if (exp.temporalOperators().size() == 0 && exp.implicitOperators().size() == 0) + // FIXME: + if (exp.temporalOperators().size() == 0 && exp.spatialOperators().size() == 0) { NF_ERROR_EXIT("No temporal or implicit terms to solve."); } @@ -45,13 +121,20 @@ 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 { - 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/dsl/spatialOperator.hpp b/include/NeoFOAM/dsl/spatialOperator.hpp new file mode 100644 index 000000000..8b0c13519 --- /dev/null +++ b/include/NeoFOAM/dsl/spatialOperator.hpp @@ -0,0 +1,258 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2025 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/operator.hpp" + +namespace la = NeoFOAM::la; + +namespace NeoFOAM::dsl +{ + +template +concept HasExplicitOperator = requires(T t) { + { + t.explicitOperation(std::declval&>()) + } -> std::same_as; // Adjust return type and arguments as needed +}; + +template +concept HasImplicitOperator = requires(T t) { + { + t.implicitOperation(std::declval&>()) + } -> 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 + * + * 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 Operators e.g Divergence, Laplacian, etc can be stored in a vector of + * Operators + * + * @ingroup dsl + */ +class SpatialOperator +{ +public: + + template + SpatialOperator(T cls) : model_(std::make_unique>(std::move(cls))) + {} + + SpatialOperator(const SpatialOperator& eqnOperator); + + SpatialOperator(SpatialOperator&& eqnOperator); + + void explicitOperation(Field& source); + + void implicitOperation(la::LinearSystem& ls); + + la::LinearSystem createEmptyLinearSystem() const; + + /* returns the fundamental type of an operator, ie explicit, implicit */ + Operator::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 OperatorConcept + { + virtual ~OperatorConcept() = default; + + virtual void explicitOperation(Field& source) = 0; + + virtual void implicitOperation(la::LinearSystem& ls) = 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 */ + virtual Operator::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 OperatorModel : OperatorConcept + { + /* @brief build with concrete operator */ + OperatorModel(ConcreteOperatorType concreteOp) : concreteOp_(std::move(concreteOp)) {} + + /* returns the name of the operator */ + std::string getName() const override { return concreteOp_.getName(); } + + virtual void explicitOperation(Field& source) override + { + if constexpr (HasExplicitOperator) + { + concreteOp_.explicitOperation(source); + } + } + + virtual void implicitOperation(la::LinearSystem& ls) override + { + if constexpr (HasImplicitOperator) + { + concreteOp_.implicitOperation(ls); + } + } + + virtual la::LinearSystem createEmptyLinearSystem() const override + { + if constexpr (HasImplicitOperator) + { + 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 */ + Operator::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); + } + + ConcreteOperatorType concreteOp_; + }; + + std::unique_ptr model_; +}; + + +SpatialOperator operator*(scalar scalarCoeff, SpatialOperator rhs); + +SpatialOperator operator*(const Field& coeffField, SpatialOperator rhs); + +SpatialOperator operator*(const Coeff& coeff, SpatialOperator rhs); + +template + requires std::invocable +SpatialOperator operator*([[maybe_unused]] CoeffFunction coeffFunc, const SpatialOperator& lhs) +{ + // TODO implement + NF_ERROR_EXIT("Not implemented"); + SpatialOperator 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; +} + +/* @class OperatorMixin + * @brief A mixin class to simplify implementations of concrete operators + * in NeoFOAMs dsl + * + * @ingroup dsl + */ +template +class OperatorMixin +{ + +public: + + OperatorMixin(const Executor exec, FieldType& field, Operator::Type type) + : exec_(exec), coeffs_(), field_(field), type_(type) {}; + + Operator::Type getType() const { return type_; } + + virtual ~OperatorMixin() = default; + + virtual const Executor& exec() const final { return exec_; } + + Coeff& getCoefficient() { return coeffs_; } + + const Coeff& getCoefficient() const { return coeffs_; } + + FieldType& getField() { return field_; } + + const FieldType& getField() const { return field_; } + + /* @brief Given an input this function reads required coeffs */ + void build([[maybe_unused]] const Input& input) {} + +protected: + + const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.) + + Coeff coeffs_; + + FieldType& field_; + + Operator::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..08c4560ff --- /dev/null +++ b/include/NeoFOAM/dsl/temporalOperator.hpp @@ -0,0 +1,231 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2025 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/operator.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 */ + Operator::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 Operator::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 */ + Operator::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/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; 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/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. * diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp index f4ad512bc..d59ea8048 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp @@ -6,9 +6,9 @@ #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" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { @@ -22,10 +22,14 @@ 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) + }; + + 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()); @@ -33,7 +37,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]) * vol[celli]; } ); } @@ -43,19 +47,24 @@ 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 / dt; + const auto vol = getField().mesh().cellVolumes().span(); 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 * vol[celli]; + rhs[celli] += operatorScaling[celli] * rDeltat * oldField[celli] * vol[celli]; } ); } @@ -66,11 +75,10 @@ class DdtOperator : public dsl::OperatorMixin> // do nothing } - std::string getName() const { return "DivOperator"; } + std::string getName() const { return "DdtOperator"; } private: - const VolumeField& coefficients_; const std::shared_ptr sparsityPattern_; }; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp index 80419b81d..f6baa9984 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" @@ -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..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 { @@ -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/include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp new file mode 100644 index 000000000..66f291db2 --- /dev/null +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/linearSystem.hpp @@ -0,0 +1,138 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2025 NeoFOAM authors + +#pragma once + +#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 +{ + +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) + { + 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(); + 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; +} + +} diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp index 1aa4d9378..0bd1eb520 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/sourceTerm.hpp @@ -6,9 +6,9 @@ #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" +#include "NeoFOAM/finiteVolume/cellCentred/operators/sparsityPattern.hpp" namespace NeoFOAM::finiteVolume::cellCentred { @@ -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]; } ); } 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/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/CMakeLists.txt b/src/CMakeLists.txt index cd31ab41f..761477ef3 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" @@ -37,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/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..69fc1ed82 --- /dev/null +++ b/src/dsl/spatialOperator.cpp @@ -0,0 +1,67 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2025 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(); +} + +Operator::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..0bba489c2 --- /dev/null +++ b/src/dsl/temporalOperator.cpp @@ -0,0 +1,70 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2025 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(); +} + +Operator::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/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index 0bf38e40c..2e2c8ef6a 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 ) @@ -136,23 +142,27 @@ 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]); - // 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); + 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 - std::size_t rowOwnStart = rowPtrs[own]; - values[rowOwnStart + ownOffs[facei]] += flux * (1 - weight); - Kokkos::atomic_sub(&values[rowOwnStart + diagOffs[own]], flux * (1 - weight)); + // add owner contribution lower + value = flux * (1 - weight); + values[rowOwnStart + ownOffs[facei]] += value; + Kokkos::atomic_sub(&values[rowNeiStart + diagOffs[nei]], value); } ); }; 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 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/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 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..5f14ffdaf 100644 --- a/test/dsl/common.hpp +++ b/test/dsl/common.hpp @@ -19,8 +19,8 @@ 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 { @@ -83,6 +83,77 @@ 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, Operator::Type::Explicit) + {} + + TemporalDummy(VolumeField& field, Operator::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/operator.cpp b/test/dsl/spatialOperator.cpp similarity index 83% rename from test/dsl/operator.cpp rename to test/dsl/spatialOperator.cpp index 1d2f1ed31..ec15b597d 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,14 +18,14 @@ 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); @@ -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,7 +73,7 @@ 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, Operator::Type::Implicit); REQUIRE(b.getName() == "Dummy"); REQUIRE(b.getType() == Operator::Type::Implicit); @@ -91,9 +93,10 @@ 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, Operator::Type::Implicit)); + dsl::SpatialOperator d = fB * dsl::SpatialOperator(Dummy(vf, Operator::Type::Implicit)); + dsl::SpatialOperator e = + 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 new file mode 100644 index 000000000..7eb444829 --- /dev/null +++ b/test/dsl/temporalOperator.cpp @@ -0,0 +1,136 @@ +// 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::Operator::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, Operator::Type::Implicit); + + REQUIRE(b.getName() == "TemporalDummy"); + REQUIRE(b.getType() == Operator::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, 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(); + [[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/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..500fed4b7 --- /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); + NeoFOAM::dsl::TemporalOperator ddtOperator = NeoFOAM::dsl::imp::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); + } +} diff --git a/test/timeIntegration/rungeKutta.cpp b/test/timeIntegration/rungeKutta.cpp index 5267153ef..dd53e0cc6 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; @@ -119,7 +120,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..7ee7c3599 100644 --- a/test/timeIntegration/timeIntegration.cpp +++ b/test/timeIntegration/timeIntegration.cpp @@ -68,10 +68,10 @@ 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; + NeoFOAM::dsl::Expression eqn = ddtOperator + dummy; double dt {2.0}; double time {1.0};