Merge pull request #147 from exasim-project/doc/dsl

Documentation of NeoFOAM DSL, see #147

.github/workflows/build_doc.yaml

@@ -5,8 +5,6 @@ env:
 
 on:
  pull_request:
-    branches:
-      - main
     types: [closed, synchronize, opened]
  push:
     tags:
@@ -21,6 +19,7 @@ jobs:
         && pip3 install furo
         && pip3 install breathe
         && pip3 install sphinx-sitemap
+        && pip3 install sphinxcontrib-mermaid
     - name: Checkout repo
       uses: actions/checkout@v4
     - name: Build docs
@@ -46,6 +45,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 }}"

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/conf.py

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

doc/finiteVolume/cellCentred/index.rst

@@ -9,7 +9,6 @@ cellCenteredFiniteVolume
     :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