Implementation of basic time integration

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<scalar>& field)
+        : coeff_(coeff), span_(field.span()), hasSpan_(true)
+    {}
+
+    Coeff(const Field<scalar>& 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<const scalar> 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<const scalar> 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<scalar>& 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

include/NeoFOAM/dsl/ddt.hpp

@@ -0,0 +1,50 @@
+// 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
+{
+
+
+// TODO add free factory function
+template<typename FieldType>
+class Ddt : public OperatorMixin<FieldType>
+{
+
+public:
+
+    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Temporal)
+    {}
+
+    std::string getName() const { return "TimeOperator"; }
+
+    void explicitOperation([[maybe_unused]] Field<scalar>& source, [[maybe_unused]] scalar scale)
+    {
+        NF_ERROR_EXIT("Not implemented");
+    }
+
+    void implicitOperation([[maybe_unused]] Field<scalar>& phi)
+    {
+        NF_ERROR_EXIT("Not implemented");
+    }
+
+private:
+};
+
+// see
+// https://github.com/exasim-project/NeoFOAM/blob/dsl/operatorIntergration/include/NeoFOAM/finiteVolume/cellCentred/operators/explicitOperators/expOp.hpp
+
+template<typename FieldType>
+Ddt<FieldType> ddt(FieldType& in)
+{
+    return Ddt(in);
+};
+
+
+} // namespace NeoFOAM

include/NeoFOAM/dsl/expression.hpp

@@ -0,0 +1,170 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <vector>
+#include <utility>
+
+#include "NeoFOAM/core/primitives/scalar.hpp"
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/dsl/operator.hpp"
+#include "NeoFOAM/core/error.hpp"
+
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.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<scalar> explicitOperation(size_t nCells) const
+    {
+        Field<scalar> source(exec_, nCells, 0.0);
+        return explicitOperation(source);
+    }
+
+    /* @brief perform all explicit operation and accumulate the result */
+    Field<scalar> explicitOperation(Field<scalar>& 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<Operator>& temporalOperators() const { return temporalOperators_; }
+
+    const std::vector<Operator>& implicitOperators() const { return implicitOperators_; }
+
+    const std::vector<Operator>& explicitOperators() const { return explicitOperators_; }
+
+    std::vector<Operator>& temporalOperators() { return temporalOperators_; }
+
+    std::vector<Operator>& implicitOperators() { return implicitOperators_; }
+
+    std::vector<Operator>& explicitOperators() { return explicitOperators_; }
+
+    const Executor& exec() const { return exec_; }
+
+private:
+
+    const Executor exec_;
+
+    std::vector<Operator> temporalOperators_;
+
+    std::vector<Operator> implicitOperators_;
+
+    std::vector<Operator> explicitOperators_;
+};
+
+Expression operator+(Expression lhs, const Expression& rhs)
+{
+    lhs.addExpression(rhs);
+    return lhs;
+}
+
+Expression operator+(Expression lhs, const Operator& rhs)
+{
+    lhs.addOperator(rhs);
+    return lhs;
+}
+
+Expression operator+(const Operator& lhs, const Operator& rhs)
+{
+    Expression expr(lhs.exec());
+    expr.addOperator(lhs);
+    expr.addOperator(rhs);
+    return expr;
+}
+
+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;
+}
+
+Expression operator-(Expression lhs, const Expression& rhs)
+{
+    lhs.addExpression(-1.0 * rhs);
+    return lhs;
+}
+
+Expression operator-(Expression lhs, const Operator& rhs)
+{
+    lhs.addOperator(-1.0 * rhs);
+    return lhs;
+}
+
+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