From 43f2aace61566eeff81041b14e407035a9c27b4c Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Wed, 26 May 2021 18:00:02 +0200
Subject: [PATCH] [MPL] Function-type Property using exprtk.
A generalized property type for (almost) arbitrary functions. Based on
the exprtk library also used in the Function-type parameters.
In this initial implementation all *scalar* variables from the variable
array passed for the value and derivative evaluation can be used. Return
values are either scalar, 2- or 3-vectors, and 2x2- or 3x3 matrices.
The evaluation is somewhat slower than the dedicated implementation; in
the case of component transport process a bilinear property was replaced
by exactly the same function-type property and the total run-time for
the ctest has increased by around 6%.
---
.../property/Function/c_Function.md | 1 +
.../property/Function/dvalue/i_dvalue.md | 1 +
.../property/Function/dvalue/t_expression.md | 1 +
.../Function/dvalue/t_variable_name.md | 2 +
.../property/Function/value/i_value.md | 1 +
.../property/Function/value/t_expression.md | 1 +
MaterialLib/CMakeLists.txt | 2 +-
MaterialLib/MPL/CreateProperty.cpp | 4 +
MaterialLib/MPL/Properties/CreateFunction.cpp | 60 +++++
MaterialLib/MPL/Properties/CreateFunction.h | 27 +++
MaterialLib/MPL/Properties/CreateProperties.h | 1 +
MaterialLib/MPL/Properties/Function.cpp | 216 ++++++++++++++++++
MaterialLib/MPL/Properties/Function.h | 59 +++++
MaterialLib/MPL/Properties/Properties.h | 1 +
14 files changed, 376 insertions(+), 1 deletion(-)
create mode 100644 Documentation/ProjectFile/properties/property/Function/c_Function.md
create mode 100644 Documentation/ProjectFile/properties/property/Function/dvalue/i_dvalue.md
create mode 100644 Documentation/ProjectFile/properties/property/Function/dvalue/t_expression.md
create mode 100644 Documentation/ProjectFile/properties/property/Function/dvalue/t_variable_name.md
create mode 100644 Documentation/ProjectFile/properties/property/Function/value/i_value.md
create mode 100644 Documentation/ProjectFile/properties/property/Function/value/t_expression.md
create mode 100644 MaterialLib/MPL/Properties/CreateFunction.cpp
create mode 100644 MaterialLib/MPL/Properties/CreateFunction.h
create mode 100644 MaterialLib/MPL/Properties/Function.cpp
create mode 100644 MaterialLib/MPL/Properties/Function.h
diff --git a/Documentation/ProjectFile/properties/property/Function/c_Function.md b/Documentation/ProjectFile/properties/property/Function/c_Function.md
new file mode 100644
index 00000000000..5a7b4408a4c
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/c_Function.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::Function
diff --git a/Documentation/ProjectFile/properties/property/Function/dvalue/i_dvalue.md b/Documentation/ProjectFile/properties/property/Function/dvalue/i_dvalue.md
new file mode 100644
index 00000000000..25a872c3e27
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/dvalue/i_dvalue.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::Function::dvalue_expressions_
diff --git a/Documentation/ProjectFile/properties/property/Function/dvalue/t_expression.md b/Documentation/ProjectFile/properties/property/Function/dvalue/t_expression.md
new file mode 100644
index 00000000000..25a872c3e27
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/dvalue/t_expression.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::Function::dvalue_expressions_
diff --git a/Documentation/ProjectFile/properties/property/Function/dvalue/t_variable_name.md b/Documentation/ProjectFile/properties/property/Function/dvalue/t_variable_name.md
new file mode 100644
index 00000000000..4a95bbbd569
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/dvalue/t_variable_name.md
@@ -0,0 +1,2 @@
+Variable name for the derivative of the function with respect to the named
+variable.
diff --git a/Documentation/ProjectFile/properties/property/Function/value/i_value.md b/Documentation/ProjectFile/properties/property/Function/value/i_value.md
new file mode 100644
index 00000000000..c3f64a50506
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/value/i_value.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::Function::value_expressions_
diff --git a/Documentation/ProjectFile/properties/property/Function/value/t_expression.md b/Documentation/ProjectFile/properties/property/Function/value/t_expression.md
new file mode 100644
index 00000000000..c3f64a50506
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/Function/value/t_expression.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::Function::value_expressions_
diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index f44b84266f2..6767d5f4dd5 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -40,5 +40,5 @@ ogs_add_library(MaterialLib ${SOURCES})
target_link_libraries(
MaterialLib PUBLIC MaterialLib_SolidModels MaterialLib_FractureModels
- PRIVATE MathLib MeshLib ParameterLib spdlog::spdlog
+ PRIVATE MathLib MeshLib ParameterLib exprtk spdlog::spdlog
)
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 4ea28761bef..d3dba9812fd 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -56,6 +56,10 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
{
return createExponential(config);
}
+ if (property_type == "Function")
+ {
+ return createFunction(config);
+ }
if (property_type == "Parameter")
{
diff --git a/MaterialLib/MPL/Properties/CreateFunction.cpp b/MaterialLib/MPL/Properties/CreateFunction.cpp
new file mode 100644
index 00000000000..e6048cb4f4e
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateFunction.cpp
@@ -0,0 +1,60 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "Function.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Function> createFunction(BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "Function");
+
+ // Second access for storage.
+ //! \ogs_file_param{properties__property__name}
+ auto property_name = config.peekConfigParameter<std::string>("name");
+
+ DBUG("Create Function property {:s}.", property_name);
+
+ std::vector<std::string> value_expressions;
+ //! \ogs_file_param{properties__property__Function__value}
+ auto const& value_config = config.getConfigSubtree("value");
+
+ //! \ogs_file_param{properties__property__Function__value__expression}
+ for (auto const& p : value_config.getConfigSubtreeList("expression"))
+ {
+ value_expressions.emplace_back(p.getValue<std::string>());
+ }
+
+ // For each derivative a name of the variable and the list of expressions.
+ std::vector<std::pair<std::string, std::vector<std::string>>>
+ dvalue_expressions;
+ //! \ogs_file_param{properties__property__Function__dvalue}
+ for (auto const& dvalue_config : config.getConfigSubtreeList("dvalue"))
+ {
+ auto variable_name =
+ //! \ogs_file_param{properties__property__Function__dvalue__variable_name}
+ dvalue_config.getConfigParameter<std::string>("variable_name");
+
+ std::vector<std::string> expressions;
+ //! \ogs_file_param{properties__property__Function__dvalue__expression}
+ for (auto const& p : dvalue_config.getConfigSubtreeList("expression"))
+ {
+ expressions.emplace_back(p.getValue<std::string>());
+ }
+ dvalue_expressions.emplace_back(std::move(variable_name),
+ std::move(expressions));
+ }
+
+ return std::make_unique<MaterialPropertyLib::Function>(
+ std::move(property_name), std::move(value_expressions),
+ std::move(dvalue_expressions));
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateFunction.h b/MaterialLib/MPL/Properties/CreateFunction.h
new file mode 100644
index 00000000000..0f86f3d02b6
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateFunction.h
@@ -0,0 +1,27 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class Function;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Function> createFunction(BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 7b45b0a11a4..de51caa2889 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -27,6 +27,7 @@
#include "CreateEffectiveThermalConductivityPorosityMixing.h"
#include "CreateEmbeddedFracturePermeability.h"
#include "CreateExponential.h"
+#include "CreateFunction.h"
#include "CreateGasPressureDependentPermeability.h"
#include "CreateIdealGasLaw.h"
#include "CreateKozenyCarmanModel.h"
diff --git a/MaterialLib/MPL/Properties/Function.cpp b/MaterialLib/MPL/Properties/Function.cpp
new file mode 100644
index 00000000000..bbdb668e633
--- /dev/null
+++ b/MaterialLib/MPL/Properties/Function.cpp
@@ -0,0 +1,216 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "MaterialLib/MPL/Properties/Function.h"
+
+#include <numeric>
+
+#include "BaseLib/Algorithm.h"
+
+namespace MaterialPropertyLib
+{
+// Passing symbol table by reference as required by register_symbol_table()
+// call.
+template <typename T>
+static std::vector<exprtk::expression<T>> compileExpressions(
+ exprtk::symbol_table<T>& symbol_table,
+ std::vector<std::string> const& string_expressions)
+{
+ exprtk::parser<T> parser;
+
+ std::vector<exprtk::expression<T>> expressions(string_expressions.size());
+ expressions.resize(string_expressions.size());
+ for (unsigned i = 0; i < string_expressions.size(); ++i)
+ {
+ expressions[i].register_symbol_table(symbol_table);
+ if (!parser.compile(string_expressions[i], expressions[i]))
+ {
+ OGS_FATAL("Error: {:s}\tExpression: {:s}\n",
+ parser.error(),
+ string_expressions[i]);
+ }
+ }
+ return expressions;
+}
+
+static void updateVariableValues(
+ std::vector<std::pair<int, double*>> const& symbol_values,
+ VariableArray const& variable_array)
+{
+ for (auto& index_value_ptr_pair : symbol_values)
+ {
+ auto const index = index_value_ptr_pair.first;
+
+ double* value_ptr = index_value_ptr_pair.second;
+ std::visit(
+ [&value_ptr, &index](auto&& v) {
+ using T = std::decay_t<decltype(v)>;
+ if constexpr (std::is_same_v<T, std::monostate>)
+ {
+ OGS_FATAL(
+ "Function property: variable {:s} value needed for "
+ "evaluation of the expression was not set by the "
+ "caller.",
+ variable_enum_to_string[index]);
+ }
+ else if constexpr (std::is_same_v<T, double>)
+ {
+ *value_ptr = v;
+ }
+ else
+ {
+ OGS_FATAL(
+ "Function property: not implemented handling for a "
+ "type {:s} of variable {:s}.",
+ typeid(T).name(), variable_enum_to_string[index]);
+ }
+ },
+ variable_array[index]);
+ }
+}
+
+static PropertyDataType evaluateExpressions(
+ std::vector<std::pair<int, double*>> const& symbol_values,
+ VariableArray const& variable_array,
+ std::vector<exprtk::expression<double>> const& expressions)
+{
+ updateVariableValues(symbol_values, variable_array);
+
+ std::vector<double> result(expressions.size());
+ std::transform(begin(expressions), end(expressions), begin(result),
+ [](auto const& e) { return e.value(); });
+
+ switch (result.size())
+ {
+ case 1:
+ {
+ return result[0];
+ }
+ case 2:
+ {
+ return Eigen::Vector2d{result[0], result[1]};
+ }
+ case 3:
+ {
+ return Eigen::Vector3d{result[0], result[1], result[2]};
+ }
+ case 4:
+ {
+ Eigen::Matrix<double, 2, 2> m;
+ m = Eigen::Map<Eigen::Matrix<double, 2, 2> const>(result.data(), 2,
+ 2);
+ return m;
+ }
+ case 9:
+ {
+ Eigen::Matrix<double, 3, 3> m;
+ m = Eigen::Map<Eigen::Matrix<double, 3, 3> const>(result.data(), 3,
+ 3);
+ return m;
+ }
+ }
+ OGS_FATAL("Cannot convert a vector of size {} to a PropertyDataType",
+ result.size());
+}
+
+static std::vector<std::string> collectVariables(
+ std::vector<std::string> const& value_string_expressions,
+ std::vector<std::pair<std::string, std::vector<std::string>>> const&
+ dvalue_string_expressions)
+{
+ std::vector<std::string> variables;
+
+ auto collect_variables = [&](auto string_expressions) {
+ for (auto const& string_expression : string_expressions)
+ {
+ if (!exprtk::collect_variables(string_expression, variables))
+ {
+ OGS_FATAL("Collecting variables didn't work.");
+ }
+ }
+ };
+
+ collect_variables(value_string_expressions);
+ for (auto const& var_name_string_expression : dvalue_string_expressions)
+ {
+ collect_variables(var_name_string_expression.second);
+ }
+
+ BaseLib::makeVectorUnique(variables);
+ return variables;
+}
+
+Function::Function(std::string name,
+ std::vector<std::string>
+ value_string_expressions,
+ std::vector<std::pair<std::string, std::vector<std::string>>>
+ dvalue_string_expressions)
+{
+ name_ = std::move(name);
+
+ // Create symbol table for used variables.
+ exprtk::symbol_table<double> symbol_table;
+
+ for (auto const& v :
+ collectVariables(value_string_expressions, dvalue_string_expressions))
+ {
+ symbol_table.create_variable(v);
+ // Store variables index in the variable array and the pointer to the
+ // value in the symbol table for fast access later.
+ int const variable_array_index =
+ static_cast<int>(convertStringToVariable(v));
+ symbol_values_.emplace_back(variable_array_index,
+ &symbol_table.get_variable(v)->ref());
+ }
+
+ // value expressions.
+ value_expressions_ =
+ compileExpressions(symbol_table, value_string_expressions);
+
+ // dValue expressions.
+ for (auto const& [variable_name, string_expressions] :
+ dvalue_string_expressions)
+ {
+ dvalue_expressions_.emplace_back(
+ convertStringToVariable(variable_name),
+ compileExpressions(symbol_table, string_expressions));
+ }
+}
+
+PropertyDataType Function::value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/, double const /*dt*/) const
+{
+ return evaluateExpressions(symbol_values_, variable_array,
+ value_expressions_);
+}
+
+PropertyDataType Function::dValue(VariableArray const& variable_array,
+ Variable const primary_variable,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/, double const /*dt*/) const
+{
+ auto const it = std::find_if(begin(dvalue_expressions_),
+ end(dvalue_expressions_),
+ [&primary_variable](auto const& v) {
+ return v.first == primary_variable;
+ });
+
+ if (it == end(dvalue_expressions_))
+ {
+ OGS_FATAL(
+ "Requested derivative with respect to the variable {:s} not "
+ "provided for Function-type property {:s}.",
+ variable_enum_to_string[static_cast<int>(primary_variable)], name_);
+ }
+
+ return evaluateExpressions(symbol_values_, variable_array, it->second);
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/Function.h b/MaterialLib/MPL/Properties/Function.h
new file mode 100644
index 00000000000..27ad4c1fa2e
--- /dev/null
+++ b/MaterialLib/MPL/Properties/Function.h
@@ -0,0 +1,59 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <exprtk.hpp>
+#include <utility>
+#include <vector>
+
+#include "MaterialLib/MPL/Property.h"
+#include "MaterialLib/MPL/VariableType.h"
+
+namespace MaterialPropertyLib
+{
+/// A function property defined by mathematical expression. For the evaluation
+/// of the expressions the exprtk library is used. In the expressions all
+/// variables defined in MaterialPropertyLib::Variable enum can be used.
+///
+/// \warning The evaluation calls are not to be used in parallel (openMP),
+/// because the values' updates are using the same space.
+class Function final : public Property
+{
+public:
+ Function(std::string name,
+ std::vector<std::string>
+ value_string_expressions,
+ std::vector<std::pair<std::string, std::vector<std::string>>>
+ dvalue_string_expressions);
+
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t,
+ double const dt) const override;
+
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const primary_variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t,
+ double const dt) const override;
+
+private:
+ using Expression = exprtk::expression<double>;
+
+ /// Mapping from variable array index to symbol table values.
+ std::vector<std::pair<int, double*>> symbol_values_;
+ /// Value expressions.
+ /// Multiple expressions are representing vector-valued functions.
+ std::vector<Expression> value_expressions_;
+ /// Derivative expressions with respect to the variable.
+ /// Multiple expressions are representing vector-valued functions.
+ std::vector<std::pair<Variable, std::vector<Expression>>>
+ dvalue_expressions_;
+};
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index f22a6061d3c..290a0e43669 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -28,6 +28,7 @@
#include "Enthalpy/LinearWaterVapourLatentHeat.h"
#include "Enthalpy/WaterVapourLatentHeatWithCriticalTemperature.h"
#include "Exponential.h"
+#include "Function.h"
#include "GasPressureDependentPermeability.h"
#include "IdealGasLaw.h"
#include "Linear.h"
--
GitLab