Skip to content

[GeoMechanicsApplication] Extend ApplyCPhiReductionProcess to use Mohr-Coulomb constitutive model#14085

Open
markelov208 wants to merge 12 commits intomasterfrom
geo/14083-extend-c-phi-reduction-with-mohr-coulomb
Open

[GeoMechanicsApplication] Extend ApplyCPhiReductionProcess to use Mohr-Coulomb constitutive model#14085
markelov208 wants to merge 12 commits intomasterfrom
geo/14083-extend-c-phi-reduction-with-mohr-coulomb

Conversation

@markelov208
Copy link
Contributor

@markelov208 markelov208 commented Dec 22, 2025

📝 Description
A brief description of the PR.

  • added usage of GEO_FRICTION_ANGLE and GEO_COHESION
  • modified checks of input data and changed the corresponding unit tests
  • modified c_phi_reduction_process.py to have two tests with implementations of Mohr-Coulomb model from MohrCoulomb64 library and internal Kratos model.
  • updated README.md file
  • added InitializeParametersForInternalMohrCoulombModel to initialize material properties for the internal Mohr-Coulomb model

Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

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

Hi Gennady,
Generally you have made the internal MC model availably for C Phi reduction. I think its good to use the names of the allowed constitutive laws for discrimination in the if, else if conditions and then do things to the parameters. The integration test files can be shrunk somewhat, by removing information that is not going to be used ( I may well be the cause of this myself as you probably started from existing files for the new tests )

I hope that Anne or Richard also will have a look.
Thank you, Wijtze Pieter

