diff --git a/include/NeoFOAM/DSL/coeff.hpp b/include/NeoFOAM/DSL/coeff.hpp
index b78c6453b..d2b2a6623 100644
--- a/include/NeoFOAM/DSL/coeff.hpp
+++ b/include/NeoFOAM/DSL/coeff.hpp
@@ -43,7 +43,6 @@ class Coeff
         return *this;
     }
 
-
     Coeff& operator*=(const Coeff& rhs)
     {
         if (hasSpan_ && rhs.hasSpan_)
@@ -73,6 +72,7 @@ class Coeff
 
 
 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
@@ -97,8 +97,8 @@ void toField(Coeff& coeff, Field<scalar>& rhs)
         fill(rhs, coeff[0]);
     }
 }
-}
 
+}
 
 inline Coeff operator*(const Coeff& lhs, const Coeff& rhs)
 {
diff --git a/include/NeoFOAM/DSL/eqnSystem.hpp b/include/NeoFOAM/DSL/eqnSystem.hpp
index f7003720f..ae34351e9 100644
--- a/include/NeoFOAM/DSL/eqnSystem.hpp
+++ b/include/NeoFOAM/DSL/eqnSystem.hpp
@@ -9,7 +9,7 @@
 
 #include "NeoFOAM/core/primitives/scalar.hpp"
 #include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/DSL/eqnTerm.hpp"
+#include "NeoFOAM/DSL/operator.hpp"
 #include "NeoFOAM/core/error.hpp"
 
 namespace NeoFOAM::DSL
@@ -20,60 +20,60 @@ class EqnSystem
 public:
 
     EqnSystem(const NeoFOAM::Executor& exec, std::size_t nCells)
-        : exec_(exec), nCells_(nCells), temporalTerms_(), implicitTerms_(), explicitTerms_(),
-          volumeField_(nullptr)
+        : exec_(exec), nCells_(nCells), temporalOperators_(), implicitOperators_(),
+          explicitOperators_(), volumeField_(nullptr)
     {}
 
     NeoFOAM::Field<NeoFOAM::scalar> explicitOperation()
     {
         NeoFOAM::Field<NeoFOAM::scalar> source(exec_, nCells_);
         NeoFOAM::fill(source, 0.0);
-        for (auto& eqnTerm : explicitTerms_)
+        for (auto& Operator : explicitOperators_)
         {
-            eqnTerm.explicitOperation(source);
+            Operator.explicitOperation(source);
         }
         return source;
     }
 
-    void addTerm(const EqnTerm& eqnTerm)
+    void addOperator(const Operator& Operator)
     {
-        switch (eqnTerm.getType())
+        switch (Operator.getType())
         {
-        case EqnTerm::Type::Temporal:
-            temporalTerms_.push_back(eqnTerm);
+        case Operator::Type::Temporal:
+            temporalOperators_.push_back(Operator);
             break;
-        case EqnTerm::Type::Implicit:
-            implicitTerms_.push_back(eqnTerm);
+        case Operator::Type::Implicit:
+            implicitOperators_.push_back(Operator);
             break;
-        case EqnTerm::Type::Explicit:
-            explicitTerms_.push_back(eqnTerm);
+        case Operator::Type::Explicit:
+            explicitOperators_.push_back(Operator);
             break;
         }
     }
 
     void addSystem(const EqnSystem& eqnSys)
     {
-        for (auto& eqnTerm : eqnSys.temporalTerms_)
+        for (auto& Operator : eqnSys.temporalOperators_)
         {
-            temporalTerms_.push_back(eqnTerm);
+            temporalOperators_.push_back(Operator);
         }
-        for (auto& eqnTerm : eqnSys.implicitTerms_)
+        for (auto& Operator : eqnSys.implicitOperators_)
         {
-            implicitTerms_.push_back(eqnTerm);
+            implicitOperators_.push_back(Operator);
         }
-        for (auto& eqnTerm : eqnSys.explicitTerms_)
+        for (auto& Operator : eqnSys.explicitOperators_)
         {
-            explicitTerms_.push_back(eqnTerm);
+            explicitOperators_.push_back(Operator);
         }
     }
 
     void solve()
     {
-        if (temporalTerms_.size() == 0 && implicitTerms_.size() == 0)
+        if (temporalOperators_.size() == 0 && implicitOperators_.size() == 0)
         {
             NF_ERROR_EXIT("No temporal or implicit terms to solve.");
         }
-        if (temporalTerms_.size() > 0)
+        if (temporalOperators_.size() > 0)
         {
             // integrate equations in time
         }
@@ -85,52 +85,58 @@ class EqnSystem
 
     size_t size() const
     {
-        return temporalTerms_.size() + implicitTerms_.size() + explicitTerms_.size();
+        return temporalOperators_.size() + implicitOperators_.size() + explicitOperators_.size();
     }
 
     // getters
-    const std::vector<EqnTerm>& temporalTerms() const { return temporalTerms_; }
+    const std::vector<Operator>& temporalOperators() const { return temporalOperators_; }
 
-    const std::vector<EqnTerm>& implicitTerms() const { return implicitTerms_; }
+    const std::vector<Operator>& implicitOperators() const { return implicitOperators_; }
 
-    const std::vector<EqnTerm>& explicitTerms() const { return explicitTerms_; }
+    const std::vector<Operator>& explicitOperators() const { return explicitOperators_; }
 
-    std::vector<EqnTerm>& temporalTerms() { return temporalTerms_; }
+    std::vector<Operator>& temporalOperators() { return temporalOperators_; }
 
-    std::vector<EqnTerm>& implicitTerms() { return implicitTerms_; }
+    std::vector<Operator>& implicitOperators() { return implicitOperators_; }
 
-    std::vector<EqnTerm>& explicitTerms() { return explicitTerms_; }
+    std::vector<Operator>& explicitOperators() { return explicitOperators_; }
 
     const NeoFOAM::Executor& exec() const { return exec_; }
 
     const std::size_t nCells() const { return nCells_; }
 
+    scalar getDt() const { return dt_; }
+
     fvcc::VolumeField<NeoFOAM::scalar>* volumeField()
     {
-        if (temporalTerms_.size() == 0 && implicitTerms_.size() == 0)
+        if (temporalOperators_.size() == 0 && implicitOperators_.size() == 0)
         {
             NF_ERROR_EXIT("No temporal or implicit terms to solve.");
         }
-        if (temporalTerms_.size() > 0)
+        if (temporalOperators_.size() > 0)
         {
-            volumeField_ = temporalTerms_[0].volumeField();
+            // FIXME
+            NF_ERROR_EXIT("Not implemented.");
+            // volumeField_ = temporalOperators_[0].volumeField();
         }
         else
         {
-            volumeField_ = implicitTerms_[0].volumeField();
+            // FIXME
+            NF_ERROR_EXIT("Not implemented.");
+            // volumeField_ = implicitOperators_[0].volumeField();
         }
         return volumeField_;
     }
 
-    NeoFOAM::scalar dt = 0;
+    NeoFOAM::scalar dt_ = 0;
 
 private:
 
     const NeoFOAM::Executor exec_;
     const std::size_t nCells_;
-    std::vector<EqnTerm> temporalTerms_;
-    std::vector<EqnTerm> implicitTerms_;
-    std::vector<EqnTerm> explicitTerms_;
+    std::vector<Operator> temporalOperators_;
+    std::vector<Operator> implicitOperators_;
+    std::vector<Operator> explicitOperators_;
     fvcc::VolumeField<NeoFOAM::scalar>* volumeField_;
 };
 
@@ -140,34 +146,35 @@ EqnSystem operator+(EqnSystem lhs, const EqnSystem& rhs)
     return lhs;
 }
 
