Merge branch 'trm-model-auto-eval' into 'master'

[TRM] More automation in constitutive setting evaluation.

See merge request ogs/ogs!5001

ProcessLib/Graph/Apply.h

@@ -16,34 +16,41 @@ namespace ProcessLib::Graph
 namespace detail
 {
 template <typename Function>
-struct GetFunctionArgumentTypesPlain  // plain, i.e., without cvref
+struct GetFunctionArgumentTypes
     /** \cond */
-    : GetFunctionArgumentTypesPlain<decltype(&Function::operator())>
+    : GetFunctionArgumentTypes<decltype(&Function::operator())>
 /** \endcond */
 {
 };
 
 // member function
 template <typename Result, typename Class, typename... Args>
-struct GetFunctionArgumentTypesPlain<Result (Class::*)(Args...)>
+struct GetFunctionArgumentTypes<Result (Class::*)(Args...)>
 {
-    using type = boost::mp11::mp_list<std::remove_cvref_t<Args>...>;
+    using type = boost::mp11::mp_list<Args...>;
 };
 
 // const member function
 template <typename Result, typename Class, typename... Args>
-struct GetFunctionArgumentTypesPlain<Result (Class::*)(Args...) const>
+struct GetFunctionArgumentTypes<Result (Class::*)(Args...) const>
 {
-    using type = boost::mp11::mp_list<std::remove_cvref_t<Args>...>;
+    using type = boost::mp11::mp_list<Args...>;
 };
 
 // standalone function
 template <typename Result, typename... Args>
-struct GetFunctionArgumentTypesPlain<Result (*)(Args...)>
+struct GetFunctionArgumentTypes<Result (*)(Args...)>
 {
-    using type = boost::mp11::mp_list<std::remove_cvref_t<Args>...>;
+    using type = boost::mp11::mp_list<Args...>;
 };
 
+// plain, i.e., without cvref
+template <typename Function>
+using GetFunctionArgumentTypesPlain =
+    std::type_identity<boost::mp11::mp_transform<
+        std::remove_cvref_t,
+        typename GetFunctionArgumentTypes<Function>::type>>;
+
 template <typename Function>
 struct GetFunctionReturnType
     /** \cond */
@@ -273,4 +280,19 @@ auto eval(Function& f, Tuples&... ts) ->
 
     return detail::applyImpl(&Function::eval, f, ts...);
 }
+
+/**
+ * Invokes the eval() method of the passed objects \c fs with arguments taken
+ * from the passed tuples.
+ *
+ * \c fs must be a tuple of objects having an eval() method. The method
+ * invocation proceeds in the order of the objects in the tuple.
+ *
+ * \see eval()
+ */
+template <typename Functions, typename... Tuples>
+void evalAllInOrder(Functions& fs, Tuples&... ts)
+{
+    boost::mp11::tuple_for_each(fs, [&ts...](auto& f) { eval(f, ts...); });
+}
 }  // namespace ProcessLib::Graph

ProcessLib/Graph/CheckEvalOrderRT.h

