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] Partially saturated flow does not work as expected #12914

Merged
merged 53 commits into from
Dec 21, 2024
Merged

[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
53 commits
Select commit Hold shift + click to select a range
952395e
Avoid zero divisions + a lot of printing
WPK4FEM Oct 15, 2024
1118c7c
Removed an unused variable and some commented out code
avdg81 Oct 25, 2024
b7aa7d8
Fix effective saturation. Commented print statements.
WPK4FEM Dec 6, 2024
5f7b856
start of unit tests for SaturatedBelowPhreaticLevelLaw
WPK4FEM Dec 6, 2024
3412224
Improved error messages for saturated below phreatic level and unit t…
WPK4FEM Dec 9, 2024
684a1c1
removed unwanted include
WPK4FEM Dec 9, 2024
bbd8a8a
ternary operator
WPK4FEM Dec 9, 2024
8397ff0
Unit test for saturated law
WPK4FEM Dec 9, 2024
cebc4ab
ternary operator
WPK4FEM Dec 9, 2024
6112a0f
Checks and unit tests for Van Genuchten added.
WPK4FEM Dec 9, 2024
ec76360
Fixes for partially saturated flow
mnabideltares Dec 9, 2024
04f9aea
Attempt to only place includes there where they are used.
WPK4FEM Dec 9, 2024
1cf2ff6
Cleanup attempts
WPK4FEM Dec 9, 2024
9def637
Trying to get rid of empty functions.
WPK4FEM Dec 9, 2024
4029c99
Remove added print statements from core
WPK4FEM Dec 9, 2024
23a5739
sonar cloud remarks and naming style corrections
WPK4FEM Dec 10, 2024
24f064e
Coupling matrix computation without potential zero division also in s…
WPK4FEM Dec 10, 2024
9594252
Also avoid zero division in RHS setup of coupling terms.
WPK4FEM Dec 10, 2024
09f4711
Added test cases for Pw 2D3N and 2D6N, with documentation
mnabideltares Dec 12, 2024
1535da2
van Genuchten gl >= 0, clean unused parameters from tests, some geome…
WPK4FEM Dec 11, 2024
231f591
IgnoreUndrained setup of LHS and RHS made symmetric.
WPK4FEM Dec 11, 2024
a3ee3e2
No water solution intended in these tests, made IGNORE_UNDRAINED cons…
WPK4FEM Dec 12, 2024
c53e206
fix for Readme
mnabideltares Dec 12, 2024
eab53df
A test xase for U-Pw diff order 2D6N is added
mnabideltares Dec 12, 2024
0ea3c26
2 tests are added for U-Pw small strain element
mnabideltares Dec 12, 2024
882574b
fix of revision 60eb250069fc89bd84a03ad2f40f113fcfcc51dc changed resu…
WPK4FEM Dec 12, 2024
c43dc78
fixes
mnabideltares Dec 12, 2024
443f675
fix
mnabideltares Dec 12, 2024
13ae2a5
Sonar Cloud unused stuff complaint resolved
WPK4FEM Dec 16, 2024
9751e1c
Removed unused retention parameters ( review Gennady )
WPK4FEM Dec 16, 2024
f38f8a0
std::transform for pressure getting ( review comment of Gennady )
WPK4FEM Dec 16, 2024
23fc002
Removed unused K0 things ( review comment Gennady )
WPK4FEM Dec 16, 2024
0559d07
Renamed kolom to the english column ( review comment Gennady )
WPK4FEM Dec 16, 2024
00e783e
avoid code repetition
WPK4FEM Dec 16, 2024
cbe888c
kp format on request of Anne
WPK4FEM Dec 16, 2024
deefa26
Set material parameters in a common directory, to avoid file duplicat…
mnabideltares Dec 16, 2024
4464da3
Review comments Anne, rule of zero. Some typografy and typo
WPK4FEM Dec 19, 2024
97f5e48
std::pow
WPK4FEM Dec 19, 2024
634b39d
style guide review comment by Anne
WPK4FEM Dec 19, 2024
4d3c153
consistency between test and text. Review Anne.
WPK4FEM Dec 19, 2024
17868b0
Text of test should match error text
WPK4FEM Dec 19, 2024
08ba4bd
const function
WPK4FEM Dec 19, 2024
452ff63
style guide issues from Review by Anne.
WPK4FEM Dec 19, 2024
e6b98c6
Removed helper in saturated below phreatic level testing.
WPK4FEM Dec 19, 2024
caf53fd
Removed helper function from saturatedLaw testing.
WPK4FEM Dec 19, 2024
ee04f2d
Removed helper function from VanGenuchtenLaw testing.
WPK4FEM Dec 19, 2024
38ceb9d
Review comments for integration testing script.
WPK4FEM Dec 19, 2024
1d9c427
README improvement and no rotation D.o.F.
WPK4FEM Dec 19, 2024
28e7cc9
Removed output of values not checked in tests.
WPK4FEM Dec 19, 2024
702fc44
Remove extrapolation process, its results are not used.
WPK4FEM Dec 19, 2024
7f6c55d
No displacement model parts needed for Pw elements
WPK4FEM Dec 19, 2024
a8cd4ca
Create members small_strain_U_Pw_diff_orders, kp format, coords in py…
WPK4FEM Dec 20, 2024
a741cd1
Merge remote-tracking branch 'origin/master' into geo/partial_saturat…
WPK4FEM Dec 20, 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
11 changes: 4 additions & 7 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,10 +135,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
}

