Skip to content

Commit

Permalink
Merge pull request #12232 from KratosMultiphysics/geo/extract-utility…
Browse files Browse the repository at this point in the history
…-function-C

[GeoMechanicsApplication] Extract a static utility function for the calculation of the Compressibility Matrix (C)
  • Loading branch information
markelov208 authored Apr 2, 2024
2 parents cac5f87 + 49713bf commit 0c77a49
Show file tree
Hide file tree
Showing 10 changed files with 241 additions and 164 deletions.

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -1271,29 +1271,18 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingMatrix(Matri
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculateCompressibilityMatrix(
BoundedMatrix<double, TNumNodes, TNumNodes>& rPMatrix, const ElementVariables& rVariables) const
{
KRATOS_TRY

noalias(rPMatrix) = -PORE_PRESSURE_SIGN_FACTOR * rVariables.DtPressureCoefficient *
rVariables.BiotModulusInverse * outer_prod(rVariables.Np, rVariables.Np) *
rVariables.IntegrationCoefficient;

KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix,
ElementVariables& rVariables)
{
KRATOS_TRY

this->CalculateCompressibilityMatrix(rVariables.PPMatrix, rVariables);
rVariables.PPMatrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

// Distribute compressibility block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix);
GeoElementUtilities::AssemblePPBlockMatrix(
rLeftHandSideMatrix, rVariables.PPMatrix * rVariables.DtPressureCoefficient);

KRATOS_CATCH("")
}
Expand Down Expand Up @@ -1431,8 +1420,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateCompressibilityFlow(
{
KRATOS_TRY

noalias(rPMatrix) = -PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) * rVariables.IntegrationCoefficient;
noalias(rPMatrix) = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

noalias(rPVector) = -prod(rPMatrix, rVariables.DtPressureVector);

Expand Down Expand Up @@ -1464,7 +1453,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculatePermeabilityFlow(
rPermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient);

noalias(rPVector) = -prod(rPermeabilityMatrix, rVariables.PressureVector);

KRATOS_CATCH("")
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -240,9 +240,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);

virtual void CalculateCompressibilityMatrix(BoundedMatrix<double, TNumNodes, TNumNodes>& rPMatrix,
const ElementVariables& rVariables) const;

virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);

virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1905,12 +1905,12 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddCompressibi
{
KRATOS_TRY

noalias(rVariables.PPMatrix) =
-PORE_PRESSURE_SIGN_FACTOR * rVariables.DtPressureCoefficient * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) * rVariables.JointWidth * rVariables.IntegrationCoefficient;
noalias(rVariables.PPMatrix) = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

// Distribute compressibility block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix);
GeoElementUtilities::AssemblePPBlockMatrix(
rLeftHandSideMatrix, rVariables.PPMatrix * rVariables.DtPressureCoefficient * rVariables.JointWidth);

KRATOS_CATCH("")
}
Expand All @@ -1924,7 +1924,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddPermeabilit
rVariables.PPMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.LocalPermeabilityMatrix,
rVariables.RelativePermeability, rVariables.JointWidth, rVariables.IntegrationCoefficient);

// Distribute permeability block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix);

Expand Down Expand Up @@ -2049,11 +2049,11 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddCompressibi
{
KRATOS_TRY

noalias(rVariables.PPMatrix) = -PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) *
rVariables.JointWidth * rVariables.IntegrationCoefficient;
noalias(rVariables.PPMatrix) = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

noalias(rVariables.PVector) = -1.0 * prod(rVariables.PPMatrix, rVariables.DtPressureVector);
noalias(rVariables.PVector) =
-1.0 * prod(rVariables.PPMatrix * rVariables.JointWidth, rVariables.DtPressureVector);

