[GeoMechanicsApplication] Implement well constitutive behaviour for line pressure element

applications/GeoMechanicsApplication/custom_elements/filter_compressibility_calculator.cpp

@@ -0,0 +1,107 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   Richard Faasse
+//
+#include "filter_compressibility_calculator.h"
+#include "custom_utilities/transport_equation_utilities.hpp"
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+FilterCompressibilityCalculator::FilterCompressibilityCalculator(InputProvider AnInputProvider)
+    : mInputProvider(std::move(AnInputProvider))
+{
+}
+
+Matrix FilterCompressibilityCalculator::LHSContribution()
+{
+    return LHSContribution(CalculateCompressibilityMatrix());
+}
+
+Vector FilterCompressibilityCalculator::RHSContribution()
+{
+    return RHSContribution(CalculateCompressibilityMatrix());
+}
+
+Vector FilterCompressibilityCalculator::RHSContribution(const Matrix& rCompressibilityMatrix) const
+{
+    return -prod(rCompressibilityMatrix, mInputProvider.GetNodalValues(DT_WATER_PRESSURE));
+}
+
+Matrix FilterCompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const
+{
+    return mInputProvider.GetMatrixScalarFactor() * rCompressibilityMatrix;
+}
+
+std::pair<Matrix, Vector> FilterCompressibilityCalculator::LocalSystemContribution()
+{
+    const auto compressibility_matrix = CalculateCompressibilityMatrix();
+    return {LHSContribution(compressibility_matrix), RHSContribution(compressibility_matrix)};
+}
+
+Matrix FilterCompressibilityCalculator::CalculateCompressibilityMatrix() const
+{
+    const auto& r_N_container            = mInputProvider.GetNContainer();
+    const auto& integration_coefficients = mInputProvider.GetIntegrationCoefficients();
+    const auto& projected_gravity_on_integration_points =
+        mInputProvider.GetProjectedGravityForIntegrationPoints();
+    auto result = Matrix{ZeroMatrix{r_N_container.size2(), r_N_container.size2()}};
+    for (unsigned int integration_point_index = 0;
+         integration_point_index < integration_coefficients.size(); ++integration_point_index) {
+        const auto N = Vector{row(r_N_container, integration_point_index)};
+        result += GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
+            N, CalculateElasticCapacity(projected_gravity_on_integration_points[integration_point_index]),
+            integration_coefficients[integration_point_index]);
+    }
+    return result;
+}
+
+double FilterCompressibilityCalculator::CalculateElasticCapacity(double ProjectedGravity) const
+{
+    const auto&  r_properties = mInputProvider.GetElementProperties();
+    const double result =
+        1.0 / (r_properties[DENSITY_WATER] * ProjectedGravity * r_properties[FILTER_LENGTH]) +
+        1.0 / r_properties[BULK_MODULUS_FLUID];
+    return result;
+}
+
+const Properties& FilterCompressibilityCalculator::InputProvider::GetElementProperties() const
+{
+    return mGetElementProperties();
+}
+
+const Matrix& FilterCompressibilityCalculator::InputProvider::GetNContainer() const
+{
+    return mGetNContainer();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetIntegrationCoefficients() const
+{
+    return mGetIntegrationCoefficients();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetProjectedGravityForIntegrationPoints() const
+{
+    return mGetProjectedGravityForIntegrationPoints();
+}
+
+double FilterCompressibilityCalculator::InputProvider::GetMatrixScalarFactor() const
+{
+    return mGetMatrixScalarFactor();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetNodalValues(const Variable<double>& rVariable) const
+{
+    return mGetNodalValues(rVariable);
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/filter_compressibility_calculator.h

@@ -0,0 +1,78 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   Richard Faasse
+
+#pragma once
+
+#include "contribution_calculator.h"
+#include "custom_retention/retention_law.h"
+#include "includes/properties.h"
+#include "includes/ublas_interface.h"
+
+#include <utility>
+#include <vector>
+
+namespace Kratos
+{
+
+class FilterCompressibilityCalculator : public ContributionCalculator
+{
+public:
+    class InputProvider
+    {
+    public:
+        InputProvider(std::function<const Properties&()> GetElementProperties,
+                      std::function<const Matrix&()>     GetNContainer,
+                      std::function<Vector()>            GetIntegrationCoefficients,
+                      std::function<Vector()>            GetProjectedGravityForIntegrationPoints,
+                      std::function<double()>            GetMatrixScalarFactor,
+                      std::function<Vector(const Variable<double>&)> GetNodalValuesOf)
+            : mGetElementProperties(std::move(GetElementProperties)),
+              mGetNContainer(std::move(GetNContainer)),
+              mGetIntegrationCoefficients(std::move(GetIntegrationCoefficients)),
+              mGetProjectedGravityForIntegrationPoints(std::move(GetProjectedGravityForIntegrationPoints)),
+              mGetMatrixScalarFactor(std::move(GetMatrixScalarFactor)),
+              mGetNodalValues(std::move(GetNodalValuesOf))
+        {
+        }
+
+        [[nodiscard]] const Properties& GetElementProperties() const;
+        [[nodiscard]] const Matrix&     GetNContainer() const;
+        [[nodiscard]] Vector            GetIntegrationCoefficients() const;
+        [[nodiscard]] Vector            GetProjectedGravityForIntegrationPoints() const;
+        [[nodiscard]] double            GetMatrixScalarFactor() const;
+        [[nodiscard]] Vector            GetNodalValues(const Variable<double>& rVariable) const;
+
+    private:
+        std::function<const Properties&()>             mGetElementProperties;
+        std::function<const Matrix&()>                 mGetNContainer;
+        std::function<Vector()>                        mGetIntegrationCoefficients;
+        std::function<Vector()>                        mGetProjectedGravityForIntegrationPoints;
+        std::function<double()>                        mGetMatrixScalarFactor;
+        std::function<Vector(const Variable<double>&)> mGetNodalValues;
+    };
+
+    explicit FilterCompressibilityCalculator(InputProvider AnInputProvider);
+
+    Matrix                    LHSContribution() override;
+    Vector                    RHSContribution() override;
+    std::pair<Matrix, Vector> LocalSystemContribution() override;
+
+private:
+    [[nodiscard]] Matrix CalculateCompressibilityMatrix() const;
+    [[nodiscard]] double CalculateElasticCapacity(double ProjectedGravity) const;
+    [[nodiscard]] Vector RHSContribution(const Matrix& rCompressibilityMatrix) const;
+    [[nodiscard]] Matrix LHSContribution(const Matrix& rCompressibilityMatrix) const;
+
+    InputProvider mInputProvider;
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h

@@ -18,6 +18,7 @@
 #include "custom_retention/retention_law_factory.h"
 #include "custom_utilities/dof_utilities.h"
 #include "custom_utilities/element_utilities.hpp"
+#include "filter_compressibility_calculator.h"
 #include "geo_mechanics_application_variables.h"
 #include "includes/cfd_variables.h"
 #include "includes/element.h"
@@ -252,57 +253,70 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
 
         const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container);
 
-        const std::size_t number_integration_points =
-            GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
-        GeometryType::JacobiansType J_container;
-        J_container.resize(number_integration_points, false);
-        for (std::size_t i = 0; i < number_integration_points; ++i) {
-            J_container[i].resize(GetGeometry().WorkingSpaceDimension(),
-                                  GetGeometry().LocalSpaceDimension(), false);
-        }
-        GetGeometry().Jacobian(J_container, this->GetIntegrationMethod());
-
         const auto&                 r_properties = GetProperties();
         BoundedMatrix<double, 1, 1> constitutive_matrix;
         GeoElementUtilities::FillPermeabilityMatrix(constitutive_matrix, r_properties);
 
         RetentionLaw::Parameters RetentionParameters(GetProperties());
 
-        array_1d<double, TNumNodes * TDim> volume_acceleration;
-        GeoElementUtilities::GetNodalVariableVector<TDim, TNumNodes>(
-            volume_acceleration, GetGeometry(), VOLUME_ACCELERATION);
+        const auto projected_gravity = CalculateProjectedGravityAtIntegrationPoints(N_container);
 
         array_1d<double, TNumNodes> fluid_body_vector = ZeroVector(TNumNodes);
         for (unsigned int integration_point_index = 0;
              integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
              ++integration_point_index) {
-            array_1d<double, TDim> body_acceleration;
-            GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
-                body_acceleration, N_container, volume_acceleration, integration_point_index);
-
-            array_1d<double, TDim> tangent_vector = column(J_container[integration_point_index], 0);
-            tangent_vector /= norm_2(tangent_vector);
-
-            array_1d<double, 1> projected_gravity = ZeroVector(1);
-            projected_gravity(0) = MathUtils<double>::Dot(tangent_vector, body_acceleration);
-            const auto N         = Vector{row(N_container, integration_point_index)};
+            const auto N = Vector{row(N_container, integration_point_index)};
             double     RelativePermeability =
                 mRetentionLawVector[integration_point_index]->CalculateRelativePermeability(RetentionParameters);
             fluid_body_vector +=
                 r_properties[DENSITY_WATER] * RelativePermeability *
-                prod(prod(shape_function_gradients[integration_point_index], constitutive_matrix), projected_gravity) *
+                prod(prod(shape_function_gradients[integration_point_index], constitutive_matrix),
+                     ScalarVector(1, projected_gravity[integration_point_index])) *
                 integration_coefficients[integration_point_index] / r_properties[DYNAMIC_VISCOSITY];
         }
         return fluid_body_vector;
     }
 
+    Vector CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) const
+    {
+        const auto number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
+        GeometryType::JacobiansType J_container{number_integration_points};
+        for (std::size_t i = 0; i < number_integration_points; ++i) {
+            J_container[i].resize(GetGeometry().WorkingSpaceDimension(),
+                                  GetGeometry().LocalSpaceDimension(), false);
+        }
+        GetGeometry().Jacobian(J_container, this->GetIntegrationMethod());
+
+        array_1d<double, TNumNodes * TDim> volume_acceleration;
+        GeoElementUtilities::GetNodalVariableVector<TDim, TNumNodes>(
+            volume_acceleration, GetGeometry(), VOLUME_ACCELERATION);
+        array_1d<double, TDim> body_acceleration;
+
+        Vector projected_gravity = ZeroVector(number_integration_points);
+
+        for (unsigned int integration_point_index = 0;
+             integration_point_index < number_integration_points; ++integration_point_index) {
+            GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
+                body_acceleration, rNContainer, volume_acceleration, integration_point_index);
+            array_1d<double, TDim> tangent_vector = column(J_container[integration_point_index], 0);
+            tangent_vector /= norm_2(tangent_vector);
+            projected_gravity(integration_point_index) =
+                MathUtils<>::Dot(tangent_vector, body_acceleration);
+        }
+        return projected_gravity;
+    }
+
     std::unique_ptr<ContributionCalculator> CreateCalculator(const CalculationContribution& rContribution,
                                                              const ProcessInfo& rCurrentProcessInfo)
     {
         switch (rContribution) {
         case CalculationContribution::Permeability:
             return std::make_unique<PermeabilityCalculator>(CreatePermeabilityInputProvider());
         case CalculationContribution::Compressibility:
+            if (GetProperties()[RETENTION_LAW] == "PressureFilterLaw") {
+                return std::make_unique<FilterCompressibilityCalculator>(
+                    CreateFilterCompressibilityInputProvider(rCurrentProcessInfo));
+            }
             return std::make_unique<CompressibilityCalculator>(CreateCompressibilityInputProvider(rCurrentProcessInfo));
         default:
             KRATOS_ERROR << "Unknown contribution" << std::endl;
@@ -317,6 +331,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
             MakeNodalVariableGetter());
     }
 
+    FilterCompressibilityCalculator::InputProvider CreateFilterCompressibilityInputProvider(const ProcessInfo& rCurrentProcessInfo)
+    {
+        return FilterCompressibilityCalculator::InputProvider(
+            MakePropertiesGetter(), MakeNContainerGetter(), MakeIntegrationCoefficientsGetter(),
+            MakeProjectedGravityForIntegrationPointsGetter(),
+            MakeMatrixScalarFactorGetter(rCurrentProcessInfo), MakeNodalVariableGetter());
+    }
+
     PermeabilityCalculator::InputProvider CreatePermeabilityInputProvider()
     {
         return PermeabilityCalculator::InputProvider(
@@ -350,6 +372,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
         };
     }
 
+    auto MakeProjectedGravityForIntegrationPointsGetter()
+    {
+        return [this]() -> Vector {
+            return CalculateProjectedGravityAtIntegrationPoints(
+                GetGeometry().ShapeFunctionsValues(GetIntegrationMethod()));
+        };
+    }
+
     static auto MakeMatrixScalarFactorGetter(const ProcessInfo& rCurrentProcessInfo)
     {
         return [&rCurrentProcessInfo]() { return rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; };

applications/GeoMechanicsApplication/custom_python/geo_mechanics_python_application.cpp

@@ -145,6 +145,8 @@ PYBIND11_MODULE(KratosGeoMechanicsApplication, m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, UMAT_PARAMETERS)
 
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DT_TEMPERATURE_COEFFICIENT)
+
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FILTER_LENGTH)
 }
 
 } // namespace Kratos::Python.

applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h

@@ -46,20 +46,22 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLawFactory
     {
         if (rMaterialProperties.Has(RETENTION_LAW)) {
             const std::string& RetentionLawName = rMaterialProperties[RETENTION_LAW];
-            if (RetentionLawName == "VanGenuchtenLaw") return make_unique<VanGenuchtenLaw>();
+            if (RetentionLawName == "VanGenuchtenLaw") return std::make_unique<VanGenuchtenLaw>();
 
-            if (RetentionLawName == "SaturatedLaw") return make_unique<SaturatedLaw>();
+            if (RetentionLawName == "SaturatedLaw") return std::make_unique<SaturatedLaw>();
 
             if (RetentionLawName == "SaturatedBelowPhreaticLevelLaw")
-                return make_unique<SaturatedBelowPhreaticLevelLaw>();
+                return std::make_unique<SaturatedBelowPhreaticLevelLaw>();
+
+            if (RetentionLawName == "PressureFilterLaw") return std::make_unique<SaturatedLaw>();
 
             KRATOS_ERROR << "Undefined RETENTION_LAW! " << RetentionLawName << std::endl;
 
             return nullptr;
         }
 
         // default is saturated law
-        return make_unique<SaturatedLaw>();
+        return std::make_unique<SaturatedLaw>();
     }
 
 }; // Class RetentionLawFactory