mRetentionLawVector.resize(number_of_integration_points);
for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) {
mRetentionLawVector[i] = RetentionLawFactory::Clone(r_properties);
mRetentionLawVector[i]->InitializeMaterial(
markelov208 marked this conversation as resolved.
Show resolved Hide resolved
r_properties, r_geometry, row(r_geometry.ShapeFunctionsValues(mThisIntegrationMethod), i));
for (auto& r_retention_law : mRetentionLawVector) {
r_retention_law = RetentionLawFactory::Clone(r_properties);
}

if (mStressVector.size() != number_of_integration_points) {
Expand All @@ -151,9 +149,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo)

mStateVariablesFinalized.resize(number_of_integration_points);
for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) {
int nStateVariables = 0;
nStateVariables = mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables);
if (nStateVariables > 0) {
if (auto nStateVariables = 0;
mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables) > 0) {
mConstitutiveLawVector[i]->SetValue(STATE_VARIABLES, mStateVariablesFinalized[i], rCurrentProcessInfo);
}
}
Expand Down
75 changes: 30 additions & 45 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 @@ -258,8 +258,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
Expand All @@ -283,8 +281,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
noalias(Variables.StressVector) = mStressVector[GPoint];
ConstitutiveParameters.SetStressVector(Variables.StressVector);
mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters);

mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters);
markelov208 marked this conversation as resolved.
Show resolved Hide resolved
}

// Reset hydraulic discharge
Expand All @@ -299,9 +295,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::ResetHydraulicDischarge()
KRATOS_TRY

// Reset hydraulic discharge
GeometryType& rGeom = this->GetGeometry();
for (unsigned int i = 0; i < TNumNodes; ++i) {
ThreadSafeNodeWrite(rGeom[i], HYDRAULIC_DISCHARGE, 0.0);
for (auto& r_node : this->GetGeometry()) {
ThreadSafeNodeWrite(r_node, HYDRAULIC_DISCHARGE, 0.0);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -357,8 +352,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeNonLinearIteration(const
{
KRATOS_TRY

const GeometryType& rGeom = this->GetGeometry();
ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, this->GetProperties(), rCurrentProcessInfo);
ConstitutiveLaw::Parameters ConstitutiveParameters(this->GetGeometry(), this->GetProperties(),
rCurrentProcessInfo);
ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_STRESS);
ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN);
ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); // Note: this is for nonlocal damage
Expand Down Expand Up @@ -403,8 +398,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::FinalizeSolutionStep(const ProcessI
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
Expand Down Expand Up @@ -434,8 +427,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::FinalizeSolutionStep(const ProcessI
mStateVariablesFinalized[GPoint] =
mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]);

mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);

