DSL implementation stack

.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 }}"

CHANGELOG.md

@@ -1,4 +1,5 @@
 # Version 0.1.0 (unreleased)
+- Implement a DSL interface []()
 - 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)
 - added launch json to debug unit test in vscode [#135](https://github.com/exasim-project/NeoFOAM/pull/135)

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.
 

doc/conf.py

@@ -23,6 +23,7 @@
 
 extensions = [
     'sphinx.ext.autodoc',
+    'sphinxcontrib.mermaid',
     'sphinx.ext.intersphinx',
     'sphinx.ext.autosectionlabel',
     'sphinx.ext.todo',

doc/dsl/equation.rst

@@ -0,0 +1,50 @@
+Equation
+---------
+
+
+The `Equation` class in NeoFOAM holds, manages, builds and solves the expressed/programmed equation.
+Its core responsibility lie 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 `Equation` stores 3 vectors:
+
+.. mermaid::
+
+    classDiagram
+        class Operator {
+            +explicitOperation(...)
+            +implicitOperation(...)
+        }
+        class DivOperator {
+            +explicitOperation(...)
+            +implicitOperation(...)
+        }
+        class TemporalOperator {
+            +explicitOperation(...)
+            +implicitOperation(...)
+        }
+        class Others["..."] {
+            +explicitOperation(...)
+            +implicitOperation(...)
+        }
+        class Equation {
+            +temporalTerms_: vector~Operator~
+            +implicitTerms_: vector~Operator~
+            +explicitTerms_: vector~Operator~
+        }
+        Operator <|-- DivOperator
+        Operator <|-- TemporalOperator
+        Operator <|-- Others
+        Equation <|-- Operator
+
+Thus, an `Equation` 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`.

doc/dsl/index.rst

@@ -0,0 +1,71 @@
+.. _fvcc:
+
+Domain Specific Language (DSL)
+==============================
+
+The concept of a Domain Specific Language (DSL) allows to dramatically 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 and 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
+    - the sparse matrix is always a 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 format 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<NeoFOAM::scalar> UEqn
+    (
+        fvcc::impOp::ddt(U)
+        + fvcc::impOp::div(phi, U)
+        - fvcc::impOp::laplacian(nu, U)
+    )
+
+    solve(UEqn == -fvcc::expOp::grad(p));
+
+
+In contrast to OpenFOAM, here the majority of the work is done in 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

doc/dsl/operator.rst

@@ -0,0 +1,118 @@
+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 a an additional coefficient, for example a scalar value or a further field.
+The `Operator` implementation uses Type Erasure (more details `[1] <https://medium.com/@gealleh/type-erasure-idiom-in-c-0d1cb4f61cf0>`_ `[2] <https://www.youtube.com/watch?v=4eeESJQk-mw>`_ `[3] <https://www.youtube.com/watch?v=qn6OqefuH08>`_) 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<NeoFOAM::scalar> divTerm =
+            Divergence(NeoFOAM::dsl::Operator<NeoFOAM::scalar>::Type::Explicit, exec, ...);
+
+        NeoFOAM::dsl::Operator<NeoFOAM::scalar> ddtTerm =
+            TimeTerm(NeoFOAM::dsl::Operator<NeoFOAM::scalar>::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<NeoFOAM::scalar> scalingField(exec, nCells, 2.0);
+        auto sF = scalingField.span();
+
+        dsl::Operator<NeoFOAM::scalar> customTerm =
+            CustomTerm(dsl::Operator<NeoFOAM::scalar>::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
+    - display: display the term
+    - getType: get the type of the term
+    - exec: get the executor
+    - nCells: get the number of cells
+    - volumeField: get the volume field
+
+An example is given below:
+
+.. code-block:: cpp
+
+    class CustomOperator : public dsl::OperatorMixin<NeoFOAM::scalar>
+    {
+
+    public:
+
+        // constructors ..
+        NeoFOAM::scalar read(const NeoFOAM::Input& input)
+        {
+            // ..
+        }
+
+        void build(const NeoFOAM::Input& input)
+        {
+            value = read(input);
+            termEvaluated = true;
+        }
+
+        std::string display() const { return "Laplacian"; }
+
+        void explicitOperation(NeoFOAM::Field<NeoFOAM::scalar>& source)
+        {
+            NeoFOAM::scalar setValue = value;
+            // scaleField is defined in OperatorMixin
+            // and accounts for the scaling of the terms
+            // and considers scaling by fields and scalars
+            auto scale = scaleField();
+            auto sourceField = source.span();
+            NeoFOAM::parallelFor(
+                source.exec(),
+                {0, source.size()},
+                KOKKOS_LAMBDA(const size_t i) { sourceField[i] += scale[i] * setValue; }
+            );
+        }
+
+        // other helper functions
+        dsl::Operator<NeoFOAM::scalar>::Type getType() const { return termType_; }
+
+        const NeoFOAM::Executor& exec() const { return exec_; }
+
+        std::size_t nCells() const { return nCells_; }
+
+        fvcc::VolumeField<NeoFOAM::scalar>* volumeField() { return nullptr; }
+
+        dsl::Operator<NeoFOAM::scalar>::Type termType_;
+
+
+        const NeoFOAM::Executor exec_;
+        std::size_t nCells_;
+        NeoFOAM::scalar value = 1.0;
+    };
+
+The required scaling of the term is handle by the `scaleField` function, provided by `OperatorMixin`. The `scaleField` function returns the 'ScalingField' class that is used to scale by fields and scalars.
+
+.. code-block:: cpp
+
+    template <typename ValueType>
+    class ScalingField
+    {
+
+        // the span is only used if it is defined
+        KOKKOS_INLINE_FUNCTION
+        ValueType operator[](const size_t i) const { return useSpan ? values[i] * value : value; }
+
+    }

doc/finiteVolume/cellCentred/index.rst

@@ -8,8 +8,6 @@ cellCenteredFiniteVolume
     :maxdepth: 2
     :glob:
 
-    DSL.rst
-    fields.rst
     boundaryConditions.rst
     operators.rst
     stencil.rst

doc/index.rst

@@ -26,6 +26,7 @@ Table of Contents
    contributing
    basics/index
    finiteVolume/cellCentred/index
+   dsl/index
    api/index
 
 Compatibility with OpenFOAM

doc/requirements.txt

@@ -2,3 +2,4 @@ sphinx
 sphinx-sitemap
 furo
 breathe
+sphinxcontrib-mermaid

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 <cxxabi.h>
-#include <string>
-#include <any>
-#include <iostream>
+
+#include <any>      // for bad_any_cast
+#include <iostream> // for operator<<, basic_ostream, cerr, endl, ostream
+#include <string>   // for operator<<, string
+#include <typeinfo> // for type_info
+
 namespace NeoFOAM
 {
 

include/NeoFOAM/core/dictionary.hpp

@@ -5,11 +5,9 @@
 #include <unordered_map>
 #include <any>
 #include <string>
-#include <iostream>
 #include <vector>
 
 #include "NeoFOAM/core/demangle.hpp"
-#include "NeoFOAM/core/error.hpp"
 
 namespace NeoFOAM
 {
@@ -162,6 +160,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_;

include/NeoFOAM/core/executor/CPUExecutor.hpp

@@ -2,8 +2,7 @@
 // SPDX-FileCopyrightText: 2023 NeoFOAM authors
 #pragma once
 
-#include <Kokkos_Core.hpp>
-#include <iostream>
+#include <Kokkos_Core.hpp> // IWYU pragma: keep
 
 namespace NeoFOAM
 {

include/NeoFOAM/core/executor/GPUExecutor.hpp

@@ -3,7 +3,6 @@
 #pragma once
 
 #include <Kokkos_Core.hpp>
-#include <iostream>
 
 namespace NeoFOAM
 {

include/NeoFOAM/core/executor/serialExecutor.hpp

@@ -3,7 +3,6 @@
 #pragma once
 
 #include <Kokkos_Core.hpp>
-#include <iostream>
 
 namespace NeoFOAM
 {