@@ -0,0 +1,136 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <boost/mp11.hpp>
+#include <typeindex>
+#include <unordered_set>
+
+#include "Apply.h"
+#include "BaseLib/Logging.h"
+
+namespace ProcessLib::Graph
+{
+namespace detail
+{
+template <typename T>
+struct IsInputArgument
+    : boost::mp11::mp_bool<std::is_lvalue_reference_v<T> &&
+                           std::is_const_v<std::remove_reference_t<T>>>
+{
+    static_assert(std::is_lvalue_reference_v<T>,
+                  "The current implementation only deals with l-value "
+                  "references as function arguments. If you want to extend it, "
+                  "test thoroughly in order to not introduce bugs.");
+};
+
+template <typename T>
+struct IsOutputArgument : boost::mp11::mp_bool<!IsInputArgument<T>::value>
+{
+    static_assert(std::is_lvalue_reference_v<T>,
+                  "The current implementation only deals with l-value "
+                  "references as function arguments. If you want to extend it, "
+                  "test thoroughly in order to not introduce bugs.");
+};
+
+template <typename Model>
+bool isEvalOrderCorrectRT(std::unordered_set<std::type_index>& computed_data)
+{
+    using namespace boost::mp11;
+
+    using ModelArgs =
+        typename GetFunctionArgumentTypes<decltype(&Model::eval)>::type;
+    using ModelInputs = mp_filter<IsInputArgument, ModelArgs>;
+    using ModelOutputs = mp_filter<IsOutputArgument, ModelArgs>;
+
+    using ModelInputsWrapped = mp_transform<mp_identity, ModelInputs>;
+
+    // Check that all inputs have already been computed before.
+    bool all_inputs_computed = true;
+    mp_for_each<ModelInputsWrapped>(
+        [&computed_data,
+         &all_inputs_computed]<typename Input>(mp_identity<Input>)
+        {
+            if (!computed_data.contains(std::type_index{typeid(Input)}))
+            {
+                ERR("Input {} of model {} has not been computed/set before the "
+                    "model evaluation.",
+                    typeid(Input).name(), typeid(Model).name());
+                all_inputs_computed = false;
+            }
+        });
+    if (!all_inputs_computed)
+    {
+        return false;
+    }
+
+    using ModelOutputsWrapped = mp_transform<mp_identity, ModelOutputs>;
+
+    // All outputs are "computed data", now.
+    bool no_output_precomputed = true;
+    mp_for_each<ModelOutputsWrapped>(
+        [&computed_data,
+         &no_output_precomputed]<typename Output>(mp_identity<Output>)
+        {
+            auto const [it, emplaced] = computed_data.emplace(typeid(Output));
+
+            if (!emplaced)
+            {
+                ERR("Output {} of model {} is computed more than once.",
+                    typeid(Output).name(),
+                    typeid(Model).name());
+                no_output_precomputed = false;
+            }
+        });
+
+    return no_output_precomputed;
+}
+
+template <typename... Models>
+bool isEvalOrderCorrectRT(boost::mp11::mp_list<Models...>,
+                          std::unordered_set<std::type_index>&& computed_data)
+{
+    return (isEvalOrderCorrectRT<Models>(computed_data) && ...);
+}
+}  // namespace detail
+
+/// Checks at runtime if the given \c Models are evaluated in the right order if
+/// evaluated in the order in which they appear in the list of \c Models.
+///
+/// I.e., all input data of a model must have been computed before that model
+/// will be evaluated and no two models must compute the same data.
+///
+/// The passed \c Inputs are data that already have been computed before the
+/// first model is evaluated.
+template <typename Models, typename Inputs>
+bool isEvalOrderCorrectRT()  // RT for runtime
+{
+    using namespace boost::mp11;
+
+    static_assert(mp_is_list<Models>::value);
+    static_assert(mp_is_list<Inputs>::value);
+
+    // Wrap inputs. The elements of InputsWrapped are default constructible.
+    using InputsWrapped = mp_transform<mp_identity, Inputs>;
+
+    // "Holds" all data that has been computed successively by the invoked
+    // models.
+    std::unordered_set<std::type_index> computed_data;
+
+    // All inputs are considered "computed data".
+    mp_for_each<InputsWrapped>(
+        [&computed_data]<typename Input>(mp_identity<Input>)
+        { computed_data.emplace(typeid(Input)); });
+
+    return detail::isEvalOrderCorrectRT(mp_rename<Models, mp_list>{},
+                                        std::move(computed_data));
+}
+}  // namespace ProcessLib::Graph

ProcessLib/Graph/ConstructModels.h

