diff --git a/.github/workflows/build_doc.yaml b/.github/workflows/build_doc.yaml index 1ee7a415f..1f9c680d4 100644 --- a/.github/workflows/build_doc.yaml +++ b/.github/workflows/build_doc.yaml @@ -5,8 +5,6 @@ env: on: pull_request: - branches: - - main types: [closed, synchronize, opened] push: tags: @@ -50,6 +48,13 @@ jobs: with: folder: docs_build target-folder: latest + - name: Comment PR + uses: thollander/actions-comment-pull-request@v2 + with: + message: | + Deployed test documentation to https://exasim-project.com/NeoFOAM/Build_PR_${{ env.PR_NUMBER }} + comment_tag: build_url + - name: Echo Build PR URL run: | echo "Deploy to: https://exasim-project.com/NeoFOAM/Build_PR_${{ env.PR_NUMBER }}" diff --git a/CHANGELOG.md b/CHANGELOG.md index 98654af6c..71fac7819 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,5 @@ # Version 0.1.0 (unreleased) +- Implement a basic DSL interface [#102](https://github.com/exasim-project/NeoFOAM/pull/102) - faster project configuration [#179](https://github.com/exasim-project/NeoFOAM/pull/179) - improved error handling and addition of tokenList and Input [#134](https://github.com/exasim-project/NeoFOAM/pull/134) - disable span from temporary objects and simplification related to fields [#139](https://github.com/exasim-project/NeoFOAM/pull/139) diff --git a/doc/basics/fields.rst b/doc/basics/fields.rst index c183955bb..44b44ad19 100644 --- a/doc/basics/fields.rst +++ b/doc/basics/fields.rst @@ -9,7 +9,7 @@ Fields Cell Centered Fields ^^^^^^^^^^^^^^^^^^^^ -The ``VolumeField`` stores the field values at cell centers and along boundaries, providing essential data for constructing the DSL (Domain Specific Language). This functionality also includes access to mesh data, integrating closely with the computational framework. +The ``VolumeField`` stores the field values at cell centers and along boundaries, providing essential data for constructing the Domain Specific Language (DSL). This functionality also includes access to mesh data, integrating closely with the computational framework. ``DomainField`` acts as the fundamental data container within this structure, offering both read and write to the ``internalField`` and ``boundaryFields`` provided by the ``DomainField``. The ``correctBoundaryConditions`` member function updates the field's boundary conditions, which are specified at construction. It does not hold the data but rather modifies the ``DomainField`` or ``BoundaryField`` container. @@ -26,7 +26,7 @@ Functionally, ``fvccVolField`` parallels several OpenFOAM classes such as ``volS Face Centered fields ^^^^^^^^^^^^^^^^^^^^ -The ``SurfaceField`` class stores the field values interpreted as face centers values. Additionally, it stores boundaries for the corresponding boundary conditions. This provides essential data for constructing the DSL (Domain Specific Language). The functionality also includes access to mesh data, integrating closely with the computational framework. +The ``SurfaceField`` class stores the field values interpreted as face centers values. Additionally, it stores boundaries for the corresponding boundary conditions. This provides essential data for constructing the DSL. The functionality also includes access to mesh data, integrating closely with the computational framework. ``DomainField`` acts as the fundamental data container within this structure, offering both read and to the ``internalField`` and ``boundaryField`` provided by the ``DomainField``. The ``correctBoundaryConditions`` member function updates field's boundary conditions, which are specified at construction. It does not hold the data, but modify the ``DomainField`` or ``BoundaryField`` container. diff --git a/doc/conf.py b/doc/conf.py index 731d51a8b..56e66dd30 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -23,6 +23,7 @@ extensions = [ 'sphinx.ext.autodoc', + 'sphinxcontrib.mermaid', 'sphinx.ext.intersphinx', 'sphinx.ext.autosectionlabel', 'sphinx.ext.todo', diff --git a/doc/dsl/equation.rst b/doc/dsl/equation.rst new file mode 100644 index 000000000..ea02c56b2 --- /dev/null +++ b/doc/dsl/equation.rst @@ -0,0 +1,50 @@ +Expression +--------- + + +The `Expression` class in NeoFOAM holds a set of operators to express an expression and ultimately helps to formulate equations. +Its core responsibility lies in the answering of the following questions: + + - How to discretize the spatial terms? + - In OpenFOAM this information is provided in **fvSchemes** + - How to integrate the system in time? + - In OpenFOAM this information is provided in **fvSchemes** + - How to solve the system? + - In OpenFOAM this information is provided in **fvSolution** + +The main difference between NeoFOAM and OpenFOAM is that the DSL is evaluated lazily, i.e. evaluation is not performed on construction by default. +Since, the evaluation is not tied to the construction but rather delayed, other numerical integration strategies (e.g. RK methods or even sub-stepping within an the equation) are possible. + +To implement lazy evaluation, the `Expression` stores 3 vectors: + +.. mermaid:: + + classDiagram + class Operator { + +explicitOperation(...) + +implicitOperation(...) + } + class DivOperator { + +explicitOperation(...) + +implicitOperation(...) + } + class TemporalOperator { + +explicitOperation(...) + +implicitOperation(...) + } + class Others["..."] { + +explicitOperation(...) + +implicitOperation(...) + } + class Expression { + +temporalTerms_: vector~Operator~ + +implicitTerms_: vector~Operator~ + +explicitTerms_: vector~Operator~ + } + Operator <|-- DivOperator + Operator <|-- TemporalOperator + Operator <|-- Others + Expression <|-- Operator + +Thus, an `Expression` consists of multiple `Operators` which are either explicit, implicit, or temporal. +Consequently, addition, subtraction, and scaling with a field needs to be handled by the `Operator`. diff --git a/doc/dsl/index.rst b/doc/dsl/index.rst new file mode 100644 index 000000000..586bf2c46 --- /dev/null +++ b/doc/dsl/index.rst @@ -0,0 +1,71 @@ +.. _fvcc: + +Domain Specific Language (DSL) +============================== + +The concept of a Domain Specific Language (DSL) allows to simplify the process of implementing and solving equations in a given programming language like C++. +Engineers can express equations in a concise and readable form, closely resembling their mathematical representation, while no or little knowledge of the numerical schemes and implementation is required. +This approach allows engineers to focus on the physics of the problem rather than the numerical implementation and helps in reducing the time and effort required to develop and maintain complex simulations. + +The Navier-Stokes equations can be expressed in the DSL in OpenFOAM as follows: + +.. code-block:: cpp + + fvVectorMatrix UEqn + ( + fvm::ddt(U) + + fvm::div(phi, U) + - fvm::laplacian(nu, U) + ); + + solve(UEqn == -fvc::grad(p)); + +To solve the continuity equation in OpenFOAM with the PISO or SIMPLE algorithm, the VectorMatrix, UEqn, needs to provide the diagonal and off-diagonal terms of the matrix. + +.. code-block:: cpp + + volScalarField rAU(1.0/UEqn.A()); + volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +This approach is readable and easy to understand for engineers familiar with OpenFOAM. However, it has several limitations due to its implementation in OpenFOAM: + + - the solution system is always a sparse matrix + - individual terms are eagerly evaluated, resulting in unnecessary use of temporaries + - the sparse matrix is always an LDU matrix that is not supported by external linear solvers + - Only cell-centred discretisation is supported + +NeoFOAM DSL tries to address these issues by providing: + + - lazy evaluation of the equations terms. This allows for better optimisation of the resulting equation and can reduce the number of required temporaries. + - a more modular design + - Support for standard matrix formats like COO and CSR, to simplify the use of external linear solvers + +The use of standard matrix formats combined with lazy evaluation allows for the use of external libraries to integrate PDEs in time and space. +The equation system can be passed to **sundials** and be integrated by **RK methods** and **BDF methods** on heterogeneous architectures. + +The NeoFOAM DSL is designed as drop in replacement for OpenFOAM DSL and the adoption should be possible with minimal effort and the same equation from above should read: + +.. code-block:: cpp + + dsl::Equation UEqn + ( + dsl::temporal::ddt(U) + + dsl::implicit::div(phi, U) + - dsl::implicit::laplacian(nu, U) + ) + + solve(UEqn == -dsl::explicit::grad(p), U, fvSchemes, fvSolution); + + +In contrast to OpenFOAM, the matrix assembly is deferred till the solve step. Hence the majority of the computational work is performed during the solve step. +That is 1. assemble the system and 2. solve the system. +After the system is assembled or solved, it provides access to the linear system for the SIMPLE and PISO algorithms. + + + +.. toctree:: + :maxdepth: 2 + :glob: + + equation.rst + operator.rst diff --git a/doc/dsl/operator.rst b/doc/dsl/operator.rst new file mode 100644 index 000000000..12adb5780 --- /dev/null +++ b/doc/dsl/operator.rst @@ -0,0 +1,51 @@ +Operator +======= + + +The `Operator` class represents a term in an equation and can be instantiated with different value types. +An `Operator` is either explicit, implicit or temporal, and can be scalable by an additional coefficient, for example a scalar value or a further field. +The `Operator` implementation uses Type Erasure (more details `[1] `_ `[2] `_ `[3] `_) to achieve polymorphism without inheritance. Consequently, the class needs only to implement the interface which is used in the DSL and which is shown in the below example: + +Example: + .. code-block:: cpp + + NeoFOAM::dsl::Operator divTerm = + Divergence(NeoFOAM::dsl::Operator::Type::Explicit, exec, ...); + + NeoFOAM::dsl::Operator ddtTerm = + TimeTerm(NeoFOAM::dsl::Operator::Type::Temporal, exec, ..); + + +To fit the specification of the EqnSystem (storage in a vector), the Operator needs to be able to be scaled: + +.. code-block:: cpp + + NeoFOAM::Field scalingField(exec, nCells, 2.0); + auto sF = scalingField.span(); + + dsl::Operator customTerm = + CustomTerm(dsl::Operator::Type::Explicit, exec, nCells, 1.0); + + auto constantScaledTerm = 2.0 * customTerm; // A constant scaling factor of 2 for the term. + auto fieldScaledTerm = scalingField * customTerm; // scalingField is used to scale the term. + + // Operator also supports a similar syntax as OpenFOAM + auto multiScaledTerm = (scale + scale + scale + scale) * customTerm; + + // Operator also supports the use of a lambda as scaling function to reduce the number of temporaries generated + auto lambdaScaledTerm = + (KOKKOS_LAMBDA(const NeoFOAM::size_t i) { return sF[i] + sF[i] + sF[i] + sF[i]; }) * customTerm; + +To add a user-defined `Operator`, a new derived class must be created, inheriting from `OperatorMixin`, + and provide the definitions of the below virtual functions that are required for the `Operator` interface: + + - build: build the term + - explicitOperation: perform the explicit operation + - implicitOperation: perform the implicit operation + - getType: get the type of the term + - exec: get the executor + - volumeField: get the volume field + +An example can be found in `test/dsl/operator.cpp`. + +The required scaling of the term is handled by the `Coeff` type which can be retrieved by the `getCoefficient` function of `Operator`. diff --git a/doc/finiteVolume/DSL.rst b/doc/finiteVolume/DSL.rst deleted file mode 100644 index 514441de9..000000000 --- a/doc/finiteVolume/DSL.rst +++ /dev/null @@ -1,4 +0,0 @@ -.. _fvcc_DSL: - -Domain Specific Language -======================== diff --git a/doc/finiteVolume/cellCentred/index.rst b/doc/finiteVolume/cellCentred/index.rst index 21ce5a6d2..c8c37828f 100644 --- a/doc/finiteVolume/cellCentred/index.rst +++ b/doc/finiteVolume/cellCentred/index.rst @@ -8,8 +8,6 @@ cellCenteredFiniteVolume :maxdepth: 2 :glob: - DSL.rst - fields.rst boundaryConditions.rst operators.rst stencil.rst diff --git a/doc/index.rst b/doc/index.rst index 0d7e3fef6..243b96f9e 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -26,6 +26,7 @@ Table of Contents contributing basics/index finiteVolume/cellCentred/index + dsl/index api/index Compatibility with OpenFOAM diff --git a/doc/requirements.txt b/doc/requirements.txt index a49f2dc6c..f7520c582 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -2,3 +2,4 @@ sphinx sphinx-sitemap furo breathe +sphinxcontrib-mermaid diff --git a/include/NeoFOAM/core/demangle.hpp b/include/NeoFOAM/core/demangle.hpp index 818d2202d..afc450037 100644 --- a/include/NeoFOAM/core/demangle.hpp +++ b/include/NeoFOAM/core/demangle.hpp @@ -1,11 +1,13 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2024 NeoFOAM authors + #pragma once -// TODO For WIN builds, needs to be ifdef'ed out. -#include -#include -#include -#include + +#include // for bad_any_cast +#include // for operator<<, basic_ostream, cerr, endl, ostream +#include // for operator<<, string +#include // for type_info + namespace NeoFOAM { diff --git a/include/NeoFOAM/core/dictionary.hpp b/include/NeoFOAM/core/dictionary.hpp index cd06b2dc9..f7230ae73 100644 --- a/include/NeoFOAM/core/dictionary.hpp +++ b/include/NeoFOAM/core/dictionary.hpp @@ -5,11 +5,9 @@ #include #include #include -#include #include #include "NeoFOAM/core/demangle.hpp" -#include "NeoFOAM/core/error.hpp" namespace NeoFOAM { @@ -162,6 +160,12 @@ class Dictionary */ const std::unordered_map& getMap() const; + /** + * @brief Checks whether the dictionary is empty + * @return A bool indicating if the dictionary is empty + */ + bool empty() const { return data_.empty(); } + private: std::unordered_map data_; diff --git a/include/NeoFOAM/core/executor/CPUExecutor.hpp b/include/NeoFOAM/core/executor/CPUExecutor.hpp index 36e0bc255..c51362fe8 100644 --- a/include/NeoFOAM/core/executor/CPUExecutor.hpp +++ b/include/NeoFOAM/core/executor/CPUExecutor.hpp @@ -2,8 +2,7 @@ // SPDX-FileCopyrightText: 2023 NeoFOAM authors #pragma once -#include -#include +#include // IWYU pragma: keep namespace NeoFOAM { diff --git a/include/NeoFOAM/core/executor/GPUExecutor.hpp b/include/NeoFOAM/core/executor/GPUExecutor.hpp index d98ac974e..7ef88a0e1 100644 --- a/include/NeoFOAM/core/executor/GPUExecutor.hpp +++ b/include/NeoFOAM/core/executor/GPUExecutor.hpp @@ -3,7 +3,6 @@ #pragma once #include -#include namespace NeoFOAM { diff --git a/include/NeoFOAM/core/executor/serialExecutor.hpp b/include/NeoFOAM/core/executor/serialExecutor.hpp index b68cfe960..c2e7b9248 100644 --- a/include/NeoFOAM/core/executor/serialExecutor.hpp +++ b/include/NeoFOAM/core/executor/serialExecutor.hpp @@ -3,7 +3,6 @@ #pragma once #include -#include namespace NeoFOAM { diff --git a/include/NeoFOAM/core/inputs.hpp b/include/NeoFOAM/core/input.hpp similarity index 100% rename from include/NeoFOAM/core/inputs.hpp rename to include/NeoFOAM/core/input.hpp diff --git a/include/NeoFOAM/core/primitives/vector.hpp b/include/NeoFOAM/core/primitives/vector.hpp index 9d6f7133a..1843adefd 100644 --- a/include/NeoFOAM/core/primitives/vector.hpp +++ b/include/NeoFOAM/core/primitives/vector.hpp @@ -1,8 +1,9 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors + #pragma once -#include +#include // IWYU pragma: keep #include "NeoFOAM/core/primitives/scalar.hpp" #include "NeoFOAM/core/primitives/label.hpp" diff --git a/include/NeoFOAM/dsl/coeff.hpp b/include/NeoFOAM/dsl/coeff.hpp new file mode 100644 index 000000000..cd91ecdb7 --- /dev/null +++ b/include/NeoFOAM/dsl/coeff.hpp @@ -0,0 +1,110 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023 NeoFOAM authors +#pragma once + +namespace NeoFOAM::dsl +{ + +/** + * @class Coeff + * @brief A class that represents a coefficient for the NeoFOAM dsl. + * + * This class stores a single scalar coefficient and optionally span of values. + * It is used to delay the evaluation of a scalar multiplication with a field to + * avoid the creation of a temporary field copy. + * It provides an indexing operator `operator[]` that returns the evaluated value at the specified + * index. + */ +class Coeff +{ + +public: + + Coeff() : coeff_(1.0), span_(), hasSpan_(false) {} + + Coeff(scalar value) : coeff_(value), span_(), hasSpan_(false) {} + + Coeff(scalar coeff, const Field& field) + : coeff_(coeff), span_(field.span()), hasSpan_(true) + {} + + Coeff(const Field& field) : coeff_(1.0), span_(field.span()), hasSpan_(true) {} + + KOKKOS_INLINE_FUNCTION + scalar operator[](const size_t i) const { return (hasSpan_) ? span_[i] * coeff_ : coeff_; } + + bool hasSpan() { return hasSpan_; } + + std::span span() { return span_; } + + Coeff& operator*=(scalar rhs) + { + coeff_ *= rhs; + return *this; + } + + Coeff& operator*=(const Coeff& rhs) + { + if (hasSpan_ && rhs.hasSpan_) + { + NF_ERROR_EXIT("Not implemented"); + } + + if (!hasSpan_ && rhs.hasSpan_) + { + // Take over the span + span_ = rhs.span_; + hasSpan_ = true; + } + + return this->operator*=(rhs.coeff_); + } + + +private: + + scalar coeff_; + + std::span span_; + + bool hasSpan_; +}; + + +namespace detail + +{ +/* @brief function to force evaluation to a field, the field will be resized to hold either a + * single value or the full field + * + * @param field to store the result + */ +void toField(Coeff& coeff, Field& rhs) +{ + if (coeff.hasSpan()) + { + rhs.resize(coeff.span().size()); + fill(rhs, 1.0); + auto rhsSpan = rhs.span(); + // otherwise we are unable to capture values in the lambda + parallelFor( + rhs.exec(), rhs.range(), KOKKOS_LAMBDA(const size_t i) { rhsSpan[i] *= coeff[i]; } + ); + } + else + { + rhs.resize(1); + fill(rhs, coeff[0]); + } +} + +} + +inline Coeff operator*(const Coeff& lhs, const Coeff& rhs) +{ + Coeff result = lhs; + result *= rhs; + return result; +} + +} // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/ddt.hpp b/include/NeoFOAM/dsl/ddt.hpp new file mode 100644 index 000000000..fda1de19c --- /dev/null +++ b/include/NeoFOAM/dsl/ddt.hpp @@ -0,0 +1,43 @@ +// SPDX-License-Identifier: MIT +// +// SPDX-FileCopyrightText: 2023 NeoFOAM authors + +#pragma once + +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/finiteVolume/cellCentred.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" + +namespace NeoFOAM::dsl::temporal +{ + +template +class Ddt : public OperatorMixin +{ + +public: + + Ddt(FieldType& field) : OperatorMixin(field.exec(), field, Operator::Type::Temporal) + {} + + std::string getName() const { return "TimeOperator"; } + + void explicitOperation([[maybe_unused]] Field& source, [[maybe_unused]] scalar scale) + { + NF_ERROR_EXIT("Not implemented"); + } + + void implicitOperation([[maybe_unused]] Field& phi) + { + NF_ERROR_EXIT("Not implemented"); + } +}; + +/* @brief factory function to create a Ddt term as ddt() */ +template +Ddt ddt(FieldType& in) +{ + return Ddt(in); +}; + +} // namespace NeoFOAM diff --git a/include/NeoFOAM/dsl/expression.hpp b/include/NeoFOAM/dsl/expression.hpp new file mode 100644 index 000000000..ac2e2ef68 --- /dev/null +++ b/include/NeoFOAM/dsl/expression.hpp @@ -0,0 +1,167 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#pragma once + +#include +#include +#include + +#include "NeoFOAM/core/primitives/scalar.hpp" +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/dsl/operator.hpp" +#include "NeoFOAM/core/error.hpp" + +namespace NeoFOAM::dsl +{ + + +class Expression +{ +public: + + Expression(const Executor& exec) + : exec_(exec), temporalOperators_(), implicitOperators_(), explicitOperators_() + {} + + /* @brief perform all explicit operation and accumulate the result */ + Field explicitOperation(size_t nCells) const + { + Field source(exec_, nCells, 0.0); + return explicitOperation(source); + } + + /* @brief perform all explicit operation and accumulate the result */ + Field explicitOperation(Field& source) const + { + for (auto& Operator : explicitOperators_) + { + Operator.explicitOperation(source); + } + return source; + } + + void addOperator(const Operator& Operator) + { + switch (Operator.getType()) + { + case Operator::Type::Temporal: + temporalOperators_.push_back(Operator); + break; + case Operator::Type::Implicit: + implicitOperators_.push_back(Operator); + break; + case Operator::Type::Explicit: + explicitOperators_.push_back(Operator); + break; + } + } + + void addExpression(const Expression& equation) + { + for (auto& Operator : equation.temporalOperators_) + { + temporalOperators_.push_back(Operator); + } + for (auto& Operator : equation.implicitOperators_) + { + implicitOperators_.push_back(Operator); + } + for (auto& Operator : equation.explicitOperators_) + { + explicitOperators_.push_back(Operator); + } + } + + + /* @brief getter for the total number of terms in the equation */ + size_t size() const + { + return temporalOperators_.size() + implicitOperators_.size() + explicitOperators_.size(); + } + + // getters + const std::vector& temporalOperators() const { return temporalOperators_; } + + const std::vector& implicitOperators() const { return implicitOperators_; } + + const std::vector& explicitOperators() const { return explicitOperators_; } + + std::vector& temporalOperators() { return temporalOperators_; } + + std::vector& implicitOperators() { return implicitOperators_; } + + std::vector& explicitOperators() { return explicitOperators_; } + + const Executor& exec() const { return exec_; } + +private: + + const Executor exec_; + + std::vector temporalOperators_; + + std::vector implicitOperators_; + + std::vector explicitOperators_; +}; + +[[nodiscard]] Expression operator+(Expression lhs, const Expression& rhs) +{ + lhs.addExpression(rhs); + return lhs; +} + +[[nodiscard]] Expression operator+(Expression lhs, const Operator& rhs) +{ + lhs.addOperator(rhs); + return lhs; +} + +[[nodiscard]] Expression operator+(const Operator& lhs, const Operator& rhs) +{ + Expression expr(lhs.exec()); + expr.addOperator(lhs); + expr.addOperator(rhs); + return expr; +} + +[[nodiscard]] Expression operator*(scalar scale, const Expression& es) +{ + Expression expr(es.exec()); + for (const auto& Operator : es.temporalOperators()) + { + expr.addOperator(scale * Operator); + } + for (const auto& Operator : es.implicitOperators()) + { + expr.addOperator(scale * Operator); + } + for (const auto& Operator : es.explicitOperators()) + { + expr.addOperator(scale * Operator); + } + return expr; +} + +[[nodiscard]] Expression operator-(Expression lhs, const Expression& rhs) +{ + lhs.addExpression(-1.0 * rhs); + return lhs; +} + +[[nodiscard]] Expression operator-(Expression lhs, const Operator& rhs) +{ + lhs.addOperator(-1.0 * rhs); + return lhs; +} + +[[nodiscard]] Expression operator-(const Operator& lhs, const Operator& rhs) +{ + Expression expr(lhs.exec()); + expr.addOperator(lhs); + expr.addOperator(Coeff(-1) * rhs); + return expr; +} + + +} // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/dsl/operator.hpp b/include/NeoFOAM/dsl/operator.hpp new file mode 100644 index 000000000..cd19db332 --- /dev/null +++ b/include/NeoFOAM/dsl/operator.hpp @@ -0,0 +1,250 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#pragma once + +#include +#include +#include +#include +#include + +#include "NeoFOAM/core/primitives/scalar.hpp" +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/core/input.hpp" +#include "NeoFOAM/dsl/coeff.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" + +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 +}; + +/* @class Operator + * @brief A class to represent a 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) : model_ {eqnOperator.model_->clone()} {} + + Operator(Operator&& eqnOperator) : model_ {std::move(eqnOperator.model_)} {} + + void explicitOperation(Field& source) const { model_->explicitOperation(source); } + + void temporalOperation(Field& field) { model_->temporalOperation(field); } + + /* returns the fundamental type of an operator, ie explicit, implicit, temporal */ + Operator::Type getType() const { return model_->getType(); } + + Coeff& getCoefficient() { return model_->getCoefficient(); } + + Coeff getCoefficient() const { return model_->getCoefficient(); } + + /* @brief Given an input this function reads required coeffs */ + void build(const Input& input) { model_->build(input); } + + /* @brief Get the executor */ + const Executor& exec() const { return model_->exec(); } + + +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) const = 0; + + virtual void temporalOperation(Field& field) = 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) const override + { + if constexpr (HasExplicitOperator) + { + concreteOp_.explicitOperation(source); + } + } + + virtual void temporalOperation(Field& field) override + { + if constexpr (HasTemporalOperator) + { + concreteOp_.temporalOperation(field); + } + } + + /* @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_; +}; + + +auto operator*(scalar scalarCoeff, const Operator& rhs) +{ + Operator result = rhs; + result.getCoefficient() *= scalarCoeff; + return result; +} + +auto operator*(const Field& coeffField, const Operator& rhs) +{ + Operator result = rhs; + result.getCoefficient() *= Coeff(coeffField); + return result; +} + +auto operator*(const Coeff& coeff, const Operator& rhs) +{ + Operator result = rhs; + result.getCoefficient() *= coeff; + return result; +} + +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 new file mode 100644 index 000000000..2f84fe57b --- /dev/null +++ b/include/NeoFOAM/dsl/solver.hpp @@ -0,0 +1,55 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#pragma once + +#include +#include +#include +#include +#include + +#include "NeoFOAM/core/primitives/scalar.hpp" +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/core/input.hpp" +#include "NeoFOAM/dsl/expression.hpp" +#include "NeoFOAM/timeIntegration/timeIntegration.hpp" + + +namespace NeoFOAM::dsl +{ + +/* @brief solve an expression + * + * @param exp - Expression which is to be solved/updated. + * @param solution - Solution field, where the solution will be 'written to'. + * @param dt - time step for the temporal integration + * @param fvSchemes - Dictionary containing spatial operator and time integration properties + * @param fvSolution - Dictionary containing linear solver properties + */ +template +void solve( + Expression& exp, + FieldType& solution, + scalar dt, + [[maybe_unused]] const Dictionary& fvSchemes, + [[maybe_unused]] const Dictionary& fvSolution +) +{ + if (exp.temporalOperators().size() == 0 && exp.implicitOperators().size() == 0) + { + NF_ERROR_EXIT("No temporal or implicit terms to solve."); + } + if (exp.temporalOperators().size() > 0) + { + // integrate equations in time + TimeIntegration timeIntegrator(fvSchemes.subDict("ddtSchemes")); + timeIntegrator.solve(exp, solution, dt); + } + else + { + NF_ERROR_EXIT("Not implemented."); + // solve sparse matrix system + } +} + +} // namespace NeoFOAM::dsl diff --git a/include/NeoFOAM/fields/domainField.hpp b/include/NeoFOAM/fields/domainField.hpp index a92d2fb08..905f28382 100644 --- a/include/NeoFOAM/fields/domainField.hpp +++ b/include/NeoFOAM/fields/domainField.hpp @@ -55,7 +55,7 @@ class DomainField {} - DomainField(const ::NeoFOAM::DomainField& rhs) + DomainField(const DomainField& rhs) : exec_(rhs.exec_), internalField_(rhs.internalField_), boundaryFields_(rhs.boundaryFields_) {} @@ -66,7 +66,7 @@ class DomainField {} - DomainField& operator=(const ::NeoFOAM::DomainField& rhs) + DomainField& operator=(const DomainField& rhs) { internalField_ = rhs.internalField_; boundaryFields_ = rhs.boundaryFields_; diff --git a/include/NeoFOAM/fields/field.hpp b/include/NeoFOAM/fields/field.hpp index c5e02ab51..fd5ffd691 100644 --- a/include/NeoFOAM/fields/field.hpp +++ b/include/NeoFOAM/fields/field.hpp @@ -79,7 +79,7 @@ class Field { void* ptr = nullptr; std::visit( - [this, &ptr, size](const auto& concreteExec) + [&ptr, size](const auto& concreteExec) { ptr = concreteExec.alloc(size * sizeof(ValueType)); }, exec_ ); @@ -413,7 +413,13 @@ class Field void validateOtherField(const Field& rhs) const { NF_DEBUG_ASSERT(size() == rhs.size(), "Fields are not the same size."); - NF_DEBUG_ASSERT(exec() == rhs.exec(), "Executors are not the same."); + + std::string execName = std::visit([](auto e) { return e.print(); }, exec()); + std::string rhsExecName = std::visit([](auto e) { return e.print(); }, rhs.exec()); + NF_DEBUG_ASSERT( + exec() == rhs.exec(), + "Executors: " + execName + " and " + rhsExecName + " are not the same." + ); } }; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/calculated.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/calculated.hpp index c59d0c7bd..e508f5e9a 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/calculated.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/calculated.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/core.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/surfaceBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/empty.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/empty.hpp index 9220ad33b..3e327396d 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/empty.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/empty.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/core.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/surfaceBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/fixedValue.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/fixedValue.hpp index 8cca98fdf..b071ef7a9 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/fixedValue.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/surface/fixedValue.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/surfaceBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/calculated.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/calculated.hpp index b85298ec8..497f77db0 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/calculated.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/calculated.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/core.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/volumeBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/empty.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/empty.hpp index c156310c7..4479a3374 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/empty.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/empty.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/core.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/volumeBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedGradient.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedGradient.hpp index f1b1d61cd..4b299d591 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedGradient.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedGradient.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/finiteVolume/cellCentred/boundary/volumeBoundaryFactory.hpp" #include "NeoFOAM/mesh/unstructured.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedValue.hpp b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedValue.hpp index f374b80ef..65580f4c2 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedValue.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/boundary/volume/fixedValue.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/finiteVolume/cellCentred/boundary/volumeBoundaryFactory.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp index 53b68a47d..7e22c8fb1 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp @@ -33,14 +33,38 @@ class GeometricFieldMixin * * @param exec The executor object. * @param mesh The unstructured mesh object. - * @param field The domain field object. + * @param domainField The domain field object. */ GeometricFieldMixin( - const Executor& exec, const UnstructuredMesh& mesh, const DomainField& field + const Executor& exec, + const UnstructuredMesh& mesh, + const DomainField& domainField ) - : exec_(exec), mesh_(mesh), field_(field) + : exec_(exec), mesh_(mesh), field_(domainField) {} + /** + * @brief Constructor for GeometricFieldMixin. + * + * @param exec The executor object. + * @param mesh The unstructured mesh object. + * @param internalField The internal field object. + * @param boundaryFields The boundary field object. + */ + GeometricFieldMixin( + const Executor& exec, + const UnstructuredMesh& mesh, + const Field& internalField, + const BoundaryFields& boundaryFields + ) + : exec_(exec), mesh_(mesh), field_({exec, internalField, boundaryFields}) + { + if (mesh.nCells() != internalField.size()) + { + NF_ERROR_EXIT("Inconsistent size of mesh and internal field detected"); + } + } + /** * @brief Returns a const reference to the internal field. * @@ -55,6 +79,13 @@ class GeometricFieldMixin */ Field& internalField() { return field_.internalField(); } + /** + * @brief Returns the size of the internal field + * + * @return The size of the internal field + */ + size_t size() const { return field_.internalField().size(); } + /** * @brief Returns a const reference to the boundary field. * diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp index 64ad685b7..200bbd18c 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/surfaceField.hpp @@ -55,13 +55,30 @@ class SurfaceField : public GeometricFieldMixin SurfaceField( const Executor& exec, const UnstructuredMesh& mesh, - const Field& internalField, + const DomainField& domainField, const std::vector>& boundaryConditions ) - : GeometricFieldMixin(exec, mesh, {exec, mesh, internalField}), + : GeometricFieldMixin(exec, mesh, domainField), boundaryConditions_(boundaryConditions) {} + /* @brief Constructor for a surfaceField with a given internal field + * + * @param exec The executor + * @param mesh The underlying mesh + * @param internalField the underlying internal field + * @param boundaryConditions a vector of boundary conditions + */ + SurfaceField( + const Executor& exec, + const UnstructuredMesh& mesh, + const Field& internalField, + const BoundaryFields& boundaryFields, + const std::vector>& boundaryConditions + ) + : GeometricFieldMixin(exec, mesh, {exec, mesh, internalField, boundaryFields}), + boundaryConditions_(boundaryConditions) + {} SurfaceField(const SurfaceField& other) : GeometricFieldMixin(other), boundaryConditions_(other.boundaryConditions_) diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp index ad590f425..3e834d470 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp @@ -47,18 +47,37 @@ class VolumeField : public GeometricFieldMixin {} /* @brief Constructor for a VolumeField with a given internal field + * + * @param mesh The underlying mesh + * @param domainField the underlying domain field i.e. combination of internal and boundary + * fields + * @param boundaryConditions a vector of boundary conditions + */ + VolumeField( + const Executor& exec, + const UnstructuredMesh& mesh, + const DomainField& domainField, + const std::vector>& boundaryConditions + ) + : GeometricFieldMixin(exec, mesh, domainField), + boundaryConditions_(boundaryConditions) + {} + + /* @brief Constructor for a VolumeField with a given internal and boundary field * * @param mesh The underlying mesh * @param internalField the underlying internal field + * @param boundaryFields the underlying boundary data fields * @param boundaryConditions a vector of boundary conditions */ VolumeField( const Executor& exec, const UnstructuredMesh& mesh, const Field& internalField, + const BoundaryFields& boundaryFields, const std::vector>& boundaryConditions ) - : GeometricFieldMixin(exec, mesh, {exec, mesh, internalField}), + : GeometricFieldMixin(exec, mesh, internalField, boundaryFields), boundaryConditions_(boundaryConditions) {} diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp index cbcb78e6e..6f75c317c 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp @@ -6,10 +6,10 @@ #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/core/executor/executor.hpp" #include "NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp" -#include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/stencil/geometryScheme.hpp" +#include "NeoFOAM/mesh/unstructured.hpp" -#include "Kokkos_Core.hpp" +#include #include diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp index 12a3d1ed1..d47b35cd6 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/surfaceInterpolation.hpp @@ -8,7 +8,7 @@ #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred.hpp" -#include "Kokkos_Core.hpp" +#include #include diff --git a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp index 3cdf42ac8..469e24c7d 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/interpolation/upwind.hpp @@ -9,7 +9,7 @@ #include "NeoFOAM/mesh/unstructured.hpp" #include "NeoFOAM/finiteVolume/cellCentred/stencil/geometryScheme.hpp" -#include "Kokkos_Core.hpp" +#include #include diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp index ce8dd5f68..203aee628 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp @@ -3,7 +3,7 @@ #pragma once -#include "Kokkos_Core.hpp" +#include #include @@ -27,6 +27,8 @@ class GaussGreenDiv div(VolumeField& divPhi, const SurfaceField& faceFlux, VolumeField& phi ); + void div(Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi); + private: SurfaceInterpolation surfaceInterpolation_; diff --git a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenGrad.hpp b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenGrad.hpp index 43629ec75..ae0c7a07d 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenGrad.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenGrad.hpp @@ -3,9 +3,7 @@ #pragma once -#include - -#include "Kokkos_Core.hpp" +#include #include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/core/executor/executor.hpp" diff --git a/include/NeoFOAM/finiteVolume/cellCentred/stencil/basicGeometryScheme.hpp b/include/NeoFOAM/finiteVolume/cellCentred/stencil/basicGeometryScheme.hpp index 640ca4d95..cb57879e0 100644 --- a/include/NeoFOAM/finiteVolume/cellCentred/stencil/basicGeometryScheme.hpp +++ b/include/NeoFOAM/finiteVolume/cellCentred/stencil/basicGeometryScheme.hpp @@ -3,11 +3,8 @@ #pragma once -#include "NeoFOAM/core/primitives/vector.hpp" #include "NeoFOAM/core/primitives/scalar.hpp" -#include "NeoFOAM/core/primitives/label.hpp" #include "NeoFOAM/core/executor/executor.hpp" -#include "NeoFOAM/fields/field.hpp" #include "NeoFOAM/finiteVolume/cellCentred/stencil/geometryScheme.hpp" namespace NeoFOAM::finiteVolume::cellCentred diff --git a/include/NeoFOAM/mesh/unstructured/boundaryMesh.hpp b/include/NeoFOAM/mesh/unstructured/boundaryMesh.hpp index 2f374e824..ca8a0ce11 100644 --- a/include/NeoFOAM/mesh/unstructured/boundaryMesh.hpp +++ b/include/NeoFOAM/mesh/unstructured/boundaryMesh.hpp @@ -2,7 +2,6 @@ // SPDX-FileCopyrightText: 2023 NeoFOAM authors #pragma once -#include #include #include "NeoFOAM/core/primitives/label.hpp" diff --git a/include/NeoFOAM/mesh/unstructured/unstructuredMesh.hpp b/include/NeoFOAM/mesh/unstructured/unstructuredMesh.hpp index 8717c60d3..25fd2b00a 100644 --- a/include/NeoFOAM/mesh/unstructured/unstructuredMesh.hpp +++ b/include/NeoFOAM/mesh/unstructured/unstructuredMesh.hpp @@ -3,10 +3,6 @@ #pragma once -#include -#include - -#include "NeoFOAM/core/primitives/vector.hpp" #include "NeoFOAM/fields/fieldTypeDefs.hpp" #include "NeoFOAM/mesh/unstructured/boundaryMesh.hpp" #include "NeoFOAM/finiteVolume/cellCentred/stencil/stencilDataBase.hpp" diff --git a/include/NeoFOAM/timeIntegration/forwardEuler.hpp b/include/NeoFOAM/timeIntegration/forwardEuler.hpp new file mode 100644 index 000000000..a02d48e46 --- /dev/null +++ b/include/NeoFOAM/timeIntegration/forwardEuler.hpp @@ -0,0 +1,54 @@ +// SPDX-License-Identifier: MIT +// +// SPDX-FileCopyrightText: 2023 NeoFOAM authors + +#pragma once + +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/timeIntegration/timeIntegration.hpp" + +namespace NeoFOAM::dsl +{ + +template +class ForwardEuler : + public TimeIntegratorBase::template Register> +{ + +public: + + using Base = TimeIntegratorBase::template Register>; + + ForwardEuler(const Dictionary& dict) : Base(dict) {} + + static std::string name() { return "forwardEuler"; } + + static std::string doc() { return "first order time integration method"; } + + static std::string schema() { return "none"; } + + void solve(Expression& eqn, SolutionType& sol, scalar dt) const override + { + auto source = eqn.explicitOperation(sol.size()); + + sol.internalField() -= source * dt; + sol.correctBoundaryConditions(); + + // check if executor is GPU + if (std::holds_alternative(eqn.exec())) + { + Kokkos::fence(); + } + }; + + std::unique_ptr> clone() const override + { + return std::make_unique(*this); + } +}; + +// unfortunately needs explicit instantiation +template class ForwardEuler>; + + +} // namespace NeoFOAM diff --git a/include/NeoFOAM/timeIntegration/timeIntegration.hpp b/include/NeoFOAM/timeIntegration/timeIntegration.hpp new file mode 100644 index 000000000..89b1feff9 --- /dev/null +++ b/include/NeoFOAM/timeIntegration/timeIntegration.hpp @@ -0,0 +1,76 @@ +// SPDX-License-Identifier: MIT +// +// SPDX-FileCopyrightText: 2023 NeoFOAM authors + +#pragma once + +#include + +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/finiteVolume/cellCentred.hpp" +#include "NeoFOAM/dsl/expression.hpp" + +namespace NeoFOAM::dsl +{ + +/* @class Factory class to create time integration method by a given name + * using NeoFOAMs runTimeFactory mechanism + */ +template +class TimeIntegratorBase : + public RuntimeSelectionFactory, Parameters> +{ + +public: + + static std::string name() { return "timeIntegrationFactory"; } + + TimeIntegratorBase(const Dictionary& dict) : dict_(dict) {} + + virtual ~TimeIntegratorBase() {} + + virtual void solve(Expression& eqn, SolutionType& sol, scalar dt) + const = 0; // Pure virtual function for solving + + // Pure virtual function for cloning + virtual std::unique_ptr clone() const = 0; + +protected: + + const Dictionary& dict_; +}; + +/* @class Factory class to create time integration method by a given name + * using NeoFOAMs runTimeFactory mechanism + * + * @tparam SolutionType Type of the solution field eg, volumeField or just a plain Field + */ +template +class TimeIntegration +{ + +public: + + TimeIntegration(const TimeIntegration& timeIntegrator) + : timeIntegratorStrategy_(timeIntegrator.timeIntegratorStrategy_->clone()) {}; + + TimeIntegration(TimeIntegration&& timeIntegrator) + : timeIntegratorStrategy_(std::move(timeIntegrator.timeIntegratorStrategy_)) {}; + + TimeIntegration(const Dictionary& dict) + : timeIntegratorStrategy_( + TimeIntegratorBase::create(dict.get("type"), dict) + ) {}; + + void solve(Expression& eqn, SolutionType& sol, scalar dt) + { + timeIntegratorStrategy_->solve(eqn, sol, dt); + } + +private: + + std::unique_ptr> timeIntegratorStrategy_; +}; + + +} // namespace NeoFOAM::dsl diff --git a/scripts/package.py b/scripts/package.py new file mode 100644 index 000000000..d69f8e89c --- /dev/null +++ b/scripts/package.py @@ -0,0 +1,21 @@ +# Copyright 2013-2024 Lawrence Livermore National Security, LLC and other +# Spack Project Developers. See the top-level COPYRIGHT file for details. +# +# SPDX-License-Identifier: MIT + +from spack.package import * + +class Neofoam(CMakePackage): + """NeoFOAM is a WIP prototype of a modern CFD core.""" + + homepage = "https://github.com/exasim-project/NeoFOAM" + git = homepage + + # maintainers("github_user1", "github_user2") + + license("UNKNOWN", checked_by="github_user1") + + version("develop", branch="main") + + depends_on("mpi") + depends_on("kokkos@4.3.0") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 34d39e564..8d8136c33 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,5 @@ # SPDX-License-Identifier: Unlicense # SPDX-FileCopyrightText: 2023 NeoFOAM authors -# add_subdirectory(DSL) add_library(NeoFOAM SHARED) @@ -13,7 +12,6 @@ target_sources( "core/dictionary.cpp" "core/demangle.cpp" "core/tokenList.cpp" - "core/kokkos.cpp" "executor/CPUExecutor.cpp" "executor/GPUExecutor.cpp" "executor/serialExecutor.cpp" diff --git a/src/core/demangle.cpp b/src/core/demangle.cpp index 9f605c79e..1d1309ae9 100644 --- a/src/core/demangle.cpp +++ b/src/core/demangle.cpp @@ -1,6 +1,9 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2024 NeoFOAM authors +#include // for __cxa_demangle +#include // for free + #include "NeoFOAM/core/demangle.hpp" diff --git a/src/core/dictionary.cpp b/src/core/dictionary.cpp index e5a8297d7..9ab0124ce 100644 --- a/src/core/dictionary.cpp +++ b/src/core/dictionary.cpp @@ -1,9 +1,11 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors +#include +#include // for operator<<, basic_ostream, endl, cerr, ostream + #include "NeoFOAM/core/dictionary.hpp" #include "NeoFOAM/core/error.hpp" -#include namespace NeoFOAM { @@ -13,6 +15,7 @@ void logOutRange( const std::unordered_map& data ) { + // TODO use NeoFOAM error here std::cerr << "Key not found: " << key << " \n" << "available keys are: \n" << std::accumulate( diff --git a/src/core/kokkos.cpp b/src/core/kokkos.cpp deleted file mode 100644 index 577935344..000000000 --- a/src/core/kokkos.cpp +++ /dev/null @@ -1,16 +0,0 @@ -// SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: 2023 NeoFOAM authors - -#include - -namespace NeoFOAM -{ - -struct HelloWorld -{ - KOKKOS_INLINE_FUNCTION - void operator()(const int i) const { printf("Hello from i = %i\n", i); } -}; - -void foo() { Kokkos::parallel_for("HelloWorld", 15, HelloWorld()); } -} // namespace NeoFOAM diff --git a/src/core/primitives/vector.cpp b/src/core/primitives/vector.cpp index 60b8054eb..7e7bbab28 100644 --- a/src/core/primitives/vector.cpp +++ b/src/core/primitives/vector.cpp @@ -1,9 +1,7 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors - #include "NeoFOAM/core/primitives/vector.hpp" -#include "NeoFOAM/core/primitives/scalar.hpp" namespace NeoFOAM { diff --git a/src/finiteVolume/cellCentred/interpolation/linear.cpp b/src/finiteVolume/cellCentred/interpolation/linear.cpp index eb0c9de64..bc057cd62 100644 --- a/src/finiteVolume/cellCentred/interpolation/linear.cpp +++ b/src/finiteVolume/cellCentred/interpolation/linear.cpp @@ -4,7 +4,6 @@ #include #include "NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp" -#include "NeoFOAM/core/error.hpp" #include "NeoFOAM/core/parallelAlgorithms.hpp" namespace NeoFOAM::finiteVolume::cellCentred diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp index d1193f4bb..b2ac722d3 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp @@ -1,11 +1,8 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors -#include - #include "NeoFOAM/core/parallelAlgorithms.hpp" #include "NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp" -#include "NeoFOAM/finiteVolume/cellCentred/stencil/geometryScheme.hpp" namespace NeoFOAM::finiteVolume::cellCentred { @@ -33,9 +30,85 @@ void computeDiv( size_t nInternalFaces = mesh.nInternalFaces(); const auto surfV = mesh.cellVolumes().span(); + // check if the executor is GPU + if (std::holds_alternative(exec)) + { + for (size_t i = 0; i < nInternalFaces; i++) + { + scalar flux = surfFaceFlux[i] * surfPhif[i]; + surfDivPhi[static_cast(surfOwner[i])] += flux; + surfDivPhi[static_cast(surfNeighbour[i])] -= flux; + } + + for (size_t i = nInternalFaces; i < surfPhif.size(); i++) + { + auto own = static_cast(surfFaceCells[i - nInternalFaces]); + scalar valueOwn = surfFaceFlux[i] * surfPhif[i]; + surfDivPhi[own] += valueOwn; + } + + for (size_t celli = 0; celli < mesh.nCells(); celli++) + { + surfDivPhi[celli] *= 1 / surfV[celli]; + } + } + else + { + parallelFor( + exec, + {0, nInternalFaces}, + KOKKOS_LAMBDA(const size_t i) { + scalar flux = surfFaceFlux[i] * surfPhif[i]; + Kokkos::atomic_add(&surfDivPhi[static_cast(surfOwner[i])], flux); + Kokkos::atomic_sub(&surfDivPhi[static_cast(surfNeighbour[i])], flux); + } + ); + + parallelFor( + exec, + {nInternalFaces, surfPhif.size()}, + KOKKOS_LAMBDA(const size_t i) { + auto own = static_cast(surfFaceCells[i - nInternalFaces]); + scalar valueOwn = surfFaceFlux[i] * surfPhif[i]; + Kokkos::atomic_add(&surfDivPhi[own], valueOwn); + } + ); + + parallelFor( + exec, + {0, mesh.nCells()}, + KOKKOS_LAMBDA(const size_t celli) { surfDivPhi[celli] *= 1 / surfV[celli]; } + ); + } +} + +void computeDiv( + const SurfaceField& faceFlux, + const VolumeField& phi, + [[maybe_unused]] const SurfaceInterpolation& surfInterp, + Field& divPhi +) +{ + const UnstructuredMesh& mesh = phi.mesh(); + const auto exec = phi.exec(); + SurfaceField phif(exec, mesh, createCalculatedBCs(mesh)); + const auto surfFaceCells = mesh.boundaryMesh().faceCells().span(); + + // TODO check if copy can be avoided + auto faceFluxCopy = SurfaceField(faceFlux); + surfInterp.interpolate(phif, phi, faceFluxCopy); + + auto surfDivPhi = divPhi.span(); + + const auto surfPhif = phif.internalField().span(); + const auto surfOwner = mesh.faceOwner().span(); + const auto surfNeighbour = mesh.faceNeighbour().span(); + const auto surfFaceFlux = faceFlux.internalField().span(); + size_t nInternalFaces = mesh.nInternalFaces(); + const auto surfV = mesh.cellVolumes().span(); // check if the executor is GPU - if (std::holds_alternative(exec)) + if (std::holds_alternative(exec)) { for (size_t i = 0; i < nInternalFaces; i++) { @@ -101,4 +174,11 @@ void GaussGreenDiv::div( computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); }; +void GaussGreenDiv::div( + Field& divPhi, const SurfaceField& faceFlux, VolumeField& phi +) +{ + computeDiv(faceFlux, phi, surfaceInterpolation_, divPhi); +}; + }; diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp index f0f8c4451..94bd88873 100644 --- a/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp +++ b/src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp @@ -1,8 +1,6 @@ // SPDX-License-Identifier: MIT // SPDX-FileCopyrightText: 2023 NeoFOAM authors -#include - #include "NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenGrad.hpp" #include "NeoFOAM/finiteVolume/cellCentred/interpolation/linear.hpp" #include "NeoFOAM/core/parallelAlgorithms.hpp" @@ -32,7 +30,7 @@ void computeGrad( const auto surfV = mesh.cellVolumes().span(); - if (std::holds_alternative(exec)) + if (std::holds_alternative(exec)) { for (size_t i = 0; i < nInternalFaces; i++) { diff --git a/src/mesh/unstructured/unstructuredMesh.cpp b/src/mesh/unstructured/unstructuredMesh.cpp index 5e172f3fb..d7b26410c 100644 --- a/src/mesh/unstructured/unstructuredMesh.cpp +++ b/src/mesh/unstructured/unstructuredMesh.cpp @@ -3,6 +3,9 @@ #include "NeoFOAM/mesh/unstructured/unstructuredMesh.hpp" +#include "NeoFOAM/core/primitives/vector.hpp" // for Vector + + namespace NeoFOAM { @@ -90,7 +93,7 @@ UnstructuredMesh createSingleCellMesh(const Executor exec) ); return UnstructuredMesh( {exec, {{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 0, 0}}}, // points, - {exec, {1}}, // cellVolumes + {exec, 1}, // cellVolumes {exec, {{0.5, 0.5, 0.0}}}, // cellCentres faceAreasVectors, faceCentresVectors, diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6a7629b87..a0dbcffef 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -6,7 +6,8 @@ list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras) # has to be first since it adds the main target add_subdirectory(catch2) -# This function creates unit tests. It allows the following keywords: +# include "NeoFOAM/core/parallelAlgorithms.hpp" This function creates unit tests. It provides the +# following keywords: # # * MPI_SIZE: the number of MPI processors to be used, defaults to 1 if not set # * COMMAND: the test command (same behavior as for CMake's add_test), defaults to the test name @@ -74,3 +75,4 @@ add_subdirectory(core) add_subdirectory(fields) add_subdirectory(finiteVolume) add_subdirectory(mesh) +add_subdirectory(dsl) diff --git a/test/catch2/test_main.cpp b/test/catch2/test_main.cpp index 35f25e599..304ac404d 100644 --- a/test/catch2/test_main.cpp +++ b/test/catch2/test_main.cpp @@ -7,7 +7,7 @@ #include "catch2/catch_test_macros.hpp" #include "catch2/generators/catch_generators_adapters.hpp" #include "catch2/reporters/catch_reporter_registrars.hpp" -#include "Kokkos_Core.hpp" +#include int main(int argc, char* argv[]) diff --git a/test/core/CMakeLists.txt b/test/core/CMakeLists.txt index e8bb2c689..55de3363d 100644 --- a/test/core/CMakeLists.txt +++ b/test/core/CMakeLists.txt @@ -9,7 +9,7 @@ add_subdirectory(primitives) neofoam_unit_test(dictionary) neofoam_unit_test(tokenList) -neofoam_unit_test(inputs) +neofoam_unit_test(input) neofoam_unit_test(executor) neofoam_unit_test(parallelAlgorithms) diff --git a/test/core/inputs.cpp b/test/core/input.cpp similarity index 98% rename from test/core/inputs.cpp rename to test/core/input.cpp index 5d9a86e54..2bef240fb 100644 --- a/test/core/inputs.cpp +++ b/test/core/input.cpp @@ -7,7 +7,7 @@ #include #include -#include "NeoFOAM/core/inputs.hpp" +#include "NeoFOAM/core/input.hpp" #include "NeoFOAM/core/primitives/label.hpp" #include "NeoFOAM/core/primitives/scalar.hpp" diff --git a/test/dsl/CMakeLists.txt b/test/dsl/CMakeLists.txt new file mode 100644 index 000000000..c0b7c92f2 --- /dev/null +++ b/test/dsl/CMakeLists.txt @@ -0,0 +1,8 @@ +# SPDX-License-Identifier: Unlicense +# SPDX-FileCopyrightText: 2024 NeoFOAM authors + +neofoam_unit_test(coeff) +neofoam_unit_test(operator) +neofoam_unit_test(expression) + +add_subdirectory(timeIntegration) diff --git a/test/dsl/coeff.cpp b/test/dsl/coeff.cpp new file mode 100644 index 000000000..f945a00ae --- /dev/null +++ b/test/dsl/coeff.cpp @@ -0,0 +1,112 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create + // a custom main +#include +#include +#include +#include + +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/dsl/coeff.hpp" + +using Field = NeoFOAM::Field; +using Coeff = NeoFOAM::dsl::Coeff; +using namespace NeoFOAM::dsl; + + +TEST_CASE("Coeff") +{ + 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.print(); }, exec); + + SECTION("Coefficient evaluation on " + execName) + { + Field fA(exec, 3, 2.0); + Field res(exec, 1); + + Coeff a {}; + Coeff b {2.0}; + Coeff c = 2 * a * b; + REQUIRE(c[0] == 4.0); + + Coeff d {3.0, fA}; + detail::toField(d, res); + auto hostResD = res.copyToHost(); + REQUIRE(hostResD.data()[0] == 6.0); + REQUIRE(hostResD.data()[1] == 6.0); + REQUIRE(hostResD.data()[2] == 6.0); + + Coeff e = d * b; + detail::toField(e, res); + auto hostResE = res.copyToHost(); + REQUIRE(hostResE.data()[0] == 12.0); + REQUIRE(hostResE.data()[1] == 12.0); + REQUIRE(hostResE.data()[2] == 12.0); + + Coeff f = b * d; + detail::toField(f, res); + auto hostResF = res.copyToHost(); + REQUIRE(hostResF.data()[0] == 12.0); + REQUIRE(hostResF.data()[1] == 12.0); + REQUIRE(hostResF.data()[2] == 12.0); + } + + SECTION("evaluation in parallelFor" + execName) + { + size_t size = 3; + + Field fieldA(exec, size, 0.0); + Field fieldB(exec, size, 1.0); + + SECTION("span") + { + Coeff coeff = fieldB; // is a span with uniform value 1.0 + { + NeoFOAM::parallelFor( + fieldA, KOKKOS_LAMBDA(const size_t i) { return coeff[i] + 2.0; } + ); + }; + auto hostFieldA = fieldA.copyToHost(); + REQUIRE(coeff.hasSpan() == true); + REQUIRE(hostFieldA[0] == 3.0); + REQUIRE(hostFieldA[1] == 3.0); + REQUIRE(hostFieldA[2] == 3.0); + } + + SECTION("scalar") + { + Coeff coeff = Coeff(2.0); + { + NeoFOAM::parallelFor( + fieldA, KOKKOS_LAMBDA(const size_t i) { return coeff[i] + 2.0; } + ); + }; + auto hostFieldA = fieldA.copyToHost(); + REQUIRE(coeff.hasSpan() == false); + REQUIRE(hostFieldA[0] == 4.0); + REQUIRE(hostFieldA[1] == 4.0); + REQUIRE(hostFieldA[2] == 4.0); + } + + SECTION("span and scalar") + { + Coeff coeff {-5.0, fieldB}; + { + NeoFOAM::parallelFor( + fieldA, KOKKOS_LAMBDA(const size_t i) { return coeff[i] + 2.0; } + ); + }; + auto hostFieldA = fieldA.copyToHost(); + REQUIRE(coeff.hasSpan() == true); + REQUIRE(hostFieldA[0] == -3.0); + REQUIRE(hostFieldA[1] == -3.0); + REQUIRE(hostFieldA[2] == -3.0); + } + } +} diff --git a/test/dsl/common.hpp b/test/dsl/common.hpp new file mode 100644 index 000000000..7a7115905 --- /dev/null +++ b/test/dsl/common.hpp @@ -0,0 +1,53 @@ +// SPDX-License-Identifier: MIT +// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors +#define CATCH_CONFIG_RUNNER // Define this before including catch.hpp to create + // a custom main +#include +#include +#include +#include + +#include "NeoFOAM/fields/field.hpp" +#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp" +#include "NeoFOAM/dsl/coeff.hpp" +#include "NeoFOAM/dsl/operator.hpp" + +namespace fvcc = NeoFOAM::finiteVolume::cellCentred; + +using Field = NeoFOAM::Field; +using Coeff = NeoFOAM::dsl::Coeff; +using Operator = NeoFOAM::dsl::Operator; +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 */ +class Dummy : public OperatorMixin +{ + +public: + + Dummy(VolumeField& field) : OperatorMixin(field.exec(), field, Operator::Type::Explicit) {} + + void explicitOperation(Field& source) const + { + 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]; } + ); + } + + std::string getName() const { return "Dummy"; } +}; + +NeoFOAM::scalar getField(const NeoFOAM::Field& source) +{ + auto sourceField = source.copyToHost(); + return sourceField.span()[0]; +} diff --git a/test/dsl/expression.cpp b/test/dsl/expression.cpp new file mode 100644 index 000000000..c1d509e36 --- /dev/null +++ b/test/dsl/expression.cpp @@ -0,0 +1,52 @@ +// 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" +#include "NeoFOAM/dsl/expression.hpp" + +using Expression = NeoFOAM::dsl::Expression; + +TEST_CASE("Expression") +{ + 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.print(); }, exec); + auto mesh = NeoFOAM::createSingleCellMesh(exec); + + const size_t size {1}; + Field fA(exec, size, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + + std::vector> bcs {}; + auto vf = VolumeField(exec, mesh, fA, bf, bcs); + auto fB = Field(exec, 1, 4.0); + + auto a = Dummy(vf); + auto b = Dummy(vf); + + SECTION("Create equation and perform explicitOperation on " + execName) + { + auto eqnA = a + b; + auto eqnB = fB * Dummy(vf) + 2 * Dummy(vf); + auto eqnC = Expression(2 * a - b); + auto eqnD = 3 * (2 * a - b); + auto eqnE = (2 * a - b) + (2 * a - b); + auto eqnF = (2 * a - b) - (2 * a - b); + + REQUIRE(eqnA.size() == 2); + REQUIRE(eqnB.size() == 2); + REQUIRE(eqnC.size() == 2); + + REQUIRE(getField(eqnA.explicitOperation(size)) == 4); + REQUIRE(getField(eqnB.explicitOperation(size)) == 12); + REQUIRE(getField(eqnC.explicitOperation(size)) == 2); + REQUIRE(getField(eqnD.explicitOperation(size)) == 6); + REQUIRE(getField(eqnE.explicitOperation(size)) == 4); + REQUIRE(getField(eqnF.explicitOperation(size)) == 0); + } +} diff --git a/test/dsl/operator.cpp b/test/dsl/operator.cpp new file mode 100644 index 000000000..5203aa9c9 --- /dev/null +++ b/test/dsl/operator.cpp @@ -0,0 +1,66 @@ +// 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" + +TEST_CASE("Operator") +{ + 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.print(); }, 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, mesh, fA, bf, bcs); + auto b = Dummy(vf); + + REQUIRE(b.getName() == "Dummy"); + REQUIRE(b.getType() == Operator::Type::Explicit); + } + + SECTION("Supports Coefficients" + execName) + { + std::vector> bcs {}; + + Field fA(exec, 1, 2.0); + Field fB(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + auto vf = VolumeField(exec, mesh, fA, bf, bcs); + + auto c = 2 * Dummy(vf); + auto d = fB * Dummy(vf); + auto e = Coeff(-3, fB) * Dummy(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); + + // 2 += 2 * 2 + auto hostSourceC = source.copyToHost(); + REQUIRE(hostSourceC.span()[0] == 6.0); + + // 6 += 2 * 2 + d.explicitOperation(source); + auto hostSourceD = source.copyToHost(); + REQUIRE(hostSourceD.span()[0] == 10.0); + + // 10 += - 6 * 2 + e.explicitOperation(source); + auto hostSourceE = source.copyToHost(); + REQUIRE(hostSourceE.span()[0] == -2.0); + } +} diff --git a/test/dsl/timeIntegration/CMakeLists.txt b/test/dsl/timeIntegration/CMakeLists.txt new file mode 100644 index 000000000..fd05e4245 --- /dev/null +++ b/test/dsl/timeIntegration/CMakeLists.txt @@ -0,0 +1,4 @@ +# SPDX-License-Identifier: Unlicense +# SPDX-FileCopyrightText: 2024 NeoFOAM authors + +neofoam_unit_test(timeIntegration) diff --git a/test/dsl/timeIntegration/timeIntegration.cpp b/test/dsl/timeIntegration/timeIntegration.cpp new file mode 100644 index 000000000..4abb44e66 --- /dev/null +++ b/test/dsl/timeIntegration/timeIntegration.cpp @@ -0,0 +1,55 @@ +// 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 "../common.hpp" + +#include "NeoFOAM/core/dictionary.hpp" +#include "NeoFOAM/core/parallelAlgorithms.hpp" +#include "NeoFOAM/timeIntegration/forwardEuler.hpp" +#include "NeoFOAM/dsl/solver.hpp" +#include "NeoFOAM/dsl/ddt.hpp" + +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.print(); }, exec); + auto mesh = NeoFOAM::createSingleCellMesh(exec); + + NeoFOAM::Dictionary fvSchemes; + NeoFOAM::Dictionary ddtSchemes; + ddtSchemes.insert("type", std::string("forwardEuler")); + fvSchemes.insert("ddtSchemes", ddtSchemes); + NeoFOAM::Dictionary fvSolution; + + Field fA(exec, 1, 2.0); + BoundaryFields bf(exec, mesh.nBoundaryFaces(), mesh.nBoundaries()); + + std::vector> bcs {}; + auto vf = VolumeField(exec, mesh, fA, bf, bcs); + + 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}; + + // 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 + solve(eqn, vf, dt, fvSchemes, fvSolution); + REQUIRE(getField(vf.internalField()) == -2.0); + } +} diff --git a/test/fields/field.cpp b/test/fields/field.cpp index b157aea32..9afde7241 100644 --- a/test/fields/field.cpp +++ b/test/fields/field.cpp @@ -274,17 +274,12 @@ TEST_CASE("getSpans") NeoFOAM::Executor(NeoFOAM::GPUExecutor {}) ); - NeoFOAM::Field a(exec, 3, 1.0); NeoFOAM::Field b(exec, 3, 2.0); NeoFOAM::Field c(exec, 3, 3.0); auto [spanA, spanB, spanC] = NeoFOAM::spans(a, b, c); - REQUIRE(spanA[0] == 1.0); - REQUIRE(spanB[0] == 2.0); - REQUIRE(spanC[0] == 3.0); - NeoFOAM::parallelFor( a, KOKKOS_LAMBDA(const NeoFOAM::size_t i) { return spanB[i] + spanC[i]; } );