[MPL] Orthotropic permeability power law.

Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/c_PermeabilityOrthotropicPowerLaw.md

+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw

Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_exponents.md

+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw::_lambda

Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_intrinsic_permeabilities.md

+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw::_k

MaterialLib/MPL/CreateProperty.cpp

@@ -65,6 +65,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createIdealGasLaw(config);
     }
 
+    if (boost::iequals(property_type, "PermeabilityOrthotropicPowerLaw"))
+    {
+        return createPermeabilityOrthotropicPowerLaw(config,
+                                                     local_coordinate_system);
+    }
+
     if (boost::iequals(property_type, "PorosityFromMassBalance"))
     {
         return createPorosityFromMassBalance(config, parameters);

MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "ParameterLib/CoordinateSystem.h"
+#include "PermeabilityOrthotropicPowerLaw.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPermeabilityOrthotropicPowerLaw(
+    BaseLib::ConfigTree const& config,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "PermeabilityOrthotropicPowerLaw");
+    DBUG("Create PermeabilityOrthotropicPowerLaw solid phase property");
+
+    auto const intrinsic_permeabilities =
+        //! \ogs_file_param{properties__property__PermeabilityOrthotropicPowerLaw__intrinsic_permeabilities}
+        config.getConfigParameter<std::vector<double>>(
+            "intrinsic_permeabilities");
+
+    if (!((intrinsic_permeabilities.size() == 3) ||
+        (intrinsic_permeabilities.size() == 2)))
+    {
+        OGS_FATAL(
+            "The number of intrinsic permeabilities must be two or three, but "
+            "%d were given.",
+            intrinsic_permeabilities.size());
+    }
+
+    auto const exponents =
+        //! \ogs_file_param{properties__property__PermeabilityOrthotropicPowerLaw__exponents}
+        config.getConfigParameter<std::vector<double>>("exponents");
+
+    if (exponents.size() != 3 && exponents.size() != 2)
+    {
+        OGS_FATAL(
+            "The number of exponents must be two or three, but %d were given.",
+            exponents.size());
+    }
+
+    if (intrinsic_permeabilities.size() != exponents.size())
+    {
+        OGS_FATAL(
+            "The number of intrinsic permeabilities and exponents must be "
+            "equal, but they are %d and %d, respectively.",
+            intrinsic_permeabilities.size(), exponents.size());
+    }
+
+    if (exponents.size() == 2)
+    {
+        return std::make_unique<PermeabilityOrthotropicPowerLaw<2>>(
+            std::array<double, 2>{intrinsic_permeabilities[0],
+                                  intrinsic_permeabilities[1]},
+            std::array<double, 2>{exponents[0], exponents[1]},
+            local_coordinate_system);
+    }
+    if (exponents.size() == 3)
+    {
+        return std::make_unique<PermeabilityOrthotropicPowerLaw<3>>(
+            std::array<double, 3>{intrinsic_permeabilities[0],
+                                  intrinsic_permeabilities[1],
+                                  intrinsic_permeabilities[2]},
+            std::array<double, 3>{exponents[0], exponents[1], exponents[2]},
+            local_coordinate_system);
+    }
+    OGS_FATAL(
+        "Could not create PermeabilityOrthotropicPowerLaw material model.");
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 ParameterLib
+{
+struct CoordinateSystem;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPermeabilityOrthotropicPowerLaw(
+    BaseLib::ConfigTree const& config,
+    ParameterLib::CoordinateSystem const* const local_coordinate_system);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

@@ -17,6 +17,7 @@
 #include "CreateIdealGasLaw.h"
 #include "CreateLinearProperty.h"
 #include "CreateParameterProperty.h"
+#include "CreatePermeabilityOrthotropicPowerLaw.h"
 #include "CreatePorosityFromMassBalance.h"
 #include "CreateRelPermBrooksCorey.h"
 #include "CreateRelPermLiakopoulos.h"

MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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/PermeabilityOrthotropicPowerLaw.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MathLib/KelvinVector.h"
+#include "ParameterLib/CoordinateSystem.h"
+
+namespace MaterialPropertyLib
+{
+template <int DisplacementDim>
+PermeabilityOrthotropicPowerLaw<DisplacementDim>::
+    PermeabilityOrthotropicPowerLaw(
+        std::array<double, DisplacementDim> intrinsic_permeabilities,
+        std::array<double, DisplacementDim>
+            exponents,
+        ParameterLib::CoordinateSystem const* const local_coordinate_system)
+    : _k(std::move(intrinsic_permeabilities)),
+      _lambda(std::move(exponents)),
+      _local_coordinate_system(local_coordinate_system)
+{
+}
+
+template <int DisplacementDim>
+void PermeabilityOrthotropicPowerLaw<DisplacementDim>::setScale(
+    std::variant<Medium*, Phase*, Component*> scale_pointer)
+{
+    if (std::holds_alternative<Phase*>(scale_pointer))
+    {
+        _phase = std::get<Phase*>(scale_pointer);
+        if (_phase->name != "Solid")
+        {
+            OGS_FATAL(
+                "The property 'PermeabilityOrthotropicPowerLaw' must be "
+                "given in the 'Solid' phase, not in '%s' phase.",
+                _phase->name.c_str());
+        }
+    }
+    else
+    {
+        OGS_FATAL(
+            "The property 'PermeabilityOrthotropicPowerLaw' is "
+            "implemented on the 'phase' scales only.");
+    }
+}
+template <int DisplacementDim>
+PropertyDataType PermeabilityOrthotropicPowerLaw<DisplacementDim>::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& pos, double const /*t*/,
+    double const /*dt*/) const
+{
+    auto const phi =
+        std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
+    // TODO (naumov) The phi0 must be evaluated once upon
+    // creation/initialization and be stored in a local state.
+    // For now assume porosity's initial value does not change with time.
+    auto const phi_0 = _phase->property(PropertyType::porosity)
+                           .template initialValue<double>(
+                               pos, std::numeric_limits<double>::quiet_NaN());
+
+    Eigen::Matrix<double, DisplacementDim, DisplacementDim> k =
+        Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+
+    Eigen::Matrix<double, DisplacementDim, DisplacementDim> const e =
+        _local_coordinate_system == nullptr
+            ? Eigen::Matrix<double, DisplacementDim,
+                            DisplacementDim>::Identity()
+            : _local_coordinate_system->transformation<DisplacementDim>(pos);
+
+    // k = \sum_i k_i (\phi / \phi_0)^{\lambda_i} e_i \otimes e_i
+    // e_i \otimes e_i = square matrix e_i,0^2 e_i,0*e_i,1 etc.
+    for (int i = 0; i < DisplacementDim; ++i)
+    {
+        Eigen::Matrix<double, DisplacementDim, DisplacementDim> const
+            ei_otimes_ei = e.col(i) * e.col(i).transpose();
+
+        k += _k[i] * std::pow(phi / phi_0, _lambda[i]) * ei_otimes_ei;
+    }
+    return k;
+}
+template <int DisplacementDim>
+PropertyDataType PermeabilityOrthotropicPowerLaw<DisplacementDim>::dValue(
+    VariableArray const& /*variable_array*/, Variable const primary_variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    (void)primary_variable;
+    assert((primary_variable == Variable::strain) &&
+           "PermeabilityOrthotropicPowerLaw::dValue is implemented for "
+           " derivatives with respect to strain only.");
+
+    return 0;
+}
+
+template class PermeabilityOrthotropicPowerLaw<2>;
+template class PermeabilityOrthotropicPowerLaw<3>;
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.h

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "MaterialLib/MPL/Property.h"
+
+
+namespace ParameterLib
+{
+struct CoordinateSystem;
+}
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+
+/// Orthotropic permeability model based on power law dependency on porosity.
+/// \details This property must be a solid phase property, it
+/// computes the permeability depending on the porosity. A local coordinate
+/// system can be given for orthotropy.
+template <int DisplacementDim>
+class PermeabilityOrthotropicPowerLaw final : public Property
+{
+private:
+    Phase* _phase = nullptr;
+    /// Intrinsic permeabilities, one for each spatial dimension.
+    std::array<double, DisplacementDim> const _k;
+    /// Exponents, one for each spatial dimension.
+    std::array<double, DisplacementDim> const _lambda;
+    ParameterLib::CoordinateSystem const* const _local_coordinate_system;
+
+public:
+    PermeabilityOrthotropicPowerLaw(
+        std::array<double, DisplacementDim> intrinsic_permeabilities,
+        std::array<double, DisplacementDim>
+            exponents,
+        ParameterLib::CoordinateSystem const* const local_coordinate_system);
+
+    void setScale(
+        std::variant<Medium*, Phase*, Component*> scale_pointer) override;
+
+    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 variable,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t, double const dt) const override;
+};
+
+}  // namespace MaterialPropertyLib