-EqnSystem operator+(EqnSystem lhs, const EqnTerm& rhs)
+EqnSystem operator+(EqnSystem lhs, const Operator& rhs)
 {
-    lhs.addTerm(rhs);
+    lhs.addOperator(rhs);
     return lhs;
 }
 
-EqnSystem operator+(const EqnTerm& lhs, const EqnTerm& rhs)
+EqnSystem operator+(const Operator& lhs, const Operator& rhs)
 {
-    EqnSystem eqnSys(lhs.exec(), lhs.nCells());
-    eqnSys.addTerm(lhs);
-    eqnSys.addTerm(rhs);
-    return eqnSys;
+    NF_ERROR_EXIT("Not implemented.");
+    //     EqnSystem eqnSys(lhs.exec(), lhs.nCells());
+    //     eqnSys.addOperator(lhs);
+    //     eqnSys.addOperator(rhs);
+    //     return eqnSys;
 }
 
 EqnSystem operator*(NeoFOAM::scalar scale, const EqnSystem& es)
 {
     EqnSystem results(es.exec(), es.nCells());
-    for (const auto& eqnTerm : es.temporalTerms())
+    for (const auto& Operator : es.temporalOperators())
     {
-        results.addTerm(scale * eqnTerm);
+        results.addOperator(scale * Operator);
     }
-    for (const auto& eqnTerm : es.implicitTerms())
+    for (const auto& Operator : es.implicitOperators())
     {
-        results.addTerm(scale * eqnTerm);
+        results.addOperator(scale * Operator);
     }
-    for (const auto& eqnTerm : es.explicitTerms())
+    for (const auto& Operator : es.explicitOperators())
     {
-        results.addTerm(scale * eqnTerm);
+        results.addOperator(scale * Operator);
     }
     return results;
 }
