temporalOperator.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2023 - 2026 NeoN authors
//
// SPDX-License-Identifier: MIT
 
#pragma once
 
#include <memory>
#include <concepts>
 
#include "NeoN/core/error.hpp"
#include "NeoN/core/primitives/scalar.hpp"
#include "NeoN/core/vector/vector.hpp"
#include "NeoN/core/input.hpp"
#include "NeoN/linearAlgebra/linearSystem.hpp"
#include "NeoN/linearAlgebra/faceToMatrixAddress.hpp"
#include "NeoN/dsl/coeff.hpp"
#include "NeoN/dsl/operator.hpp"
#include "NeoN/finiteVolume/cellCentred/operators/ddtOperator.hpp"
 
namespace NeoN::dsl
{
 
template<typename T>
concept HasTemporalExplicitOperator = requires(T t) {
    {
        t.explicitOperation(
            std::declval<Vector<typename T::VectorValueType>&>(),
            std::declval<NeoN::scalar>(),
            std::declval<NeoN::scalar>()
        )
    } -> std::same_as<void>; // Adjust return type and arguments as needed
};
 
template<typename T>
concept HasTemporalImplicitOperator = requires(T t) {
    {
        t.implicitOperation(
            std::declval<la::LinearSystem<typename T::VectorValueType>&>(),
            std::declval<NeoN::scalar>(),
            std::declval<NeoN::scalar>()
        )
    } -> std::same_as<void>; // Adjust return type and arguments as needed
};
 
/* @brief Concept satisfied when T can assemble its temporal contribution into a
 *        LinearSystem whose matrix coefficients are scalar while the RHS holds T's
 *        field value type (segregated vector-solve form).
 */
template<typename T>
concept HasTemporalImplicitOperatorScalarMtx = requires(T t) {
    {
        t.implicitOperation(
            std::declval<la::LinearSystem<scalar, typename T::VectorValueType>&>(),
            std::declval<NeoN::scalar>(),
            std::declval<NeoN::scalar>()
        )
    } -> std::same_as<void>;
};
 
template<typename T>
concept HasTemporalOperator = HasTemporalExplicitOperator<T> || HasTemporalImplicitOperator<T>
                           || HasTemporalImplicitOperatorScalarMtx<T>;
 
/* @class TemporalOperator
 * @brief A class to represent a TemporalOperator in NeoNs 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 TemporalOperator e.g ddt, d2dt2, etc can be stored in a vector of
 * TemporalOperator
 *
 * @ingroup dsl
 */
template<typename ValueType>
class TemporalOperator
{
public:
 
    using VectorValueType = ValueType;
 
    template<HasTemporalOperator T>
    TemporalOperator(T cls) : model_(std::make_unique<TemporalOperatorModel<T>>(std::move(cls)))
    {}
 
    TemporalOperator(const TemporalOperator& eqnOperator) : model_ {eqnOperator.model_->clone()} {}
 
    TemporalOperator(TemporalOperator&& eqnOperator) : model_ {std::move(eqnOperator.model_)} {}
 
    TemporalOperator& operator=(const TemporalOperator& eqnOperator)
    {
        model_ = eqnOperator.model_->clone();
        return *this;
    }
 
    void explicitOperation(Vector<ValueType>& source, scalar t, scalar dt) const
    {
        model_->explicitOperation(source, t, dt);
    }
 
    void implicitOperation(la::LinearSystem<ValueType>& ls, scalar t, scalar dt) const
    {
        model_->implicitOperation(ls, t, dt);
    }
 
    /* @brief Implicit temporal assembly into a scalar-matrix / ValueType-rhs linear system
     *        (segregated vector-solve form). Disabled when ValueType == scalar to avoid
     *        colliding with the same-type overload above.
     */
    template<typename U = ValueType>
        requires(!std::is_same_v<U, scalar>)
    void implicitOperation(la::LinearSystem<scalar, ValueType>& ls, scalar t, scalar dt) const
    {
        model_->implicitOperationScalarMtx(ls, t, dt);
    }
 
    /* returns the fundamental type of an operator, ie explicit, implicit */
    Operator::Type getType() const { return model_->getType(); }
 
    std::string getName() const { return model_->getName(); }
 
    Coeff& getCoefficient() { return model_->getCoefficient(); }
 
    Coeff getCoefficient() const { return model_->getCoefficient(); }
 
    /* @brief Given an input this function reads required properties */
    void read(const Input& input) { model_->read(input); }
 
    /* @brief Get the executor */
    const Executor& exec() const { return model_->exec(); }
 
