diff --git a/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/c_PermeabilityOrthotropicPowerLaw.md b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/c_PermeabilityOrthotropicPowerLaw.md
new file mode 100644
index 0000000000000000000000000000000000000000..e45d28d54534f8e9b60b38e6cb6cbdcc42f52701
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/c_PermeabilityOrthotropicPowerLaw.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_exponents.md b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_exponents.md
new file mode 100644
index 0000000000000000000000000000000000000000..af8a83d9d2fecfdbf9652df9f836f9e876a03d9c
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_exponents.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw::_lambda
diff --git a/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_intrinsic_permeabilities.md b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_intrinsic_permeabilities.md
new file mode 100644
index 0000000000000000000000000000000000000000..ca5090d0f2ad6129b74de4f5b373c7e025bb80bd
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/PermeabilityOrthotropicPowerLaw/t_intrinsic_permeabilities.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PermeabilityOrthotropicPowerLaw::_k
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 8455f1822bd5adaeef7200e83444f1047208b818..f92f5c8f4ec3d3834c0c284e74964afa6c36a346 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/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);
diff --git a/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.cpp b/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..092d9b4f8b111ca26ee094e8758db6e090ad7736
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.cpp
@@ -0,0 +1,77 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.h b/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ab258340b0b143d7713ec3d7a239916ca3d15ea
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreatePermeabilityOrthotropicPowerLaw.h
@@ -0,0 +1,28 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 83c339efbecc7721878d19d024a91675a992a01b..61fabacdd4ed09219211cfc74d52cacabbda6698 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/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"
diff --git a/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp b/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da261bf716808df7405854385d4d906f38b6adb7
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp
@@ -0,0 +1,107 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.h b/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.h
new file mode 100644
index 0000000000000000000000000000000000000000..175a4c873e51ea274e869da2e0a058f018f72b47
--- /dev/null
+++ b/MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.h
@@ -0,0 +1,64 @@
+/**
+ * \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