Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[GeoMechanicsApplication] Clean up retention response #12392

Merged
merged 10 commits into from
May 22, 2024
Merged

[GeoMechanicsApplication] Clean up retention response #12392

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
7ddb7b5
Replaced two calls to `CalculateRetentionResponse` by more specific ones
avdg81 May 17, 2024
565f47e
Removed two redundant `using` statements
avdg81 May 17, 2024
fdf5ca2
Use a local variable for the Bishop coefficient when possible
avdg81 May 17, 2024
2eb434e
Eliminated member `FluidPressure` of struct `ElementVariables`
avdg81 May 17, 2024
63a6b30
Removed member `EffectiveSaturation` of struct `ElementVariables`
avdg81 May 17, 2024
f650827
Replaced a call to `CalculateRetentionResponse` by a more specific one
avdg81 May 17, 2024
17afc1b
Removed member `FluidPressure` from struct `ElementVariables`
avdg81 May 17, 2024
191efc5
Removed member `Density` from struct `ElementVariables`
avdg81 May 17, 2024
9cf87eb
Fixed GCC builds by adding `this->` to `CalculateRetentionResponse`
avdg81 May 21, 2024
85da166
Processed Wijtze Pieter's review comments
avdg81 May 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 26 additions & 29 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -520,15 +520,21 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const
// Compute Np, GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);

this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = Variables.DegreeOfSaturation;
if (rVariable == EFFECTIVE_SATURATION) rOutput[GPoint] = Variables.EffectiveSaturation;
if (rVariable == BISHOP_COEFFICIENT) rOutput[GPoint] = Variables.BishopCoefficient;
if (rVariable == DERIVATIVE_OF_SATURATION)
rOutput[GPoint] = Variables.DerivativeOfSaturation;
if (rVariable == RELATIVE_PERMEABILITY)
rOutput[GPoint] = Variables.RelativePermeability;
RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector));

if (rVariable == DEGREE_OF_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters);
else if (rVariable == EFFECTIVE_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters);
else if (rVariable == BISHOP_COEFFICIENT)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters);
else if (rVariable == DERIVATIVE_OF_SATURATION)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters);
else if (rVariable == RELATIVE_PERMEABILITY)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
}
} else if (rVariable == HYDRAULIC_HEAD) {
const PropertiesType& rProp = this->GetProperties();
Expand Down Expand Up @@ -693,8 +699,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const

const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT);

Vector TotalStressVector(mStressVector[0].size());

const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = CalculateDeformationGradients();
auto strain_vectors = StressStrainUtilities::CalculateStrains(
Expand All @@ -714,16 +718,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const

Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient);

this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

noalias(TotalStressVector) = mStressVector[GPoint];
noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient *
Variables.BishopCoefficient * Variables.FluidPressure * VoigtVector;
const auto fluid_pressure =
GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector);
RetentionParameters.SetFluidPressure(fluid_pressure);
const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters);

if (rOutput[GPoint].size() != TotalStressVector.size())
rOutput[GPoint].resize(TotalStressVector.size(), false);

rOutput[GPoint] = TotalStressVector;
rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient *
bishop_coefficient * fluid_pressure * VoigtVector;
}
} else if (rVariable == ENGINEERING_STRAIN_VECTOR) {
// Definition of variables
Expand Down Expand Up @@ -1088,9 +1089,9 @@ std::vector<array_1d<double, TDim>> UPwSmallStrainElement<TDim, TNumNodes>::Calc

GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
Variables.FluidPressure =
GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector);
RetentionParameters.SetFluidPressure(Variables.FluidPressure);

RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector));

Variables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
Expand Down Expand Up @@ -1189,7 +1190,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeElementVariables(ElementV
rVariables.UVoigtMatrix.resize(TNumNodes * TDim, VoigtSize, false);

// Retention law
rVariables.FluidPressure = 0.0;
rVariables.DegreeOfSaturation = 1.0;
rVariables.DerivativeOfSaturation = 0.0;
rVariables.RelativePermeability = 1.0;
Expand Down Expand Up @@ -1628,18 +1628,15 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateRetentionResponse(ElementV
{
KRATOS_TRY

rVariables.FluidPressure =
GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector);
rRetentionParameters.SetFluidPressure(rVariables.FluidPressure);
rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
rVariables.Np, rVariables.PressureVector));

rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters);
rVariables.DerivativeOfSaturation =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters);
rVariables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters);
rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters);
rVariables.EffectiveSaturation =
mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(rRetentionParameters);

KRATOS_CATCH("")
}
Expand Down
2 changes: 0 additions & 2 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,12 +179,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
GeometryType::ShapeFunctionsGradientsType DN_DXContainer;

/// Retention Law parameters
double FluidPressure;
double DegreeOfSaturation;
double DerivativeOfSaturation;
double RelativePermeability;
double BishopCoefficient;
double EffectiveSaturation;

// needed for updated Lagrangian:
double detJ;
Expand Down
45 changes: 18 additions & 27 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -923,20 +923,19 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable
// Compute Np, GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);

Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector);
RetentionParameters.SetFluidPressure(Variables.FluidPressure);
RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector));

