[MPL] Porosity from mass balance model.

Documentation/ProjectFile/properties/property/PorosityFromMassBalance/c_PorosityFromMassBalance.md

+\copydoc MaterialPropertyLib::PorosityFromMassBalance

Documentation/ProjectFile/properties/property/PorosityFromMassBalance/t_initial_porosity.md

+\copydoc MaterialPropertyLib::PorosityFromMassBalance::_phi0

MaterialLib/MPL/CreateProperty.cpp

@@ -59,6 +59,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createIdealGasLaw(config);
     }
 
+    if (boost::iequals(property_type, "PorosityFromMassBalance"))
+    {
+        return createPorosityFromMassBalance(config, parameters);
+    }
+
     if (boost::iequals(property_type, "SaturationBrooksCorey"))
     {
         return createSaturationBrooksCorey(config);

MaterialLib/MPL/Properties/CreatePorosityFromMassBalance.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 "CreatePorosityFromMassBalance.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "ParameterLib/Utils.h"
+#include "PorosityFromMassBalance.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<PorosityFromMassBalance> createPorosityFromMassBalance(
+    BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "PorosityFromMassBalance");
+    DBUG("Create PorosityFromMassBalance medium property");
+
+    std::string const& parameter_name =
+        //! \ogs_file_param{properties__property__PorosityFromMassBalance__initial_porosity}
+        config.getConfigParameter<std::string>("initial_porosity");
+    auto const& initial_porosity = ParameterLib::findParameter<double>(
+        parameter_name, parameters, 0, nullptr);
+
+    return std::make_unique<PorosityFromMassBalance>(initial_porosity);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreatePorosityFromMassBalance.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>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class PorosityFromMassBalance;
+}
+
+namespace ParameterLib
+{
+struct ParameterBase;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<PorosityFromMassBalance> createPorosityFromMassBalance(
+    BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+        parameters);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

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

MaterialLib/MPL/Properties/PorosityFromMassBalance.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 "PorosityFromMassBalance.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+PropertyDataType PorosityFromMassBalance::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& pos,
+    double const t, double const dt) const
+{
+    double const K_SR = _phase->property(PropertyType::bulk_modulus)
+                            .template value<double>(variable_array, pos, t, dt);
+    auto const alpha_b =
+        _phase->property(PropertyType::biot_coefficient)
+            .template value<double>(variable_array, pos, t, dt);
+
+    double const e_dot = std::get<double>(
+        variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
+
+    double const p_dot = std::get<double>(
+        variable_array[static_cast<int>(Variable::phase_pressure_rate)]);
+
+    double const phi = std::get<double>(
+        variable_array[static_cast<int>(Variable::porosity)]);
+
+    double const w = dt * (e_dot + p_dot / K_SR);
+    return (phi + alpha_b * w) / (1 + w);
+}
+
+PropertyDataType PorosityFromMassBalance::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::liquid_saturation) &&
+           "PorosityFromMassBalance::dValue is implemented for "
+           "derivatives with respect to liquid saturation only.");
+    return 0;
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/PorosityFromMassBalance.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 "MaterialLib/MPL/Property.h"
+
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+
+/// Porosity depending on the volumetric strain rate and phase pressure rate.
+///
+/// This property must be a solid phase.
+class PorosityFromMassBalance final : public Property
+{
+private:
+    Phase* _phase = nullptr;
+
+    /// Parameter, which is used by FEM to set the initial porosity value.
+    ParameterLib::Parameter<double> const& _phi0;
+
+public:
+    PorosityFromMassBalance(
+        ParameterLib::Parameter<double> const& initial_porosity)
+        : _phi0(initial_porosity)
+    {
+    }
+
+    void setScale(
+        std::variant<Medium*, Phase*, Component*> scale_pointer) override
+    {
+        if (std::holds_alternative<Phase*>(scale_pointer))
+        {
+            _phase = std::get<Phase*>(scale_pointer);
+            // TODO (naumov) Test phase is solid phase.
+        }
+        else
+        {
+            OGS_FATAL(
+                "The property 'PorosityFromMassBalance' is "
+                "implemented on the 'phase' scales only.");
+        }
+    }
+
+    PropertyDataType initialValue(ParameterLib::SpatialPosition const& pos,
+                                  double const t) const override
+    {
+        return fromVector(_phi0(t, pos));
+    }
+
+    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

MaterialLib/MPL/Properties/Properties.h

@@ -16,6 +16,7 @@
 #include "IdealGasLaw.h"
 #include "LinearProperty.h"
 #include "ParameterProperty.h"
+#include "PorosityFromMassBalance.h"
 #include "RelPermBrooksCorey.h"
 #include "RelPermLiakopoulos.h"
 #include "RelPermVanGenuchten.h"