Skip to content

[GeomechanicsApplication] Geo/mc adaptive substepping - #14572

Merged
aronnoordam merged 9 commits into
masterfrom
geo/mc_adaptive_substepping
Aug 4, 2026
Merged

[GeomechanicsApplication] Geo/mc adaptive substepping#14572
aronnoordam merged 9 commits into
masterfrom
geo/mc_adaptive_substepping

Conversation

@aronnoordam

Copy link
Copy Markdown
Member

Added substepping to Mohr coulomb law

@aronnoordam aronnoordam added the GeoMechanics Issues related to the GeoMechanicsApplication label Jul 13, 2026
@aronnoordam
aronnoordam marked this pull request as ready for review July 14, 2026 11:44
@aronnoordam
aronnoordam requested a review from a team as a code owner July 14, 2026 11:44

@markelov208 markelov208 left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Hi Aron, thank you for this substep model. I have a few comments now.

Comment thread applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_law.cpp Outdated
Comment on lines +140 to +145
if (rMaterialProperties.Has(GEO_MAX_RELATIVE_OVERSHOOT)) {
mMaxRelativeOvershoot = rMaterialProperties[GEO_MAX_RELATIVE_OVERSHOOT];
}
if (rMaterialProperties.Has(GEO_MAX_NUMBER_OF_SUB_STEPS)) {
mMaxNumberOfSubSteps = rMaterialProperties[GEO_MAX_NUMBER_OF_SUB_STEPS];
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The user may make mistakes. To report it, please add these lines to Check function

if (rMaterialProperties.Has(GEO_MAX_RELATIVE_OVERSHOOT)) {
    check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive)
        .Check(GEO_MAX_RELATIVE_OVERSHOOT);
}

if (rMaterialProperties.Has(GEO_MAX_NUMBER_OF_SUB_STEPS)) {
    check_properties.Check(GEO_MAX_NUMBER_OF_SUB_STEPS, 1);
}

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

i altered it slightly to check for the correct bounds

Comment on lines +173 to +178
if (rMaterialProperties.Has(GEO_MAX_RELATIVE_OVERSHOOT)) {
mMaxRelativeOvershoot = rMaterialProperties[GEO_MAX_RELATIVE_OVERSHOOT];
}
if (rMaterialProperties.Has(GEO_MAX_NUMBER_OF_SUB_STEPS)) {
mMaxNumberOfSubSteps = rMaterialProperties[GEO_MAX_NUMBER_OF_SUB_STEPS];
}

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

the same as for interface_coulomb file. please add these lines to Check function

if (rMaterialProperties.Has(GEO_MAX_RELATIVE_OVERSHOOT)) {
    check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive)
        .Check(GEO_MAX_RELATIVE_OVERSHOOT);
}

if (rMaterialProperties.Has(GEO_MAX_NUMBER_OF_SUB_STEPS)) {
    check_properties.Check(GEO_MAX_NUMBER_OF_SUB_STEPS, 1);
}

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

i altered it slightly to check for the correct bounds

@markelov208 markelov208 left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Hi Aron, I've added two more comments. Thank you very much for the unit tests. I think they may help to improve the current implementation of MC model. Unfortunately, I do not have time for this activity. My opinion is that it is very useful to ask @mnabideltares , @WPK4FEM and @avdg81 to check why the current MC model is not accurate for these tests.

Comment on lines +603 to +607
const auto error_no_sub_step = norm_2(Vector{traction_no_sub_stepping - reference_traction});
const auto error_with_sub_step = norm_2(Vector{traction_with_sub_stepping - reference_traction});

// Check that the error without sub-stepping is at least 10x larger than the error with sub-stepping.
KRATOS_EXPECT_GT(error_no_sub_step, error_with_sub_step * 10);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

My opinion is better to compare the obtained solutions with expected values and not the solutions. It is more visible what is going on.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

i can add the check for expected value, but i find it very relevant to check if it also much more accurate than not using sub stepping. Because i want to show the added value of sub stepping

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I think the numbers speak for themselves. When I printed values without sub-stepping, the value was very very far from the expected one.

Comment on lines +1097 to +1101
const auto error_no_sub_step = norm_2(Vector{stress_no_sub_stepping - reference_stress});
const auto error_with_sub_step = norm_2(Vector{stress_with_sub_stepping - reference_stress});

// Check that the error without sub-stepping is at least 10x larger than the error with sub-stepping
KRATOS_EXPECT_GT(error_no_sub_step, error_with_sub_step * 10);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

as for interface_coulomb_law I think it better to compare solutions with values. I printed solutions obtained for the first and second legs separately using zero_state as initial and the solutions are almost the same with and without sub-stepping. It may be that something is missed when the second leg is applied after the first leg. I would propose @mnabideltares , @WPK4FEM and @avdg81 take a look. Perhaps, something shall be improved in the current MC model. Many thanks for this test.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

hi @markelov208, the issue in the current mohr coulomb is that it doesnt recognize well when the yield surface changes, while already being on a yield surface. Because in that case the flow direction changes, but the calculation is still done with the initial flow direction. There might be a better way to do this, but substepping is a relative easy solution. Also when there is hardening involved, sub stepping helps

