-
Notifications
You must be signed in to change notification settings - Fork 248
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.
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
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi Mohamed and Wijtze-Pieter, the code looks very professional! Sorry I missed what was done to fix the problem. Based on the recent discussion, a higher order shall be used; however, I do not see it here. Does it mean that this PR's goal is to make the code professional? Many thanks.
applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
Show resolved
Hide resolved
applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
Outdated
Show resolved
Hide resolved
...ly_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage2.json
Outdated
Show resolved
Hide resolved
...y_saturated/test_saturated_below_phreatic_level_pw_triangle6N/MaterialParameters_stage1.json
Outdated
Show resolved
Hide resolved
...y_saturated/test_saturated_below_phreatic_level_pw_triangle6N/MaterialParameters_stage2.json
Outdated
Show resolved
Hide resolved
...ly_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage1.json
Outdated
Show resolved
Hide resolved
...ly_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage2.json
Outdated
Show resolved
Hide resolved
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Gentlemen, you two have done a lot of work to clean up the code and to fix a few issues. I'm very satisfied with the result. It's yet another step towards code that is easier to understand and more maintainable.
I do have quite a few suggestions and some questions, but none of them is blocking in my opinion.
else if (rVariable == RELATIVE_PERMEABILITY) | ||
rOutput[GPoint] = | ||
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); | ||
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue(RetentionParameters, rVariable, rOutput[GPoint]); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have noticed that the switch that has been taken out here, is repeated for each retention law (in their overrides of CalculateValue
). We shouldn't try to address that in this PR, but I think it would be nice to pick that up soon after merging this one. I would propose to implement the switch in class RetentionLaw
and to not make CalculateValue
virtual. In that way, the switch is implemented only once, in a place where you would expect it to be.
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); |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Perhaps we can try to apply the Rule of zero here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Tried with some help.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Perhaps we can try to apply the Rule of zero here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Tried.
SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw() : RetentionLaw() {} | ||
|
||
SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther) | ||
: RetentionLaw(rOther) | ||
{ | ||
} | ||
SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther) = default; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We should be able to remove these when we apply the Rule of zero.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
tried.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
tried.
Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT | ||
Begin SubModelPartNodes | ||
1 | ||
2 | ||
43 | ||
End SubModelPartNodes | ||
End SubModelPart |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is an unused sub-model part.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed
End SubModelPartNodes | ||
End SubModelPart | ||
|
||
Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
And also this sub-model part is not being used.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed
End SubModelPartNodes | ||
End SubModelPart | ||
|
||
Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also this sub-model part is not being used...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am missing displacement boundary conditions here. I had expected to see that the bottom edge would be fixed, and sliders along the sides of the column. Shouldn't we better add these?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed the bottom boundary, as it is not necessary
End SubModelPartNodes | ||
End SubModelPart | ||
|
||
Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This sub-model part should be referenced from the ProjectParameters.json
files to fix the displacements at the bottom, but that is currently not being done.
UPDATE: Now I see that you fix all displacement DOFs, by prescribing zeroes for them. In that case, we don't need this sub-model part.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Removed
…est for the errors.
6ace9cb
to
a8cd4ca
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for all the hard work. I believe this PR is ready to be merged.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let me approve it once more :-)
This issue is time-boxed to 10 developer days.
22/10/2024 - JDN - Comment
I need more information to able to prioritize this for this sprint or beyond. i.e. whats the urgency ? who is it for ? also needs refinement ? I wasn't aware we were using partial saturation anywhere and we are currently awaiting for 2025 to do so, but do I interpret correctly that the mechanism also fails when the zone above the phreatic line is fully unsaturated
Description
The Partially saturated flow is implemented but it has been never tested. Moreover, instability arrised when we tried to apply it.
The equation for partially saturated flow is:
where
Here we start from a simple test case, with saturated flow below the phreatic line:$S_r = 0$ , the pressure above the phreatic line needs to be zero. However, from this equation, the left hand side includes the saturation $S$ which is zero. It leads to a zero $C$ matrix. Morover, $k_r = 0$ too, which leads to a zero $K$ matrix. Therefore, a zero matrix is added which makes the diagonal of the generic matrix partly to be zero, and it lead to a singular matrix.$C$ matrix manually as a user defined parameter.
In the case of
This can be probably avoided to adding "BIOT_COEFFICIENT" to the material parameters to set the
In the case of$S_r > 0$ , then $K_r =0$ `but the left hand side gets a nonzero value. it means $\frac{dp}{dt} = 0$ . This is proposed to lead to unchanged pressure above the phreatic line.
However, here we see a change in pressure. It indicates there is driving force, which leads the water pressure goes up. Which is an unexpected behaviour,