Towards element permeability matrices (#12408)

Precalculate relative permeability values as preparation for calculating element permeability matrices. When needed, modify these values by accounting for the permeability update factors (i.e. to account for volume changes in geometrically linear elements).

Other changes include:
- Members `CalculateRetentionResponse` no longer compute the relative permeability. It is now pre-calculated.
- Extracted members that calculate the relative permeability values for all integration points of an element.
- Moved the calculation of permeability update factors to the transport equation utilities.
- Extracted a utility function that calculates fluid pressures at all integration points of an element.
- Inlined the wrapper member functions that calculate and add the permeability matrix, to avoid needless levels of indirection.
- Removed member `PermeabilityUpdateFactor` from two `ElementVariables` data structures.
- Removed the unused `ProcessInfo` data member from the retention laws.
- Removed several comments that had no additional value.
- Miscellaneous cleanup:
  * Reduced the scope of some variables.
  * Don't use `noalias` if no performance benefit is to be expected.
  * Prefer to use `this->` to refer to base members rather than a `using` statement.
  * Removed some unnecessary `KRATOS_TRY` and `KRATOS_CATCH` statements.
  * Removed several empty statements.
  * Unnamed some unused function parameters.

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp

@@ -26,12 +26,11 @@ GeoThermalDispersionLaw::GeoThermalDispersionLaw(std::size_t NumberOfDimensions)
         << "Got invalid number of dimensions: " << mNumberOfDimensions << std::endl;
 }
 
-Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Properties& rProp,
-                                                                 const ProcessInfo& rProcessInfo) const
+Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Properties& rProp) const
 {
     Matrix result = ZeroMatrix(mNumberOfDimensions, mNumberOfDimensions);
 
-    RetentionLaw::Parameters parameters(rProp, rProcessInfo);
+    RetentionLaw::Parameters parameters(rProp);
     auto                     retention_law  = RetentionLawFactory::Clone(rProp);
     const double             saturation     = retention_law->CalculateSaturation(parameters);
     const double             water_fraction = rProp[POROSITY] * saturation;

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h

@@ -30,7 +30,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalDispersionLaw : public Geo
 
     explicit GeoThermalDispersionLaw(SizeType NumberOfDimensions);
 
-    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const override;
+    [[nodiscard]] Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const override;
 
 private:
     std::size_t mNumberOfDimensions = 2;

applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp

@@ -17,7 +17,7 @@
 namespace Kratos
 {
 
-Matrix GeoThermalFilterLaw::CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const
+Matrix GeoThermalFilterLaw::CalculateThermalDispersionMatrix(const Properties& rProp) const
 {
     return ScalarMatrix(1, 1, rProp[THERMAL_CONDUCTIVITY_WATER]);
 }

applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h

@@ -23,7 +23,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalFilterLaw : public GeoTher
 public:
     KRATOS_CLASS_POINTER_DEFINITION(GeoThermalFilterLaw);
 
-    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const override;
+    [[nodiscard]] Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const override;
 
 private:
     friend class Serializer;

applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h

@@ -27,8 +27,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalLaw
 
     virtual ~GeoThermalLaw() = default;
 
-    virtual Matrix CalculateThermalDispersionMatrix(const Properties&  rProp,
-                                                    const ProcessInfo& rProcessInfo) const = 0;
+    [[nodiscard]] virtual Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const = 0;
 
 private:
     friend class Serializer;

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp

@@ -446,8 +446,7 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
     FICElementVariables FICVariables;
     this->InitializeFICElementVariables(FICVariables, Variables.DN_DXContainer, Geom, Prop, CurrentProcessInfo);
 
-    // create general parameters of retention law
-    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo);
+    RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
     const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
     const auto integration_coefficients =
@@ -462,6 +461,8 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
                                          strain_vectors, mStressVector, constitutive_matrices);
     const auto biot_coefficients =
         GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop);
+    const auto relative_permeability_values = this->CalculateRelativePermeabilityValues(
+        GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector));
 
     for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
         this->CalculateKinematics(Variables, GPoint);
@@ -476,6 +477,7 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
 
         this->CalculateShapeFunctionsSecondOrderGradients(FICVariables, Variables);
         this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
+        Variables.RelativePermeability = relative_permeability_values[GPoint];
 
         FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix);
 

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

@@ -149,8 +149,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
     ElementVariables Variables;
     this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
-    // Create general parameters of retention law
-    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+    RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
     const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
     const auto deformation_gradients = CalculateDeformationGradients();
