Merge pull request #259 from exasim-project/dsl/temporalOperators

Dsl/temporal operators

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)

include/NeoFOAM/core/parallelAlgorithms.hpp

@@ -124,7 +124,7 @@ void parallelReduce(
     const NeoFOAM::Executor& exec, std::pair<size_t, size_t> 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<typename ValueType, typename Kernel, typename T>
 void parallelReduce(Field<ValueType>& 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<typename Executor, typename Kernel>
@@ -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);
 }
 
 

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"

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<FieldType>
 
 public:
 
-    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Temporal)
+    Ddt(FieldType& field) : OperatorMixin<FieldType>(field.exec(), field, Operator::Type::Implicit)
     {}
 
     std::string getName() const { return "TimeOperator"; }
 
-    void explicitOperation([[maybe_unused]] Field<scalar>& source, [[maybe_unused]] scalar scale)
+    void explicitOperation(
+        [[maybe_unused]] Field<scalar>& source,
+        [[maybe_unused]] scalar t,
+        [[maybe_unused]] scalar dt
+    )
     {
         NF_ERROR_EXIT("Not implemented");
     }
 
-    void implicitOperation([[maybe_unused]] Field<scalar>& phi)
+    void implicitOperation( [[maybe_unused]] la::LinearSystem<scalar, localIdx>& ls,[[maybe_unused]] scalar t, [[maybe_unused]] scalar dt)
     {
         NF_ERROR_EXIT("Not implemented");
     }

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<NeoFOAM::scalar>& faceFlux, fvcc::VolumeField<NeoFOAM::scalar>& 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<NeoFOAM::scalar>& coeff, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Explicit, coeff, phi));
 }
 
 

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,91 +52,92 @@ class Expression
     /* @brief perform all explicit operation and accumulate the result */
     Field<scalar> explicitOperation(Field<scalar>& 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<scalar, localIdx> implicitOperation()
+    Field<scalar> explicitOperation(Field<scalar>& 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<scalar, localIdx> 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<scalar, localIdx>& 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<Operator>& temporalOperators() const { return temporalOperators_; }
-
-    const std::vector<Operator>& implicitOperators() const { return implicitOperators_; }
+    const std::vector<TemporalOperator>& temporalOperators() const { return temporalOperators_; }
 
-    const std::vector<Operator>& explicitOperators() const { return explicitOperators_; }
+    const std::vector<SpatialOperator>& spatialOperators() const { return spatialOperators_; }
 
-    std::vector<Operator>& temporalOperators() { return temporalOperators_; }
+    std::vector<TemporalOperator>& temporalOperators() { return temporalOperators_; }
 
-    std::vector<Operator>& implicitOperators() { return implicitOperators_; }
-
-    std::vector<Operator>& explicitOperators() { return explicitOperators_; }
+    std::vector<SpatialOperator>& spatialOperators() { return spatialOperators_; }
 
     const Executor& exec() const { return exec_; }
 
 private:
 
     const Executor exec_;
 
-    std::vector<Operator> temporalOperators_;
-
-    std::vector<Operator> implicitOperators_;
+    std::vector<TemporalOperator> temporalOperators_;
 
-    std::vector<Operator> explicitOperators_;
+    std::vector<SpatialOperator> 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<typename leftOperator, typename rightOperator>
+[[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);

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<NeoFOAM::scalar>& phi) { return dsl::temporal::ddt(phi); }
+TemporalOperator ddt(fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return fvcc::DdtOperator(dsl::Operator::Type::Implicit, phi);
+}
+
+SpatialOperator
+Source(fvcc::VolumeField<NeoFOAM::scalar>& coeff, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return SpatialOperator(fvcc::SourceTerm(dsl::Operator::Type::Implicit, coeff, phi));
+}
+
+SpatialOperator
+div(fvcc::SurfaceField<NeoFOAM::scalar>& faceFlux, fvcc::VolumeField<NeoFOAM::scalar>& phi)
+{
+    return SpatialOperator(fvcc::DivOperator(dsl::Operator::Type::Implicit, faceFlux, phi));
+}
 
 } // namespace NeoFOAM