if (rCurrentProcessInfo[NODAL_SMOOTHING])
this->SaveGPStress(StressContainer, mStressVector[GPoint], GPoint);
}
Expand Down Expand Up @@ -622,19 +613,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const

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);
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue(
RetentionParameters, rVariable, rOutput[GPoint]);
}
} else if (rVariable == HYDRAULIC_HEAD) {
const PropertiesType& rProp = this->GetProperties();
Expand Down Expand Up @@ -1246,16 +1226,17 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddLHS(MatrixType& rLef
KRATOS_TRY

this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix, rVariables);
this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables);

if (!rVariables.IgnoreUndrained) {
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);

if (!rVariables.IgnoreUndrained) {
const auto permeability_matrix =
GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
rVariables.RelativePermeability, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix);

this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables);
}

KRATOS_CATCH("")
Expand All @@ -1280,17 +1261,18 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingMatrix(Matri
{
KRATOS_TRY

const auto coupling_matrix = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
const Matrix coupling_matrix_up = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix);
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix_up);
Comment on lines +1264 to +1267
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am confused here: I had expected a minus sign in front of the coupling matrix, but I haven't found it yet... Perhaps @WPK4FEM can shed some light here, please?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As we discussed, I have the same doubt and have been looking for the same. Trying to change it leads at least to some tests failing. We need to sort this out building up a really validated test suite.


if (!rVariables.IgnoreUndrained) {
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
const BoundedMatrix<double, TNumNodes, TNumNodes * TDim> transposed_coupling_matrix =
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * rVariables.VelocityCoefficient *
trans(coupling_matrix);
GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, transposed_coupling_matrix);
const Matrix coupling_matrix_pu = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssemblePUBlockMatrix(
rLeftHandSideMatrix,
PORE_PRESSURE_SIGN_FACTOR * rVariables.VelocityCoefficient * trans(coupling_matrix_pu));
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -1389,15 +1371,18 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingTerms(Vector
rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient);

const array_1d<double, TNumNodes * TDim> coupling_terms = prod(coupling_matrix, rVariables.PressureVector);
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_terms);
const array_1d<double, TNumNodes * TDim> coupling_force = prod(coupling_matrix, rVariables.PressureVector);
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_force);

if (!rVariables.IgnoreUndrained) {
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
const array_1d<double, TNumNodes> transposed_coupling_matrix =
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient *
prod(trans(coupling_matrix), rVariables.VelocityVector);
GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, transposed_coupling_matrix);
const Matrix p_coupling_matrix =
(-1.0) * GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation,
rVariables.IntegrationCoefficient);
const array_1d<double, TNumNodes> coupling_flow =
PORE_PRESSURE_SIGN_FACTOR * prod(trans(p_coupling_matrix), rVariables.VelocityVector);
GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, coupling_flow);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -1493,7 +1478,7 @@ array_1d<double, TNumNodes> UPwSmallStrainElement<TDim, TNumNodes>::CalculateFlu
const BoundedMatrix<double, TNumNodes, TDim> temp_matrix =
prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient;

return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] *
return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] * rVariables.BishopCoefficient *
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know how I can verify that this change is correct... Perhaps you could try to explain it to me in-person? Thanks.

rVariables.RelativePermeability * prod(temp_matrix, rVariables.BodyAcceleration);

KRATOS_CATCH("")
Expand Down
10 changes: 0 additions & 10 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::InitializeSolutionStep(con
unsigned int NumGPoints = mConstitutiveLawVector.size();
std::vector<double> JointWidthContainer(NumGPoints);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint);
Expand All @@ -206,9 +204,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::InitializeSolutionStep(con
noalias(StressVector) = mStressVector[GPoint];
ConstitutiveParameters.SetStressVector(StressVector);
mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters);

// Initialize retention law
mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters);
markelov208 marked this conversation as resolved.
Show resolved Hide resolved
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -253,8 +248,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const
unsigned int NumGPoints = mConstitutiveLawVector.size();
std::vector<double> JointWidthContainer(NumGPoints);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint);
Expand All @@ -279,9 +272,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const

mStateVariablesFinalized[GPoint] =
mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]);

// retention law
mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);
markelov208 marked this conversation as resolved.
Show resolved Hide resolved
}

if (rCurrentProcessInfo[NODAL_SMOOTHING]) this->ExtrapolateGPValues(JointWidthContainer);
Expand Down
Loading
Loading