@@ -304,8 +303,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::FinalizeSolutionStep(const ProcessI
     ElementVariables Variables;
     this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
-    // Create general parameters of retention law
-    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+    RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
     const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
     const auto deformation_gradients = CalculateDeformationGradients();
@@ -517,8 +515,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const
         ElementVariables Variables;
         this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
-        // create general parameters of retention law
-        RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+        RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
         for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
             // Compute Np, GradNpT, B and StrainVector
@@ -641,12 +638,10 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(
             deformation_gradients, b_matrices, Variables.DisplacementVector,
             Variables.UseHenckyStrain, this->GetStressStatePolicy().GetVoigtSize());
 
-        std::vector<double> permeability_update_factors;
-        for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
-            permeability_update_factors.push_back(this->CalculatePermeabilityUpdateFactor(strain_vectors[GPoint]));
-        }
-
-        const auto fluid_fluxes = CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo);
+        const auto fluid_fluxes =
+            CalculateFluidFluxes(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors(
+                                     strain_vectors, this->GetProperties()),
+                                 rCurrentProcessInfo);
         for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
             GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]);
         }
@@ -694,8 +689,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const
         ElementVariables Variables;
         this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
-        // Create general parameters of retention law
-        RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+        RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
         const Vector& VoigtVector = this->GetStressStatePolicy().GetVoigtVector();
 
@@ -920,8 +914,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMassMatrix(MatrixType& rMa
     const auto N_container = r_geom.ShapeFunctionsValues(integration_method);
 
     const auto degrees_saturation = GeoTransportEquationUtilities::CalculateDegreesSaturation(
-        this->GetPressureSolutionVector(), N_container, mRetentionLawVector, this->GetProperties(),
-        rCurrentProcessInfo);
+        this->GetPressureSolutionVector(), N_container, mRetentionLawVector, this->GetProperties());
 
     const auto solid_densities =
         GeoTransportEquationUtilities::CalculateSoilDensities(degrees_saturation, this->GetProperties());
@@ -968,8 +961,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType&        rLe
     ElementVariables Variables;
     this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
-    // Create general parameters of retention law
-    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+    RetentionLaw::Parameters RetentionParameters(this->GetProperties());
 
     const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
     const auto integration_coefficients =
@@ -984,6 +976,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType&        rLe
                                          strain_vectors, mStressVector, constitutive_matrices);
     const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(
         constitutive_matrices, this->GetProperties());
+    auto relative_permeability_values = this->CalculateRelativePermeabilityValues(
+        GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector));
+    const auto permeability_update_factors = GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors(
+        strain_vectors, this->GetProperties());
+    std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(),
+                   permeability_update_factors.cbegin(), relative_permeability_values.begin(),
+                   std::multiplies{});
 
     for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
         this->CalculateKinematics(Variables, GPoint);
@@ -997,13 +996,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType&        rLe
         GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
             Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
 
-        CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
+        this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
+        Variables.RelativePermeability = relative_permeability_values[GPoint];
 
         Variables.BiotCoefficient    = biot_coefficients[GPoint];
         Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
             Variables.BiotCoefficient, Variables.DegreeOfSaturation,
             Variables.DerivativeOfSaturation, this->GetProperties());
-        Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector);
 
         Variables.IntegrationCoefficient = integration_coefficients[GPoint];
 
@@ -1034,54 +1033,31 @@ std::vector<array_1d<double, TDim>> UPwSmallStrainElement<TDim, TNumNodes>::Calc
 
     const PropertiesType& rProp = this->GetProperties();
 
-    array_1d<double, TDim> GradPressureTerm;
-
-    // Create general parameters of retention law
-    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);
+    auto relative_permeability_values = this->CalculateRelativePermeabilityValues(
+        GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector));
+    std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(),
+                   rPermeabilityUpdateFactors.cbegin(), relative_permeability_values.begin(),
+                   std::multiplies<>{});
 
     for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
         this->CalculateKinematics(Variables, GPoint);
-        Variables.PermeabilityUpdateFactor = rPermeabilityUpdateFactors[GPoint];
 
         GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
             Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
 
-        RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
-            Variables.Np, Variables.PressureVector));
-
-        Variables.RelativePermeability =
-            mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
+        Variables.RelativePermeability = relative_permeability_values[GPoint];
 
-        noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector);
-        noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration;
+        array_1d<double, TDim> GradPressureTerm = prod(trans(Variables.GradNpT), Variables.PressureVector);
+        GradPressureTerm += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration;
 
         FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse *
-                              Variables.RelativePermeability * Variables.PermeabilityUpdateFactor *
+                              Variables.RelativePermeability *
                               prod(Variables.PermeabilityMatrix, GradPressureTerm));
     }
 
     return FluidFluxes;
 }
 