if (rVariable == DEGREE_OF_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters);
if (rVariable == EFFECTIVE_SATURATION)
else if (rVariable == EFFECTIVE_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters);
if (rVariable == BISHOP_COEFFICIENT)
else if (rVariable == BISHOP_COEFFICIENT)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters);
if (rVariable == DERIVATIVE_OF_SATURATION)
else if (rVariable == DERIVATIVE_OF_SATURATION)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters);
if (rVariable == RELATIVE_PERMEABILITY)
else if (rVariable == RELATIVE_PERMEABILITY)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
}
Expand Down Expand Up @@ -1030,9 +1029,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable
BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++];
}

Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector);
RetentionParameters.SetFluidPressure(Variables.FluidPressure);
RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector));

const double RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
Expand Down Expand Up @@ -1143,8 +1141,6 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable

const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT);

Vector TotalStressVector(mStressVector[0].size());

const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer);
const auto deformation_gradients = CalculateDeformationGradients();
auto strain_vectors = StressStrainUtilities::CalculateStrains(
Expand All @@ -1164,16 +1160,13 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable
Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];
Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient);

this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

noalias(TotalStressVector) = mStressVector[GPoint];
noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient *
Variables.BishopCoefficient * Variables.FluidPressure * VoigtVector;
const auto fluid_pressure =
GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector);
RetentionParameters.SetFluidPressure(fluid_pressure);
const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters);

if (rOutput[GPoint].size() != TotalStressVector.size()) {
rOutput[GPoint].resize(TotalStressVector.size(), false);
}
rOutput[GPoint] = TotalStressVector;
rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient *
bishop_coefficient * fluid_pressure * VoigtVector;
}
} else {
for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i)
Expand Down Expand Up @@ -1467,7 +1460,6 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables
rVariables.DtPressureCoefficient = rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT];

// Retention law
rVariables.FluidPressure = 0.0;
rVariables.DegreeOfSaturation = 1.0;
rVariables.DerivativeOfSaturation = 0.0;
rVariables.RelativePermeability = 1.0;
Expand Down Expand Up @@ -1830,7 +1822,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi
const SizeType Dim = rGeom.WorkingSpaceDimension();
const SizeType NumUNodes = rGeom.PointsNumber();

rVariables.Density = GeoTransportEquationUtilities::CalculateSoilDensity(
const auto soil_density = GeoTransportEquationUtilities::CalculateSoilDensity(
rVariables.DegreeOfSaturation, this->GetProperties());

Vector BodyAcceleration = ZeroVector(Dim);
Expand All @@ -1845,7 +1837,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi
Index = i * Dim;
for (SizeType idim = 0; idim < Dim; ++idim) {
rRightHandSideVector[Index + idim] +=
rVariables.Nu[i] * rVariables.Density * BodyAcceleration[idim] *
rVariables.Nu[i] * soil_density * BodyAcceleration[idim] *
rVariables.IntegrationCoefficientInitialConfiguration;
}
}
Expand Down Expand Up @@ -2052,9 +2044,8 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables
{
KRATOS_TRY

rVariables.FluidPressure =
GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector);
rRetentionParameters.SetFluidPressure(rVariables.FluidPressure);
rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
rVariables.Np, rVariables.PressureVector));

rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters);
rVariables.DerivativeOfSaturation =
Expand Down
2 changes: 0 additions & 2 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,12 +193,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
Vector PressureDtVector;

/// Retention Law parameters
double FluidPressure;
double DegreeOfSaturation;
double DerivativeOfSaturation;
double RelativePermeability;
double BishopCoefficient;
double Density;

// Properties and processinfo variables
bool IgnoreUndrained;
Expand Down
2 changes: 1 addition & 1 deletion applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,7 @@ void SteadyStatePwElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Expand Down
1 change: 0 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwElement : public Transi

/// The definition of the sizetype
using SizeType = std::size_t;
using BaseType::CalculateRetentionResponse;
using BaseType::mConstitutiveLawVector;
using BaseType::mIsInitialised;
using BaseType::mRetentionLawVector;
Expand Down
3 changes: 1 addition & 2 deletions applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -468,7 +468,7 @@ void TransientPwElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftH
GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

Expand Down Expand Up @@ -522,7 +522,6 @@ void TransientPwElement<TDim, TNumNodes>::InitializeElementVariables(ElementVari
rVariables.DN_DXContainer, rVariables.detJContainer, this->GetIntegrationMethod());

// Retention law
rVariables.FluidPressure = 0.0;
rVariables.DegreeOfSaturation = 1.0;
rVariables.DerivativeOfSaturation = 0.0;
rVariables.RelativePermeability = 1.0;
Expand Down
1 change: 0 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwElement : public UPwSmall

/// The definition of the sizetype
using SizeType = std::size_t;
using BaseType::CalculateRetentionResponse;
using BaseType::mConstitutiveLawVector;
using BaseType::mIsInitialised;
using BaseType::mRetentionLawVector;
Expand Down
Loading