applications/GeoMechanicsApplication/custom_utilities/element_utilities.cpp

@@ -35,7 +35,12 @@ void GeoElementUtilities::FillPermeabilityMatrix(BoundedMatrix<double, 1, 1>&
                                                  const Element::PropertiesType& Prop)
 {
     // 1D
-    rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX];
+    if (Prop[RETENTION_LAW] == "PressureFilterLaw") {
+        const double equivalent_radius_square = Prop[CROSS_AREA] / Globals::Pi;
+        rPermeabilityMatrix(0, 0)             = equivalent_radius_square * 0.125;
+    } else {
+        rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX];
+    }
 }
 
 void GeoElementUtilities::FillPermeabilityMatrix(BoundedMatrix<double, 2, 2>&   rPermeabilityMatrix,

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -528,6 +528,8 @@ void KratosGeoMechanicsApplication::Register()
     KRATOS_REGISTER_VARIABLE(PIPE_EROSION)
     KRATOS_REGISTER_VARIABLE(PIPE_ACTIVE)
 
+    KRATOS_REGISTER_VARIABLE(FILTER_LENGTH)
+
     // UDSM
     KRATOS_REGISTER_VARIABLE(UDSM_NAME) // Also for UMAT
     KRATOS_REGISTER_VARIABLE(UDSM_NUMBER)

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

@@ -173,6 +173,8 @@ KRATOS_CREATE_VARIABLE(bool, PIPE_EROSION)
 KRATOS_CREATE_VARIABLE(bool, PIPE_IN_EQUILIBRIUM)
 KRATOS_CREATE_VARIABLE(bool, PIPE_ACTIVE)
 
+KRATOS_CREATE_VARIABLE(double, FILTER_LENGTH)
+
 // UDSM
 KRATOS_CREATE_VARIABLE(std::string, UDSM_NAME) // Also for UMAT
 KRATOS_CREATE_VARIABLE(int, UDSM_NUMBER)

applications/GeoMechanicsApplication/geo_mechanics_application_variables.h

@@ -190,6 +190,9 @@ KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, DIFF_PIPE_
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, bool, PIPE_EROSION)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, bool, PIPE_ACTIVE)
 
+// Filter
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, FILTER_LENGTH)
+
 // UDSM
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, std::string, UDSM_NAME) // also for umat
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, int, UDSM_NUMBER)