-template <unsigned int TDim, unsigned int TNumNodes>
-double UPwSmallStrainElement<TDim, TNumNodes>::CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const
-{
-    KRATOS_TRY
-
-    if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) {
-        const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR];
-        const double epsV      = StressStrainUtilities::CalculateTrace(rStrainVector);
-        const double ePrevious = r_prop[POROSITY] / (1.0 - r_prop[POROSITY]);
-        const double eCurrent  = (1.0 + ePrevious) * std::exp(epsV) - 1.0;
-        const double permLog10 = (eCurrent - ePrevious) * InverseCK;
-        return std::pow(10.0, permLog10);
-    }
-
-    return 1.0;
-
-    KRATOS_CATCH("")
-}
-
 template <unsigned int TDim, unsigned int TNumNodes>
 void UPwSmallStrainElement<TDim, TNumNodes>::InitializeElementVariables(ElementVariables& rVariables,
                                                                         const ProcessInfo& rCurrentProcessInfo)
@@ -1172,7 +1148,12 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddLHS(MatrixType& rLef
 
     if (!rVariables.IgnoreUndrained) {
         this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);
-        this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix, rVariables);
+
+        const auto permeability_matrix =
+            GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
+                rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
+                rVariables.RelativePermeability, rVariables.IntegrationCoefficient);
+        GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix);
     }
 
     KRATOS_CATCH("")
@@ -1232,22 +1213,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCompressibilityMatri
     KRATOS_CATCH("")
 }
 
-template <unsigned int TDim, unsigned int TNumNodes>
-void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix,
-                                                                               const ElementVariables& rVariables)
-{
-    KRATOS_TRY
-
-    const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
-        rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
-        rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient);
-
-    // Distribute permeability block matrix into the elemental matrix
-    GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix);
-
-    KRATOS_CATCH("")
-}
-
 template <unsigned int TDim, unsigned int TNumNodes>
 void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddRHS(VectorType& rRightHandSideVector,
                                                                 ElementVariables& rVariables,
@@ -1385,13 +1350,30 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculatePermeabilityFlow(
 
     rPermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
         rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
-        rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient);
+        rVariables.RelativePermeability, rVariables.IntegrationCoefficient);
 
     noalias(rPVector) = -prod(rPermeabilityMatrix, rVariables.PressureVector);
 
     KRATOS_CATCH("")
 }
 
+template <unsigned int TDim, unsigned int TNumNodes>
+std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateRelativePermeabilityValues(
+    const std::vector<double>& rFluidPressures) const
+{
+    KRATOS_ERROR_IF_NOT(rFluidPressures.size() == mRetentionLawVector.size());
+
+    auto retention_law_params = RetentionLaw::Parameters{this->GetProperties()};
+
+    auto result = std::vector<double>{};
+    std::transform(mRetentionLawVector.begin(), mRetentionLawVector.end(), rFluidPressures.begin(),
+                   std::back_inserter(result), [&retention_law_params](auto pRetentionLaw, auto FluidPressure) {
+        retention_law_params.SetFluidPressure(FluidPressure);
+        return pRetentionLaw->CalculateRelativePermeability(retention_law_params);
+    });
+    return result;
+}
+
 template <unsigned int TDim, unsigned int TNumNodes>
 void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector,
                                                                              ElementVariables& rVariables)
@@ -1416,8 +1398,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateFluidBodyFlow(BoundedMatri
     noalias(rPDimMatrix) = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient;
 
     noalias(rPVector) = rVariables.DynamicViscosityInverse * rVariables.FluidDensity *
-                        rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor *
-                        prod(rPDimMatrix, rVariables.BodyAcceleration);
+                        rVariables.RelativePermeability * prod(rPDimMatrix, rVariables.BodyAcceleration);
 
     KRATOS_CATCH("")
 }
@@ -1543,8 +1524,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeProperties(ElementVariabl
     rVariables.Porosity                = rProp[POROSITY];
     GeoElementUtilities::FillPermeabilityMatrix(rVariables.PermeabilityMatrix, rProp);
 
-    rVariables.PermeabilityUpdateFactor = 1.0;
-
     KRATOS_CATCH("")
 }
 
@@ -1579,8 +1558,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateRetentionResponse(ElementV
     rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters);
     rVariables.DerivativeOfSaturation =
         mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters);
-    rVariables.RelativePermeability =
-        mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters);
     rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters);
 
     KRATOS_CATCH("")