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

Open
wants to merge 39 commits into
base: master
Choose a base branch
from
+4,452 −939
Open

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

Show file tree
Hide file tree
Changes from 30 commits
Commits
Show all changes
39 commits
Select commit Hold shift + click to select a range
93e1263
Avoid zero divisions + a lot of printing
WPK4FEM Oct 15, 2024
6fcc7fa
Removed an unused variable and some commented out code
avdg81 Oct 25, 2024
37758a7
Fix effective saturation. Commented print statements.
WPK4FEM Dec 6, 2024
409851f
start of unit tests for SaturatedBelowPhreaticLevelLaw
WPK4FEM Dec 6, 2024
bf5bd8d
Improved error messages for saturated below phreatic level and unit t…
WPK4FEM Dec 9, 2024
188099f
removed unwanted include
WPK4FEM Dec 9, 2024
7e10e65
ternary operator
WPK4FEM Dec 9, 2024
85cb890
Unit test for saturated law
WPK4FEM Dec 9, 2024
d31ed60
ternary operator
WPK4FEM Dec 9, 2024
c4412d9
Checks and unit tests for Van Genuchten added.
WPK4FEM Dec 9, 2024
60eb250
Fixes for partially saturated flow
mnabideltares Dec 9, 2024
e203f4c
Attempt to only place includes there where they are used.
WPK4FEM Dec 9, 2024
124346e
Cleanup attempts
WPK4FEM Dec 9, 2024
62f388d
Trying to get rid of empty functions.
WPK4FEM Dec 9, 2024
28ff169
Remove added print statements from core
WPK4FEM Dec 9, 2024
6574b78
sonar cloud remarks and naming style corrections
WPK4FEM Dec 10, 2024
c5871a2
Coupling matrix computation without potential zero division also in s…
WPK4FEM Dec 10, 2024
ef05087
Also avoid zero division in RHS setup of coupling terms.
WPK4FEM Dec 10, 2024
2a98e5c
van Genuchten gl >= 0, clean unused parameters from tests, some geome…
WPK4FEM Dec 11, 2024
83debf8
IgnoreUndrained setup of LHS and RHS made symmetric.
WPK4FEM Dec 11, 2024
04ce103
No water solution intended in these tests, made IGNORE_UNDRAINED cons…
WPK4FEM Dec 12, 2024
dd719eb
Added test cases for Pw 2D3N and 2D6N, with documentation
mnabideltares Dec 12, 2024
021cc95
Merge branch 'geo/partial_saturation2' of https://github.com/KratosMu…
mnabideltares Dec 12, 2024
0a1cdba
fix for Readme
mnabideltares Dec 12, 2024
03c1ba7
A test xase for U-Pw diff order 2D6N is added
mnabideltares Dec 12, 2024
fb3b09e
fix of revision 60eb250069fc89bd84a03ad2f40f113fcfcc51dc changed resu…
WPK4FEM Dec 12, 2024
eb14f80
2 tests are added for U-Pw small strain element
mnabideltares Dec 12, 2024
5825e04
Merge branch 'geo/partial_saturation2' of https://github.com/KratosMu…
mnabideltares Dec 12, 2024
437df1f
fixes
mnabideltares Dec 12, 2024
bd43cc4
fix
mnabideltares Dec 12, 2024
162c3f4
Sonar Cloud unused stuff complaint resolved
WPK4FEM Dec 16, 2024
f7c0261
Removed unused retention parameters ( review Gennady )
WPK4FEM Dec 16, 2024
be608a3
std::transform for pressure getting ( review comment of Gennady )
WPK4FEM Dec 16, 2024
83d65fe
Removed unused K0 things ( review comment Gennady )
WPK4FEM Dec 16, 2024
35b9872
Renamed kolom to the english column ( review comment Gennady )
WPK4FEM Dec 16, 2024
27b69bc
avoid code repetition
WPK4FEM Dec 16, 2024
856256e
kp format on request of Anne
WPK4FEM Dec 16, 2024
5e78ddd
Set material parameters in a common directory, to avoid file duplicat…
mnabideltares Dec 16, 2024
eec6c84
Merge remote-tracking branch 'refs/remotes/origin/geo/partial_saturat…
WPK4FEM Dec 16, 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(
Copy link
Contributor

Choose a reason for hiding this comment

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

As far as I understand, InitializeMaterial and other retention law functions were doing nothing after certain changes. Nice to remove them finally.

Copy link
Contributor

Choose a reason for hiding this comment

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

Thank you.

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
57 changes: 28 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 @@ -283,8 +283,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);
Copy link
Contributor

Choose a reason for hiding this comment

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

RetentionParameters is not used any longer.

Copy link
Contributor

Choose a reason for hiding this comment

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

Removed

}

// Reset hydraulic discharge
Expand All @@ -299,9 +297,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 +354,7 @@ 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 @@ -434,8 +430,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 @@ -1246,17 +1240,18 @@ 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 +1275,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);

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 +1385,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 +1492,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 *
rVariables.RelativePermeability * prod(temp_matrix, rVariables.BodyAcceleration);

KRATOS_CATCH("")
Expand Down
6 changes: 0 additions & 6 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 @@ -206,9 +206,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);
Copy link
Contributor

Choose a reason for hiding this comment

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

RetentionParameters should be removed as well.

Copy link
Contributor

Choose a reason for hiding this comment

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

Removed.

}

KRATOS_CATCH("")
Expand Down Expand Up @@ -279,9 +276,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const

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

// retention law
mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);
Copy link
Contributor

Choose a reason for hiding this comment

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

Sorry, again about RetentionParameters

Copy link
Contributor

Choose a reason for hiding this comment

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

Removed.

}

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