@@ -178,18 +185,19 @@ EqnSystem operator-(EqnSystem lhs, const EqnSystem& rhs)
     return lhs;
 }
 
-EqnSystem operator-(EqnSystem lhs, const EqnTerm& rhs)
+EqnSystem operator-(EqnSystem lhs, const Operator& rhs)
 {
-    lhs.addTerm(-1.0 * rhs);
+    lhs.addOperator(-1.0 * rhs);
     return lhs;
 }
 
-EqnSystem operator-(const EqnTerm& lhs, const EqnTerm& rhs)
+EqnSystem operator-(const Operator& lhs, const Operator& rhs)
 {
-    EqnSystem results(lhs.exec(), lhs.nCells());
-    results.addTerm(lhs);
-    results.addTerm(-1.0 * rhs);
-    return results;
+    NF_ERROR_EXIT("Not implemented.");
+    // EqnSystem results(lhs.exec(), lhs.nCells());
+    // results.addOperator(lhs);
+    // results.addOperator(-1.0 * rhs);
+    // return results;
 }
 
 
diff --git a/include/NeoFOAM/DSL/eqnTerm.hpp b/include/NeoFOAM/DSL/eqnTerm.hpp
deleted file mode 100644
index 3e2eab1e2..000000000
--- a/include/NeoFOAM/DSL/eqnTerm.hpp
+++ /dev/null
@@ -1,170 +0,0 @@
-// SPDX-License-Identifier: MIT
-// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
-#pragma once
-
-#include <iostream>
-#include <memory>
-#include <type_traits>
-#include <utility>
-
-#include "NeoFOAM/core/primitives/scalar.hpp"
-#include "NeoFOAM/fields/field.hpp"
-#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
-
-namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
-
-namespace NeoFOAM::DSL
-{
-
-template<typename T>
-concept HasTemporalTerm = requires(T t) {
-    {
-        t.temporalOperation(
-            std::declval<NeoFOAM::Field<NeoFOAM::scalar>&>(), std::declval<NeoFOAM::scalar>()
-        )
-    } -> std::same_as<void>; // Adjust return type and arguments as needed
-};
-
-template<typename T>
-concept HasExplicitTerm = requires(T t) {
-    {
-        t.explicitOperation(
-            std::declval<NeoFOAM::Field<NeoFOAM::scalar>&>(), std::declval<NeoFOAM::scalar>()
-        )
-    } -> std::same_as<void>; // Adjust return type and arguments as needed
-};
-// EqnTerm class that uses type erasure without inheritance
-class EqnTerm
-{
-public:
-
-    enum class Type
-    {
-        Temporal,
-        Implicit,
-        Explicit
-    };
-
-    template<typename T>
-    EqnTerm(T cls) : model_(std::make_unique<Model<T>>(std::move(cls)))
-    {}
-
-    EqnTerm(const EqnTerm& eqnTerm) : model_ {eqnTerm.model_->clone()} {}
-
-    EqnTerm(EqnTerm&& eqnTerm) : model_ {std::move(eqnTerm.model_)} {}
-
-    EqnTerm& operator=(const EqnTerm& eqnTerm)
-    {
-        model_ = eqnTerm.model_->clone();
-        return *this;
-    }
-
-    EqnTerm& operator=(EqnTerm&& eqnTerm)
-    {
-        model_ = std::move(eqnTerm.model_);
-        return *this;
-    }
-
-    std::string display() const { return model_->display(); }
-
-    void explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source)
-    {
-        model_->explicitOperation(source, model_->scaleCoeff);
-    }
-
-    void temporalOperation(NeoFOAM::Field<NeoFOAM::scalar>& field)
-    {
-        model_->temporalOperation(field, model_->scaleCoeff);
-    }
-
-    EqnTerm::Type getType() const { return model_->getType(); }
-
-    void setScale(NeoFOAM::scalar scale) { model_->scaleCoeff *= scale; }
-
-
-    const NeoFOAM::Executor& exec() const { return model_->exec(); }
-
-    std::size_t nCells() const { return model_->nCells(); }
-
-
-    fvcc::VolumeField<NeoFOAM::scalar>* volumeField() { return model_->volumeField(); }
-
-
-private:
-
-    // Base class to hold the type-erased value and the display function
-    struct Concept
-    {
-        virtual ~Concept() = default;
-        virtual std::string display() const = 0;
-        virtual void
-        explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source, NeoFOAM::scalar scale) = 0;
-        virtual void
-        temporalOperation(NeoFOAM::Field<NeoFOAM::scalar>& field, NeoFOAM::scalar scale) = 0;
-        NeoFOAM::scalar scaleCoeff = 1.0;
-        virtual EqnTerm::Type getType() const = 0;
-
-        virtual const NeoFOAM::Executor& exec() const = 0;
-        virtual std::size_t nCells() const = 0;
-        virtual fvcc::VolumeField<NeoFOAM::scalar>* volumeField() = 0;
-
-        // The Prototype Design Pattern
-        virtual std::unique_ptr<Concept> clone() const = 0;
-    };
-
-    // Templated derived class to implement the type-specific behavior
-    template<typename T>
-    struct Model : Concept
-    {
-        Model(T cls) : cls_(std::move(cls)) {}
-
-        std::string display() const override { return cls_.display(); }
-
-        virtual void
-        explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source, NeoFOAM::scalar scale) override
-        {
-            if constexpr (HasExplicitTerm<T>)
-            {
-                cls_.explicitOperation(source, scale);
-            }
-        }
-
-        virtual void
-        temporalOperation(NeoFOAM::Field<NeoFOAM::scalar>& field, NeoFOAM::scalar scale) override
-        {
-            if constexpr (HasTemporalTerm<T>)
-            {
-                cls_.temporalOperation(field, scale);
-            }
-        }
-
-        virtual fvcc::VolumeField<NeoFOAM::scalar>* volumeField() override
-        {
-            return cls_.volumeField();
-        }
-
-        EqnTerm::Type getType() const override { return cls_.getType(); }
-
-        const NeoFOAM::Executor& exec() const override { return cls_.exec(); }
-
-        std::size_t nCells() const override { return cls_.nCells(); }
-
-        // The Prototype Design Pattern
-        std::unique_ptr<Concept> clone() const override { return std::make_unique<Model>(*this); }
-
-        T cls_;
-    };
-
-    std::unique_ptr<Concept> model_;
-};
-
-
-// add multiply operator to EqnTerm
-EqnTerm operator*(NeoFOAM::scalar scale, const EqnTerm& lhs)
-{
-    EqnTerm result = lhs;
-    result.setScale(scale);
-    return result;
-}
-
-} // namespace NeoFOAM::DSL
diff --git a/include/NeoFOAM/DSL/operator.hpp b/include/NeoFOAM/DSL/operator.hpp
new file mode 100644
index 000000000..8ac92d543
--- /dev/null
+++ b/include/NeoFOAM/DSL/operator.hpp
@@ -0,0 +1,246 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: 2023-2024 NeoFOAM authors
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <concepts>
+
+#include "NeoFOAM/core/primitives/scalar.hpp"
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/core/parallelAlgorithms.hpp"
+#include "NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp"
+#include "NeoFOAM/core/input.hpp"
+#include "NeoFOAM/DSL/coeff.hpp"
+
+
+namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
+
+namespace NeoFOAM::DSL
+{
+
+template<typename T>
+concept HasTemporalOperator = requires(T t) {
+    {
+        t.temporalOperation(std::declval<Field<scalar>&>(), std::declval<scalar>())
+    } -> std::same_as<void>; // Adjust return type and arguments as needed
+};
+
+template<typename T>
+concept HasExplicitOperator = requires(T t) {
+    {
+        t.explicitOperation(std::declval<Field<scalar>&>())
+    } -> std::same_as<void>; // Adjust return type and arguments as needed
+};
+
+
+/* @class OperatorMixin
+ * @brief A mixin class to represent simplify implementations of concrete operators
+ * in NeoFOAMs DSL
+ *
+ *
+ * @ingroup DSL
+ */
+class OperatorMixin
+{
+
+public:
+
+    OperatorMixin(const Executor exec) : exec_(exec), coeffs_() {};
+
+    virtual ~OperatorMixin() = default;
+
+    const Executor& exec() const { return exec_; }
+
+    Coeff& getCoefficient() { return coeffs_; }
+
+    const Coeff& getCoefficient() const { return coeffs_; }
+
+    /* @brief Given an input this function reads required coeffs */
+    void build(const Input& input) {}
+
+protected:
+
+    const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)
+
+    Coeff coeffs_;
+};
+
+
+/* @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<typename T>
+    Operator(T cls) : model_(std::make_unique<OperatorModel<T>>(std::move(cls)))
+    {}
+
+    Operator(const Operator& eqnOperator) : model_ {eqnOperator.model_->clone()} {}
+
+    Operator(Operator&& eqnOperator) : model_ {std::move(eqnOperator.model_)} {}
+
+    void explicitOperation(Field<scalar>& source) { model_->explicitOperation(source); }
+
+    void temporalOperation(Field<scalar>& 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<scalar>& source) = 0;
+
+        virtual void temporalOperation(Field<scalar>& 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<OperatorConcept> clone() const = 0;
+    };
+
+    // Templated derived class to implement the type-specific behavior
+    template<typename ConcreteOperatorType>
+    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<scalar>& source) override
+        {
+            if constexpr (HasExplicitOperator<ConcreteOperatorType>)
+            {
+                concreteOp_.explicitOperation(source);
+            }
+        }
+
+        virtual void temporalOperation(Field<scalar>& field) override
+        {
+            if constexpr (HasTemporalOperator<ConcreteOperatorType>)
+            {
+                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<OperatorConcept> clone() const override
+        {
+            return std::make_unique<OperatorModel>(*this);
+        }
+
+        ConcreteOperatorType concreteOp_;
+    };
+
+    std::unique_ptr<OperatorConcept> model_;
+};
+
+
+auto operator*(scalar scalarCoeff, const Operator& rhs)
+{
+    Operator result = rhs;
+    result.getCoefficient() *= scalarCoeff;
+    return result;
+}
+
+auto operator*(const Field<scalar>& 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<typename CoeffFunction>
+    requires std::invocable<CoeffFunction&, size_t>
+Operator operator*(CoeffFunction coeffFunc, const Operator& lhs)
+{
+    // TODO implement
+    NF_ERROR_EXIT("Not implemented");
+    Operator result = lhs;
+    // if (!result.getCoefficient().useSpan)
+    // {
+    //     result.setField(std::make_shared<Field<scalar>>(result.exec(), result.nCells(), 1.0));
+    // }
+    // map(result.exec(), result.getCoefficient().values, scaleFunc);
+    return result;
+}
+
+} // namespace NeoFOAM::DSL
diff --git a/include/NeoFOAM/core/dictionary.hpp b/include/NeoFOAM/core/dictionary.hpp
index cd06b2dc9..94cec83fe 100644
--- a/include/NeoFOAM/core/dictionary.hpp
+++ b/include/NeoFOAM/core/dictionary.hpp
@@ -162,6 +162,12 @@ class Dictionary
      */
     const std::unordered_map<std::string, std::any>& 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<std::string, std::any> data_;
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/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<ValueType>& rhs)
+    DomainField(const DomainField<ValueType>& rhs)
         : exec_(rhs.exec_), internalField_(rhs.internalField_), boundaryFields_(rhs.boundaryFields_)
     {}
 