// Distribute compressibility block vector into elemental vector
GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, rVariables.PVector);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1793,12 +1793,12 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCompressibilityMatrix(Matrix
{
KRATOS_TRY

Matrix CompressibilityMatrix =
-PORE_PRESSURE_SIGN_FACTOR * rVariables.DtPressureCoefficient * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) * rVariables.IntegrationCoefficient;
Matrix CompressibilityMatrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

// Distribute compressibility block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, CompressibilityMatrix);
GeoElementUtilities::AssemblePPBlockMatrix(
rLeftHandSideMatrix, CompressibilityMatrix * rVariables.DtPressureCoefficient);

KRATOS_CATCH("")
}
Expand All @@ -1811,7 +1811,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddPermeabilityMatrix(MatrixTyp
Matrix PermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix(
rVariables.DNp_DX, rVariables.DynamicViscosityInverse, rVariables.IntrinsicPermeability,
rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient);

// Distribute permeability block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, PermeabilityMatrix);

Expand Down Expand Up @@ -1939,8 +1939,8 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCompressibilityFlow(VectorTy
{
KRATOS_TRY

Matrix CompressibilityMatrix = -PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) * rVariables.IntegrationCoefficient;
Matrix CompressibilityMatrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

Vector CompressibilityFlow = -prod(CompressibilityMatrix, rVariables.PressureDtVector);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -541,10 +541,11 @@ void TransientPwElement<TDim, TNumNodes>::CalculateAndAddCompressibilityMatrix(M
{
KRATOS_TRY;

this->CalculateCompressibilityMatrix(rVariables.PPMatrix, rVariables);
rVariables.PPMatrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

// Distribute compressibility block matrix into the elemental matrix
rLeftHandSideMatrix += rVariables.PPMatrix;
rLeftHandSideMatrix += (rVariables.PPMatrix * rVariables.DtPressureCoefficient);

KRATOS_CATCH("");
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -581,12 +581,11 @@ void TransientPwInterfaceElement<TDim, TNumNodes>::CalculateAndAddCompressibilit
{
KRATOS_TRY;

noalias(rVariables.PPMatrix) =
-PORE_PRESSURE_SIGN_FACTOR * rVariables.DtPressureCoefficient * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) * rVariables.JointWidth * rVariables.IntegrationCoefficient;
noalias(rVariables.PPMatrix) = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

// Distribute compressibility block matrix into the elemental matrix
rLeftHandSideMatrix += rVariables.PPMatrix;
rLeftHandSideMatrix += (rVariables.PPMatrix * rVariables.DtPressureCoefficient * rVariables.JointWidth);

KRATOS_CATCH("")
}
Expand All @@ -600,7 +599,7 @@ void TransientPwInterfaceElement<TDim, TNumNodes>::CalculateAndAddPermeabilityMa
rVariables.PPMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.LocalPermeabilityMatrix,
rVariables.RelativePermeability, rVariables.JointWidth, rVariables.IntegrationCoefficient);

// Distribute permeability block matrix into the elemental matrix
rLeftHandSideMatrix += rVariables.PPMatrix;

Expand Down Expand Up @@ -629,11 +628,11 @@ void TransientPwInterfaceElement<TDim, TNumNodes>::CalculateAndAddCompressibilit
{
KRATOS_TRY;

noalias(rVariables.PPMatrix) = -PORE_PRESSURE_SIGN_FACTOR * rVariables.BiotModulusInverse *
outer_prod(rVariables.Np, rVariables.Np) *
rVariables.JointWidth * rVariables.IntegrationCoefficient;
noalias(rVariables.PPMatrix) = GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient);

noalias(rVariables.PVector) = -1.0 * prod(rVariables.PPMatrix, rVariables.DtPressureVector);
noalias(rVariables.PVector) =
-1.0 * prod(rVariables.PPMatrix * rVariables.JointWidth, rVariables.DtPressureVector);

// Distribute compressibility block vector into elemental vector
rRightHandSideVector += rVariables.PVector;
Expand Down
10 changes: 10 additions & 0 deletions applications/GeoMechanicsApplication/custom_utilities/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,19 @@ The mathematical definition of the permeability matrix is:
$$H = \int_\Omega (\nabla N_p)^T \frac{1}{\mu} k \nabla N_p d\Omega$$
where $\nabla N_p$ is the gradient of the pressure shape function, $\mu$ is the dynamic viscosity (material parameter) and $k$ is material permeability matrix. The k matrix allows one to take into account, for example, an anisotropic permeability.

### Compressibility matrix (C)

The mathematical definition is:
$$C = \int_\Omega N_{p}^T \frac{1}{Q} N_p d\Omega$$

Where $\Omega$ is the domain, $N_p$ is the pressure shape function and $1/Q$ is the inverse Biot modulus.



File transport_equation_utilities.hpp includes

- CalculatePermeabilityMatrix function
- CalculateCompressibilityMatrix function



Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -49,5 +49,16 @@ class GeoTransportEquationUtilities
RelativePermeability * PermeabilityUpdateFactor * IntegrationCoefficient;
}

template <unsigned int TNumNodes>
static inline BoundedMatrix<double, TNumNodes, TNumNodes> CalculateCompressibilityMatrix(
const Vector& rNp, double BiotModulusInverse, double IntegrationCoefficient)
{
return CalculateCompressibilityMatrix(rNp, BiotModulusInverse, IntegrationCoefficient);
}

static inline Matrix CalculateCompressibilityMatrix(const Vector& rNp, double BiotModulusInverse, double IntegrationCoefficient)
{
return -PORE_PRESSURE_SIGN_FACTOR * BiotModulusInverse * outer_prod(rNp, rNp) * IntegrationCoefficient;
}
}; /* Class GeoTransportEquationUtilities*/
} /* namespace Kratos.*/
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
// KRATOS___
// // ) )
// // ___ ___
// // ____ //___) ) // ) )
// // / / // // / /
// ((____/ / ((____ ((___/ / MECHANICS
//
// License: geo_mechanics_application/license.txt
//
// Main authors: Gennady Markelov
//