Comment on lines +101 to +110
if (r_properties.Has(GEO_COHESION)) {
check_properties.Check(GEO_COHESION);
check_properties.Check(GEO_FRICTION_ANGLE);
} else {
check_properties.CheckAvailability(UMAT_PARAMETERS);
check_properties.Check(INDEX_OF_UMAT_PHI_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
check_properties.Check(INDEX_OF_UMAT_C_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
}
Copy link
Contributor

Choose a reason for hiding this comment

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

GEO_COHESION may be there due to a request for capacity output. I think that the if and the else if condition should use the constitutive law to discriminate. So check for the MohrCoulomb ( continuum and interface law ) and the UDSM law, and then check the properties of their parameters.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thank you for the info. Now it checks for UDSM model and else for others.

p_new_properties->SetValue(GEO_COHESION, ReducedC);
} else {
// Overwrite C and Phi in the UMAT_PARAMETERS
auto Umat_parameters = r_properties[UMAT_PARAMETERS];
Copy link
Contributor

Choose a reason for hiding this comment

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

umat_parameters is a local variable starting with lower case according to the stylesheet.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

changed

const auto& r_properties = rElement.GetProperties();
Properties::Pointer p_new_properties = Kratos::make_shared<Properties>(r_properties);

if (r_properties.Has(GEO_FRICTION_ANGLE)) {
Copy link
Contributor

Choose a reason for hiding this comment

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

Also here its probably better to distinguish based on the constitutive law.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Agreed, now it is like for Check function

"SATURATED_SATURATION" : 1.0,
"K0_MAIN_DIRECTION" : 1,
"K0_VALUE_XX" : 0.42642356364895390389196808717384,
"K0_VALUE_YY" : 0.42642356364895390389196808717384,
Copy link
Contributor

Choose a reason for hiding this comment

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

The K0 main direction is 1 ( Y direction ). Then the K0_VALUE_YY is not really used, but interpreted as 1. If possible remove K0_VALUE_YY, otherwise please set it to 1.

The values for XX and ZZ are a surprise too. In the next stage the material has a friction angle of 30 deg. That would give a K0 = 1 - sin 30 = 0.5

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Now K0_VALUE_YY=1. the angle is 35 degs, not 30 degs.

Comment on lines +53 to +55
"rayleigh_k": 0.0,
"rayleigh_m": 0.0,
"reduction_factor": 0.5,
Copy link
Contributor

Choose a reason for hiding this comment

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

for a quasistatic analysis, Rayleigh parameters are not used.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

removed

"plane_output": [],
"nodal_results": ["DISPLACEMENT"]
},
"point_data_configuration": []
Copy link
Contributor

Choose a reason for hiding this comment

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

the empty point_data_configuration can probably be removed.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

removed

Comment on lines +108 to +109
"active": [true,true,true],
"is_fixed": [true,false,true],
Copy link
Contributor

Choose a reason for hiding this comment

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

Please set the Z components to false

Copy link
Contributor Author

Choose a reason for hiding this comment

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

false

"model_part_name": "PorousDomain.All_Nodes",
"out_of_plane_direction": 2,
"second_reference_coordinate": [ 40.0, -10.0, 0.0 ],
"specific_weight": 10000.0,
Copy link
Contributor

Choose a reason for hiding this comment

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

conflicting with the water density in the material parameters.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

assigned to 0

Comment on lines +120 to +121
"active": [true,true,true],
"is_fixed": [true,true,true],
Copy link
Contributor

Choose a reason for hiding this comment

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

Please set the Z components to false.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

false


KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
(ApplyCPhiReductionProcess{model, parameters}.Check()),
Copy link
Contributor

Choose a reason for hiding this comment

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

can the ExecuteInitializeSolutionStep really be omitted?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, because I moved checks from various runtime functions to Check. ExecuteInitializeSolutionStep was used here just to check the input.

Copy link
Contributor Author

@markelov208 markelov208 left a comment

Choose a reason for hiding this comment

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

Hi Wijtze Pieter, many thanks for the thorough review.


KRATOS_EXPECT_EXCEPTION_IS_THROWN(
(ApplyCPhiReductionProcess{model, parameters}.ExecuteInitializeSolutionStep()),
(ApplyCPhiReductionProcess{model, parameters}.Check()),
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, because I moved checks from various runtime functions to Check. ExecuteInitializeSolutionStep was used here just to check the input.

Comment on lines +101 to +110
if (r_properties.Has(GEO_COHESION)) {
check_properties.Check(GEO_COHESION);
check_properties.Check(GEO_FRICTION_ANGLE);
} else {
check_properties.CheckAvailability(UMAT_PARAMETERS);
check_properties.Check(INDEX_OF_UMAT_PHI_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
check_properties.Check(INDEX_OF_UMAT_C_PARAMETER, 1,
static_cast<int>(r_properties[UMAT_PARAMETERS].size()));
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thank you for the info. Now it checks for UDSM model and else for others.

p_new_properties->SetValue(GEO_COHESION, ReducedC);
} else {
// Overwrite C and Phi in the UMAT_PARAMETERS
auto Umat_parameters = r_properties[UMAT_PARAMETERS];
Copy link
Contributor Author

Choose a reason for hiding this comment

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

changed

const auto& r_properties = rElement.GetProperties();
Properties::Pointer p_new_properties = Kratos::make_shared<Properties>(r_properties);

if (r_properties.Has(GEO_FRICTION_ANGLE)) {
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Agreed, now it is like for Check function

"SATURATED_SATURATION" : 1.0,
"K0_MAIN_DIRECTION" : 1,
"K0_VALUE_XX" : 0.42642356364895390389196808717384,
"K0_VALUE_YY" : 0.42642356364895390389196808717384,
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Now K0_VALUE_YY=1. the angle is 35 degs, not 30 degs.

Comment on lines +52 to +54
"desired_iterations": 4,
"max_radius_factor": 10.0,
"min_radius_factor": 0.1,
Copy link
Contributor Author

Choose a reason for hiding this comment

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

removed

"plane_output": [],
"nodal_results": ["DISPLACEMENT"]
},
"point_data_configuration": []
Copy link
Contributor Author

Choose a reason for hiding this comment

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

removed

Comment on lines +108 to +109
"active": [true,true,true],
"is_fixed": [true,false,true],
Copy link
Contributor Author

Choose a reason for hiding this comment

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

false

"model_part_name": "PorousDomain.All_Nodes",
"out_of_plane_direction": 2,
"second_reference_coordinate": [ 40.0, -10.0, 0.0 ],
"specific_weight": 10000.0,
Copy link
Contributor Author

Choose a reason for hiding this comment

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

assigned to 0

Comment on lines +120 to +121
"active": [true,true,true],
"is_fixed": [true,true,true],
Copy link
Contributor Author

Choose a reason for hiding this comment

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

false

@markelov208 markelov208 requested a review from WPK4FEM March 10, 2026 16:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[GeoMechanicsApplication] Extend ApplyCPhiReductionProcess to use Mohr-Coulomb constitutive model

2 participants