@avdg81 avdg81 left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Hi Aron,

First of all, thank you very much for preparing this nice extension that will help to stabilize analyses involving the Mohr-Coulomb material model. I haven't looked into the PR in much detail yet, but based on what I've seen so far, I have a few ideas that I would like to share with you, just to see how you feel about them.

I'm really curious whether we are able to successfully run the CROW case model with sand only using the sub-stepping mechanism. It's just another variant of the CROW case, but so far, we haven't been able to obtain a converged solution for it (that is, using linear iteration only, I didn't try it yet with a quasi-Newton method). Perhaps with the sub-stepping mechanism, we can. For your convenience, I have prepared two material input files for that case, and attached them below. I realize that having another set of initial materials (sand only), requires some changes in the test script to nicely fit it in. But I am very willing to help you with that.

initial_materials.json
MaterialParameters.json

The other thing that I noticed is that the sub-stepping mechanism has been implemented twice now: once for the interface constitutive law and another time for the continuum law. However, at least conceptually, it should be the same thing, I think. Of course, the considered stress states are different (traction vectors for interfaces and stress vectors for continua), but the mechanism itself should be identical. Or am I overlooking something here?
In the past, to reduce the amount of duplication, we had introduced a class for the Coulomb implementation (CoulombImpl), which is shared between the two constitutive laws. Intuitively, I would find it a natural place for the sub-stepping mechanism. Have you considered this option (e.g., we could perhaps fit it into function CoulombImpl::DoReturnMapping)? If yes, what was keeping you from putting it there? If no, may I ask you to have a look at it, to see if my idea would be feasible? Please note that I'm not asking you to implement this idea right away. Just see whether you think it can be done. Of course, if you would prefer me to explore that in more detail (or perhaps together), I/we can. It's just that I would like to give you the opportunity to do this yourself, since it's your work. Please, let me know what you prefer. Thank you.

@aronnoordam

Copy link
Copy Markdown
Member Author

Hi Aron,

First of all, thank you very much for preparing this nice extension that will help to stabilize analyses involving the Mohr-Coulomb material model. I haven't looked into the PR in much detail yet, but based on what I've seen so far, I have a few ideas that I would like to share with you, just to see how you feel about them.

I'm really curious whether we are able to successfully run the CROW case model with sand only using the sub-stepping mechanism. It's just another variant of the CROW case, but so far, we haven't been able to obtain a converged solution for it (that is, using linear iteration only, I didn't try it yet with a quasi-Newton method). Perhaps with the sub-stepping mechanism, we can. For your convenience, I have prepared two material input files for that case, and attached them below. I realize that having another set of initial materials (sand only), requires some changes in the test script to nicely fit it in. But I am very willing to help you with that.

initial_materials.json MaterialParameters.json

The other thing that I noticed is that the sub-stepping mechanism has been implemented twice now: once for the interface constitutive law and another time for the continuum law. However, at least conceptually, it should be the same thing, I think. Of course, the considered stress states are different (traction vectors for interfaces and stress vectors for continua), but the mechanism itself should be identical. Or am I overlooking something here? In the past, to reduce the amount of duplication, we had introduced a class for the Coulomb implementation (CoulombImpl), which is shared between the two constitutive laws. Intuitively, I would find it a natural place for the sub-stepping mechanism. Have you considered this option (e.g., we could perhaps fit it into function CoulombImpl::DoReturnMapping)? If yes, what was keeping you from putting it there? If no, may I ask you to have a look at it, to see if my idea would be feasible? Please note that I'm not asking you to implement this idea right away. Just see whether you think it can be done. Of course, if you would prefer me to explore that in more detail (or perhaps together), I/we can. It's just that I would like to give you the opportunity to do this yourself, since it's your work. Please, let me know what you prefer. Thank you.

Hi Anne,

the substepping is indeed the same for both cases except the mpConstitutiveDimension and stress states. But indeed perhaps it can be added in CoulombWithTensionCutOffImpl.

But, i think the substepping method as implemented is 1 method but not the only method, more methods are available, and are not limited to be used in Mohr coulomb, but also in other constitutive models. Maybe you have an idea, how to make this even more general, than just adding it to CoulombWithTensionCutOffImpl

@avdg81 avdg81 left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Hi Aron,
Thanks for taking care of the review suggestions. This PR is good to go.

@github-project-automation github-project-automation Bot moved this from In progress to Reviewer approved in GeoMechanicsApplication (Deltares) Aug 4, 2026
@aronnoordam
aronnoordam merged commit aadb866 into master Aug 4, 2026
10 checks passed
@aronnoordam
aronnoordam deleted the geo/mc_adaptive_substepping branch August 4, 2026 14:07
@github-project-automation github-project-automation Bot moved this from Reviewer approved to Done in GeoMechanicsApplication (Deltares) Aug 4, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

GeoMechanics Issues related to the GeoMechanicsApplication

Development

Successfully merging this pull request may close these issues.

[GeoMechanicsApplication] Add sub-stepping to the (Mohr/interface) Coulomb law

3 participants