#include "containers/model.h"
#include "custom_utilities/transport_equation_utilities.hpp"
#include "includes/checks.h"
#include "testing/testing.h"
#include <boost/numeric/ublas/assignment.hpp>

using namespace Kratos;

namespace Kratos::Testing
{

KRATOS_TEST_CASE_IN_SUITE(Calculatecompressibility_matrix2D3NGivesCorrectResults, KratosGeoMechanicsFastSuite)
{
Vector n_p(3);
n_p <<= 1.0, 2.0, 3.0;

BoundedMatrix<double, 3, 3> compressibility_matrix = ZeroMatrix(3, 3);
const double integration_coefficient = 1.0;
const double biot_modulus_inverse = 0.02;
compressibility_matrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix<3>(
n_p, biot_modulus_inverse, integration_coefficient);

BoundedMatrix<double, 3, 3> expected_compressibility_matrix;
// clang-format off
expected_compressibility_matrix <<= -0.02,-0.04,-0.06,
-0.04,-0.08,-0.12,
-0.06,-0.12,-0.18;
// clang-format on

KRATOS_CHECK_MATRIX_NEAR(compressibility_matrix, expected_compressibility_matrix, 1e-12)
}

KRATOS_TEST_CASE_IN_SUITE(Calculatecompressibility_matrix3D4NGivesCorrectResults, KratosGeoMechanicsFastSuite)
{
Vector n_p(4);
n_p <<= 1.0, 2.0, 3.0, 3.0;

BoundedMatrix<double, 4, 4> compressibility_matrix = ZeroMatrix(4, 4);
BoundedMatrix<double, 3, 3> material_compressibility_matrix;

const double integration_coefficient = 1.0;
const double biot_modulus_inverse = 0.1;
compressibility_matrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix<4>(
n_p, biot_modulus_inverse, integration_coefficient);

BoundedMatrix<double, 4, 4> expected_compressibility_matrix;
// clang-format off
expected_compressibility_matrix <<= -0.1,-0.2,-0.3,-0.3,
-0.2,-0.4,-0.6,-0.6,
-0.3,-0.6,-0.9,-0.9,
-0.3,-0.6,-0.9,-0.9;
// clang-format on

KRATOS_CHECK_MATRIX_NEAR(compressibility_matrix, expected_compressibility_matrix, 1e-12)
}

} // namespace Kratos::Testing

0 comments on commit 0c77a49

Please sign in to comment.