    /* @brief Get the ddtScheme */
    NeoN::finiteVolume::cellCentred::DdtScheme ddtScheme() const { return model_->ddtScheme(); }
 
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 TemporalOperatorConcept
    {
        virtual ~TemporalOperatorConcept() = default;
 
        virtual void explicitOperation(Vector<ValueType>& source, scalar t, scalar dt) = 0;
 
        virtual void implicitOperation(la::LinearSystem<ValueType>& ls, scalar t, scalar dt) = 0;
 
        /* @brief Temporal assembly into LinearSystem<scalar, ValueType> for the
         *        scalar-matrix / ValueType-rhs (segregated vector-solve) form.
         *        Concrete operators that don't support this form leave it as a no-op.
         */
        virtual void implicitOperationScalarMtx(
            la::LinearSystem<scalar, ValueType>& ls, scalar t, scalar dt
        ) = 0;
 
        /* @brief Given an input this function reads required properties */
        virtual void read(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;
 
        /* @brief Get the ddtScheme */
        virtual NeoN::finiteVolume::cellCentred::DdtScheme ddtScheme() const
        {
            return NeoN::finiteVolume::cellCentred::DdtScheme::None;
        }
 
        // The Prototype Design Pattern
        virtual std::unique_ptr<TemporalOperatorConcept> clone() const = 0;
    };
 
    // Templated derived class to implement the type-specific behavior
    template<typename ConcreteTemporalOperatorType>
    struct TemporalOperatorModel : TemporalOperatorConcept
    {
        /* @brief build with concrete TemporalOperator */
        TemporalOperatorModel(ConcreteTemporalOperatorType concreteOp)
            : concreteOp_(std::move(concreteOp))
        {}
 
        /* returns the name of the operator */
        std::string getName() const override { return concreteOp_.getName(); }
 
        virtual void explicitOperation(Vector<ValueType>& source, scalar t, scalar dt) override
        {
            if constexpr (HasTemporalExplicitOperator<ConcreteTemporalOperatorType>)
            {
                concreteOp_.explicitOperation(source, t, dt);
            }
        }
 
        virtual void
        implicitOperation(la::LinearSystem<ValueType>& ls, scalar t, scalar dt) override
        {
            if constexpr (HasTemporalImplicitOperator<ConcreteTemporalOperatorType>)
            {
                concreteOp_.implicitOperation(ls, t, dt);
            }
        }
 
        virtual void implicitOperationScalarMtx(
            la::LinearSystem<scalar, ValueType>& ls, scalar t, scalar dt
        ) override
        {
            if constexpr (HasTemporalImplicitOperatorScalarMtx<ConcreteTemporalOperatorType>)
            {
                concreteOp_.implicitOperation(ls, t, dt);
            }
            else
            {
                // Reached only for an implicit temporal operator that lacks the scalar-matrix
                // (segregated vector-solve) overload. Silently skipping it would drop its
                // contribution (e.g. an implicit ddt term) and yield a wrong system, so fail
                // fast instead.
                NF_ERROR_EXIT(
                    "Temporal operator '" << getName()
                                          << "' does not support scalar-matrix (segregated) "
                                             "assembly."
                );
            }
        }
 
        /* @brief Given an input this function reads required coeffs */
        virtual void read(const Input& input) override { concreteOp_.read(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(); }
 
        /* @brief return the ddtScheme read by the ddtOperator */
        NeoN::finiteVolume::cellCentred::DdtScheme ddtScheme() const override
        {
            if constexpr (requires { concreteOp_.scheme(); })
            {
                return concreteOp_.scheme();
            }
            else
            {
                return NeoN::finiteVolume::cellCentred::DdtScheme::None;
            }
        }
        // The Prototype Design Pattern
        std::unique_ptr<TemporalOperatorConcept> clone() const override
        {
            return std::make_unique<TemporalOperatorModel>(*this);
        }
 
        ConcreteTemporalOperatorType concreteOp_;
    };
 
    std::unique_ptr<TemporalOperatorConcept> model_;
};
 
 
template<typename ValueType>
TemporalOperator<ValueType> operator*(scalar scalarCoeff, TemporalOperator<ValueType> rhs)
{
    TemporalOperator<ValueType> result = rhs;
    result.getCoefficient() *= scalarCoeff;
    return result;
}
 
template<typename ValueType>
TemporalOperator<ValueType>
operator*(const Vector<scalar>& coeffVector, TemporalOperator<ValueType> rhs)
{
    TemporalOperator<ValueType> result = rhs;
    result.getCoefficient() *= Coeff {coeffVector};
    return result;
}
 
template<typename ValueType>
TemporalOperator<ValueType> operator*(const Coeff& coeff, TemporalOperator<ValueType> rhs)
{
    TemporalOperator<ValueType> result = rhs;
    result.getCoefficient() *= coeff;
    return result;
}
 
 
} // namespace NeoN::dsl