@@ -66,7 +66,7 @@ class DomainField
     {}
 
 
-    DomainField<ValueType>& operator=(const ::NeoFOAM::DomainField<ValueType>& rhs)
+    DomainField<ValueType>& operator=(const DomainField<ValueType>& rhs)
     {
         internalField_ = rhs.internalField_;
         boundaryFields_ = rhs.boundaryFields_;
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
index 53b68a47d..c02ee804a 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/geometricField.hpp
@@ -36,11 +36,34 @@ class GeometricFieldMixin
      * @param field The domain field object.
      */
     GeometricFieldMixin(
-        const Executor& exec, const UnstructuredMesh& mesh, const DomainField<ValueType>& field
+        const Executor& exec,
+        const UnstructuredMesh& mesh,
+        const DomainField<ValueType>& 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 field The domain field object.
+     */
+    GeometricFieldMixin(
+        const Executor& exec,
+        const UnstructuredMesh& mesh,
+        const Field<ValueType>& internalField,
+        const BoundaryFields<ValueType>& 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.
      *
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<ValueType>
     SurfaceField(
         const Executor& exec,
         const UnstructuredMesh& mesh,
-        const Field<ValueType>& internalField,
+        const DomainField<ValueType>& domainField,
         const std::vector<SurfaceBoundary<ValueType>>& boundaryConditions
     )
-        : GeometricFieldMixin<ValueType>(exec, mesh, {exec, mesh, internalField}),
+        : GeometricFieldMixin<ValueType>(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<ValueType>& internalField,
+        const BoundaryFields<ValueType>& boundaryFields,
+        const std::vector<SurfaceBoundary<ValueType>>& boundaryConditions
+    )
+        : GeometricFieldMixin<ValueType>(exec, mesh, {exec, mesh, internalField, boundaryFields}),
+          boundaryConditions_(boundaryConditions)
+    {}
 
     SurfaceField(const SurfaceField& other)
         : GeometricFieldMixin<ValueType>(other), boundaryConditions_(other.boundaryConditions_)
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp b/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp
index ad590f425..f9ba40463 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/fields/volumeField.hpp
@@ -46,6 +46,22 @@ class VolumeField : public GeometricFieldMixin<ValueType>
           boundaryConditions_(boundaryConditions)
     {}
 
+    /* @brief Constructor for a VolumeField with a given internal field
+     *
+     * @param mesh The underlying mesh
+     * @param internalField the underlying internal field
+     * @param boundaryConditions a vector of boundary conditions
+     */
+    VolumeField(
+        const Executor& exec,
+        const UnstructuredMesh& mesh,
+        const DomainField<ValueType>& domainField,
+        const std::vector<VolumeBoundary<ValueType>>& boundaryConditions
+    )
+        : GeometricFieldMixin<ValueType>(exec, mesh, domainField),
+          boundaryConditions_(boundaryConditions)
+    {}
+
     /* @brief Constructor for a VolumeField with a given internal field
      *
      * @param mesh The underlying mesh
@@ -56,9 +72,10 @@ class VolumeField : public GeometricFieldMixin<ValueType>
         const Executor& exec,
         const UnstructuredMesh& mesh,
         const Field<ValueType>& internalField,
+        const BoundaryFields<ValueType>& boundaryFields,
         const std::vector<VolumeBoundary<ValueType>>& boundaryConditions
     )
-        : GeometricFieldMixin<ValueType>(exec, mesh, {exec, mesh, internalField}),
+        : GeometricFieldMixin<ValueType>(exec, mesh, internalField, boundaryFields),
           boundaryConditions_(boundaryConditions)
     {}
 
diff --git a/include/NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp b/include/NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp
index 63fcd57aa..9709f0746 100644
--- a/include/NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp
+++ b/include/NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp
@@ -8,7 +8,7 @@
 #include "NeoFOAM/mesh/unstructured.hpp"
 #include "NeoFOAM/finiteVolume/cellCentred.hpp"
 
-#include "NeoFOAM/DSL/eqnTerm.hpp"
+#include "NeoFOAM/DSL/operator.hpp"
 #include "NeoFOAM/DSL/eqnSystem.hpp"
 
 #include <functional>
diff --git a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
index d2b4ca994..a43d204d0 100644
--- a/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
+++ b/src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
@@ -33,7 +33,6 @@ void computeDiv(
     size_t nInternalFaces = mesh.nInternalFaces();
     const auto surfV = mesh.cellVolumes().span();
 
-
     // check if the executor is GPU
     if (std::holds_alternative<SerialExecutor>(exec))
     {
@@ -51,7 +50,6 @@ void computeDiv(
             surfDivPhi[own] += valueOwn;
         }
 
-
         for (size_t celli = 0; celli < mesh.nCells(); celli++)
         {
             surfDivPhi[celli] *= 1 / surfV[celli];
diff --git a/src/finiteVolume/cellCentred/timeIntegration/forwardEuler.cpp b/src/finiteVolume/cellCentred/timeIntegration/forwardEuler.cpp
index b67d42cae..08ce009f2 100644
--- a/src/finiteVolume/cellCentred/timeIntegration/forwardEuler.cpp
+++ b/src/finiteVolume/cellCentred/timeIntegration/forwardEuler.cpp
@@ -20,7 +20,7 @@ ForwardEuler::ForwardEuler(const dsl::EqnSystem& eqnSystem, const Dictionary& di
 void ForwardEuler::solve()
 {
     std::cout << "Solving using Forward Euler" << std::endl;
-    scalar dt = eqnSystem_.dt;
+    scalar dt = eqnSystem_.getDt();
     fvcc::VolumeField<scalar>* refField = eqnSystem_.volumeField();
     // Field<scalar> Phi(eqnSystem_.exec(), eqnSystem_.nCells());
     // NeoFOAM::fill(Phi, 0.0);
diff --git a/test/core/CMakeLists.txt b/test/core/CMakeLists.txt
index 4bdb9b185..7c57a79ab 100644
--- a/test/core/CMakeLists.txt
+++ b/test/core/CMakeLists.txt
@@ -6,7 +6,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 <catch2/catch_test_macros.hpp>
 #include <catch2/generators/catch_generators_all.hpp>
 
-#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
index 06ec5792d..5616fabe1 100644
--- a/test/dsl/CMakeLists.txt
+++ b/test/dsl/CMakeLists.txt
@@ -2,4 +2,6 @@
 # SPDX-FileCopyrightText: 2024 NeoFOAM authors
 
 neofoam_unit_test(coeff)
+neofoam_unit_test(operator)
+
 # FIXME neofoam_unit_test(dsl)
diff --git a/test/dsl/coeff.cpp b/test/dsl/coeff.cpp
index 1d3cbd0c1..e7c33958e 100644
--- a/test/dsl/coeff.cpp
+++ b/test/dsl/coeff.cpp
@@ -12,6 +12,7 @@
 
 using Field = NeoFOAM::Field<NeoFOAM::scalar>;
 using Coeff = NeoFOAM::DSL::Coeff;
+
 namespace detail = NeoFOAM::DSL::detail;
 
 
@@ -100,7 +101,7 @@ TEST_CASE("Coeff")
         {
             NeoFOAM::Field<NeoFOAM::scalar> fieldA(exec, size, 0.0);
             NeoFOAM::Field<NeoFOAM::scalar> fieldB(exec, size, 1.0);
-            Coeff coeff {-5.0,fieldB};
+            Coeff coeff {-5.0, fieldB};
             {
                 NeoFOAM::parallelFor(
                     fieldA, KOKKOS_LAMBDA(const size_t i) { return coeff[i] + 2.0; }
@@ -112,9 +113,5 @@ TEST_CASE("Coeff")
             REQUIRE(hostFieldA[1] == -3.0);
             REQUIRE(hostFieldA[2] == -3.0);
         }
-
     }
-
-        
-    
 }
diff --git a/test/dsl/dsl.cpp b/test/dsl/dsl.cpp
index 0aa1d02e7..48fafa000 100644
--- a/test/dsl/dsl.cpp
+++ b/test/dsl/dsl.cpp
@@ -64,7 +64,6 @@ class Divergence
 
     dsl::EqnTerm::Type getType() const { return termType_; }
 
-
     const NeoFOAM::Executor& exec() const { return exec_; }
 
     const std::size_t nCells() const { return nCells_; }
diff --git a/test/dsl/operator.cpp b/test/dsl/operator.cpp
new file mode 100644
index 000000000..ec8870b0a
--- /dev/null
+++ b/test/dsl/operator.cpp
@@ -0,0 +1,116 @@
+// 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 <catch2/catch_session.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators_all.hpp>
+#include <catch2/benchmark/catch_benchmark.hpp>
+
+#include "NeoFOAM/fields/field.hpp"
+#include "NeoFOAM/DSL/coeff.hpp"
+#include "NeoFOAM/DSL/operator.hpp"
+
+using Field = NeoFOAM::Field<NeoFOAM::scalar>;
+using Coeff = NeoFOAM::DSL::Coeff;
+using Operator = NeoFOAM::DSL::Operator;
+using OperatorMixin = NeoFOAM::DSL::OperatorMixin;
+using Executor = NeoFOAM::Executor;
+using VolumeField = fvcc::VolumeField<NeoFOAM::scalar>;
+using BoundaryFields = NeoFOAM::BoundaryFields<NeoFOAM::scalar>;
+namespace fvcc = NeoFOAM::finiteVolume::cellCentred;
+
+/* A dummy implementation of a Operator
+ * following the Operator interface */
+class Dummy : public OperatorMixin
+{
+
+public:
+
+    Dummy(const Executor& exec, VolumeField& field) : OperatorMixin(exec), field_(field) {}
+
+    void explicitOperation(Field& source)
+    {
+        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"; }
+
+    const VolumeField& volumeField() const { return field_; }
+
+    VolumeField& volumeField() { return field_; }
+
+    Operator::Type getType() const { return Operator::Type::Explicit; }
+
+private:
+
+    VolumeField& field_;
+};
+
+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<fvcc::VolumeBoundary<NeoFOAM::scalar>> bcs {};
+        auto vf = VolumeField(exec, mesh, fA, bf, bcs);
+        auto b = Dummy(exec, vf);
+
+        REQUIRE(b.getName() == "Dummy");
+        REQUIRE(b.getType() == Operator::Type::Explicit);
+    }
+
+    SECTION("Supports Coefficients" + execName)
+    {
+        std::vector<fvcc::VolumeBoundary<NeoFOAM::scalar>> 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(exec, vf);
+        auto d = fB * Dummy(exec, vf);
+        auto e = Coeff(-3, fB) * Dummy(exec, vf);
+
+        auto coeffc = c.getCoefficient();
+        auto coeffd = d.getCoefficient();
+        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/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp b/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp
index e276213f1..b1f3b64ce 100644
--- a/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp
+++ b/test/finiteVolume/cellCentred/timeIntegration/timeIntegration.cpp
@@ -10,15 +10,12 @@
 #include "NeoFOAM/finiteVolume/cellCentred/timeIntegration/timeIntegration.hpp"
 #include "NeoFOAM/core/dictionary.hpp"
 #include "NeoFOAM/core/parallelAlgorithms.hpp"
-
-
-#include "NeoFOAM/DSL/eqnTerm.hpp"
+#include "NeoFOAM/DSL/operator.hpp"
 #include "NeoFOAM/DSL/eqnSystem.hpp"
 
-
 namespace dsl = NeoFOAM::DSL;
 
-class Divergence
+class Divergence : public NeoFOAM::DSL::OperatorMixin
 {
 
 public:
@@ -35,7 +32,7 @@ class Divergence
         );
     }
 
-    dsl::EqnTerm::Type getType() const { return termType_; }
+    dsl::Operator::Type getType() const { return termType_; }
 
     fvcc::VolumeField<NeoFOAM::scalar>* volumeField() const { return nullptr; }
 
@@ -43,17 +40,19 @@ class Divergence
 
     std::size_t nCells() const { return nCells_; }
 
-    dsl::EqnTerm::Type termType_;
+    dsl::Operator::Type termType_;
+
     NeoFOAM::Executor exec_;
+
     std::size_t nCells_;
 };
 
-class TimeTerm
+class TimeOperator
 {
 
 public:
 
-    std::string display() const { return "TimeTerm"; }
+    std::string display() const { return "TimeOperator"; }
 
     void explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source, NeoFOAM::scalar scale)
     {
@@ -65,7 +64,7 @@ class TimeTerm
         );
     }
 
-    dsl::EqnTerm::Type getType() const { return termType_; }
+    dsl::Operator::Type getType() const { return termType_; }
 
     fvcc::VolumeField<NeoFOAM::scalar>* volumeField() const { return nullptr; }
 
@@ -73,7 +72,7 @@ class TimeTerm
 
     std::size_t nCells() const { return nCells_; }
 
-    dsl::EqnTerm::Type termType_;
+    dsl::Operator::Type termType_;
     const NeoFOAM::Executor exec_;
     std::size_t nCells_;
 };
@@ -87,11 +86,11 @@ TEST_CASE("TimeIntegration")
     NeoFOAM::Dictionary dict;
     dict.insert("type", std::string("forwardEuler"));
 
-    dsl::EqnTerm divTerm = Divergence(dsl::EqnTerm::Type::Explicit, exec, 1);
+    // dsl::Operator divOperator = Divergence(dsl::Operator::Type::Explicit, exec, 1);
 
-    dsl::EqnTerm ddtTerm = TimeTerm(dsl::EqnTerm::Type::Temporal, exec, 1);
+    // dsl::Operator ddtOperator = TimeOperator(dsl::Operator::Type::Temporal, exec, 1);
 
-    dsl::EqnSystem eqnSys = ddtTerm + divTerm;
+    // dsl::EqnSystem eqnSys = ddtOperator + divOperator;
 
-    fvcc::TimeIntegration timeIntergrator(eqnSys, dict);
+    // fvcc::TimeIntegration timeIntergrator(eqnSys, dict);
 }