@@ -0,0 +1,75 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <tuple>
+#include <utility>
+
+#include "Apply.h"
+
+namespace ProcessLib::Graph
+{
+namespace detail
+{
+template <typename T>
+concept HasCreate = requires
+{
+    T::create;
+};
+
+template <typename Model, typename TupleOfArgs>
+Model constructModel(TupleOfArgs& args)
+{
+    if constexpr (HasCreate<Model>)
+    {
+        return applyImpl(&Model::create, args);
+    }
+    else if constexpr (std::is_default_constructible_v<Model>)
+    {
+        return Model{};
+    }
+    else
+    {
+        static_assert(std::is_default_constructible_v<
+                          Model> /* This is definitely false, here. */,
+                      "The model to be constructed has neither a static "
+                      "create() function nor is it default constructible.");
+    }
+}
+template <template <typename...> typename Tuple,
+          typename... Models,
+          typename TupleOfArgs>
+Tuple<Models...> constructModels(std::type_identity<Tuple<Models...>>,
+                                 TupleOfArgs&& args)
+{
+    return Tuple{constructModel<Models>(args)...};
+}
+}  // namespace detail
+
+/**
+ * Constructs a tuple of models.
+ *
+ * Each model in the tuple is either
+ *
+ * 1. constructed by calling a static create() method of the model's class,
+ * 2. or by default construction if the former does not exist.
+ *
+ * In case (1.) the arguments are passed via their data types from the passed
+ * \c args to the create() method.
+ */
+template <typename TupleOfModels, typename... Args>
+TupleOfModels constructModels(Args&&... args)
+{
+    return detail::constructModels(
+        std::type_identity<TupleOfModels>{},
+        std::forward_as_tuple(std::forward<Args>(args)...));
+}
+}  // namespace ProcessLib::Graph

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h

@@ -13,6 +13,7 @@
 #include "LiquidDensity.h"
 #include "LiquidViscosity.h"
 #include "PermeabilityData.h"
+#include "SpecificBodyForceData.h"
 #include "ThermoOsmosis.h"
 
 namespace ProcessLib::ThermoRichardsMechanics
@@ -47,6 +48,12 @@ struct DarcyLawModel
               ThermoOsmosisData<DisplacementDim> const& th_osmosis_data,
               DarcyLawData<DisplacementDim>& out) const;
 
+    static DarcyLawModel create(
+        SpecificBodyForceData<DisplacementDim> const& specific_body_force_data)
+    {
+        return DarcyLawModel{specific_body_force_data.specific_body_force};
+    }
+
 private:
     /// Gravity vector (specific body force).
     Eigen::Vector<double, DisplacementDim> const b_;

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h

@@ -16,6 +16,7 @@
 #include "LiquidViscosity.h"
 #include "PermeabilityData.h"
 #include "Saturation.h"
+#include "SpecificBodyForceData.h"
 #include "TRMStorage.h"
 #include "TRMVaporDiffusion.h"
 
@@ -59,6 +60,12 @@ struct EqPModel
               TRMStorageData const& storage_data,
               EqPData<DisplacementDim>& out) const;
 
+    static EqPModel create(
+        SpecificBodyForceData<DisplacementDim> const& specific_body_force_data)
+    {
+        return EqPModel{specific_body_force_data.specific_body_force};
+    }
+
 private:
     /// Gravity vector (specific body force).
     Eigen::Vector<double, DisplacementDim> const b_;

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Gravity.h

@@ -13,6 +13,7 @@
 #include "LiquidDensity.h"
 #include "Saturation.h"
 #include "SolidDensity.h"
+#include "SpecificBodyForceData.h"
 
 namespace ProcessLib::ThermoRichardsMechanics
 {
@@ -39,6 +40,12 @@ struct GravityModel
               SaturationDataDeriv const& dS_L_data,
               GravityData<DisplacementDim>& out) const;
 
+    static GravityModel create(
+        SpecificBodyForceData<DisplacementDim> const& specific_body_force_data)
+    {
+        return GravityModel{specific_body_force_data.specific_body_force};
+    }
+
 private:
     Eigen::Vector<double, DisplacementDim> const specific_body_force_;
 };

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SpecificBodyForceData.h

@@ -0,0 +1,21 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <Eigen/Core>
+
+namespace ProcessLib::ThermoRichardsMechanics
+{
+template <int DisplacementDim>
+struct SpecificBodyForceData
+{
+    Eigen::Vector<double, DisplacementDim> specific_body_force;
+};
+}  // namespace ProcessLib::ThermoRichardsMechanics

ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveData.h

@@ -16,6 +16,7 @@
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqT.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EquivalentPlasticStrainData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Gravity.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidDensity.h"