Skip to content

[GeomechanicsApplication] Geo/mc adaptive substepping #14572

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

Merged
aronnoordam merged 9 commits into from
Aug 4, 2026
Merged

[GeomechanicsApplication] Geo/mc adaptive substepping #14572

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0ed5237
added mc substepping and test
aronnoordam Jul 10, 2026
5a4068e
cleanup
aronnoordam Jul 13, 2026
1fa6b11
calculates number of substeps only when required
aronnoordam Jul 13, 2026
28886dc
cleanup
aronnoordam Jul 13, 2026
bcc41e2
added substepping to interface coulomb
aronnoordam Jul 14, 2026
c549b6c
cleanup
aronnoordam Jul 14, 2026
ad56f65
applied comments
aronnoordam Jul 24, 2026
108aa99
formatted clang format
aronnoordam Aug 3, 2026
75105e0
formatted clang format
aronnoordam Aug 3, 2026
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
115 changes: 102 additions & 13 deletions applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_law.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,9 @@ ConstitutiveLaw::Pointer InterfaceCoulombLaw::Clone() const
p_result->mRelativeDisplacementVectorFinalized = mRelativeDisplacementVectorFinalized;
p_result->mpCoulombImpl = mpCoulombImpl->Clone();
p_result->mIsModelInitialized = mIsModelInitialized;
p_result->mMaxRelativeOvershoot = mMaxRelativeOvershoot;
p_result->mMaxNumberOfSubSteps = mMaxNumberOfSubSteps;
p_result->mCalculatedNumberOfSubSteps = mCalculatedNumberOfSubSteps;
return p_result;
}

Expand Down Expand Up @@ -102,6 +105,16 @@ int InterfaceCoulombLaw::Check(const Properties& rMaterialProperties,
}
check_properties.Check(INTERFACE_NORMAL_STIFFNESS);
check_properties.Check(INTERFACE_SHEAR_STIFFNESS);

if (rMaterialProperties.Has(GEO_MAX_RELATIVE_OVERSHOOT)) {
check_properties.SingleUseBounds(CheckProperties::Bounds::ExclusiveLowerAndInclusiveUpper)
.Check(GEO_MAX_RELATIVE_OVERSHOOT, 0.0, 1.0);
}

if (rMaterialProperties.Has(GEO_MAX_NUMBER_OF_SUB_STEPS)) {
constexpr auto max_value_number_of_substeps = std::numeric_limits<int>::max();
check_properties.Check(GEO_MAX_NUMBER_OF_SUB_STEPS, 1, max_value_number_of_substeps);
}
return result;
}

Expand Down Expand Up @@ -133,6 +146,13 @@ void InterfaceCoulombLaw::InitializeMaterial(const Properties& rMaterialProperti
HasInitialState() ? GetInitialState().GetInitialStrainVector() : ZeroVector{GetStrainSize()};
mTractionVectorFinalized =
HasInitialState() ? GetInitialState().GetInitialStressVector() : ZeroVector{GetStrainSize()};

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];
}
Comment on lines +150 to +155

@markelov208 markelov208 Jul 15, 2026

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);
}

@aronnoordam aronnoordam Jul 24, 2026

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

}

void InterfaceCoulombLaw::InitializeMaterialResponseCauchy(Parameters& rConstitutiveLawParameters)
Expand All @@ -151,26 +171,90 @@ void InterfaceCoulombLaw::CalculateMaterialResponseCauchy(Parameters& rConstitut
}

const auto& r_properties = rConstitutiveLawParameters.GetMaterialProperties();
const auto elastic_matrix = mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties);
const auto normal_stiffness = r_properties[INTERFACE_NORMAL_STIFFNESS];
const auto shear_stiffness = r_properties[INTERFACE_SHEAR_STIFFNESS];

// Elastic predictor over the full increment starting from the last finalized state.
const auto full_trial_sigma_tau = CalculateTrialTractionVector(
rConstitutiveLawParameters.GetStrainVector(), normal_stiffness, shear_stiffness);

// If the whole step stays elastic, there is no need to sub-step.
auto aux_sigma_tau = full_trial_sigma_tau;
aux_sigma_tau.Tau() = std::abs(aux_sigma_tau.Tau());
if (mpCoulombImpl->IsAdmissibleStressState(aux_sigma_tau)) {
mTractionVector = full_trial_sigma_tau.CopyTo<Vector>();
rConstitutiveLawParameters.GetStressVector() = mTractionVector;
return;
}

std::size_t number_of_sub_steps = 1;
if (mMaxNumberOfSubSteps > 1) {
// only calculate the number of sub-steps once per solution step for stability
if (mCalculatedNumberOfSubSteps == 0) {
mCalculatedNumberOfSubSteps = CalculateAdaptiveNumberOfSubSteps(aux_sigma_tau, elastic_matrix);
}
number_of_sub_steps = mCalculatedNumberOfSubSteps;
}

auto trial_sigma_tau = CalculateTrialTractionVector(rConstitutiveLawParameters.GetStrainVector(),
r_properties[INTERFACE_NORMAL_STIFFNESS],
r_properties[INTERFACE_SHEAR_STIFFNESS]);
auto mapped_sigma_tau = trial_sigma_tau;
// Running committed state for the sub-stepping (start from last finalized state).
Vector committed_traction = mTractionVectorFinalized;
Vector committed_displacement = mRelativeDisplacementVectorFinalized;

const auto negative = std::signbit(trial_sigma_tau.Tau());
trial_sigma_tau.Tau() = std::abs(trial_sigma_tau.Tau());
const Vector total_displacement_increment =
rConstitutiveLawParameters.GetStrainVector() - mRelativeDisplacementVectorFinalized;

if (!mpCoulombImpl->IsAdmissibleStressState(trial_sigma_tau)) {
mapped_sigma_tau = mpCoulombImpl->DoReturnMapping(
trial_sigma_tau, mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties),
Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
if (negative) mapped_sigma_tau.Tau() *= -1.0;
constexpr auto number_of_normal_components = std::size_t{1};
const auto elastic_tensor = ConstitutiveLawUtilities::MakeInterfaceElasticConstitutiveTensor(
normal_stiffness, shear_stiffness, GetStrainSize(), number_of_normal_components);

for (std::size_t sub = 1; sub <= number_of_sub_steps; ++sub) {
const Vector sub_displacement =
mRelativeDisplacementVectorFinalized +
(static_cast<double>(sub) / static_cast<double>(number_of_sub_steps)) * total_displacement_increment;

auto trial_sigma_tau = Geo::SigmaTau{
committed_traction + prod(elastic_tensor, sub_displacement - committed_displacement)};
auto mapped_sigma_tau = trial_sigma_tau;

const auto negative = std::signbit(trial_sigma_tau.Tau());
trial_sigma_tau.Tau() = std::abs(trial_sigma_tau.Tau());

if (!mpCoulombImpl->IsAdmissibleStressState(trial_sigma_tau)) {
mapped_sigma_tau = mpCoulombImpl->DoReturnMapping(
trial_sigma_tau, elastic_matrix, Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
if (negative) mapped_sigma_tau.Tau() *= -1.0;
}

mTractionVector = mapped_sigma_tau.CopyTo<Vector>();
committed_traction = mTractionVector;
committed_displacement = sub_displacement;
}

mTractionVector = mapped_sigma_tau.CopyTo<Vector>();
rConstitutiveLawParameters.GetStressVector() = mTractionVector;
}

std::size_t InterfaceCoulombLaw::CalculateAdaptiveNumberOfSubSteps(const Geo::SigmaTau& rTrialTraction,
const Matrix& rElasticMatrix) const
{
// make sure that kappa is not updated while calculating the number of required sub steps
mpCoulombImpl->SaveKappaOfCoulombYieldSurface();
const auto mapped_sigma_tau = mpCoulombImpl->DoReturnMapping(
rTrialTraction, rElasticMatrix, Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
mpCoulombImpl->RestoreKappaOfCoulombYieldSurface();

const Vector trial_values = rTrialTraction.CopyTo<Vector>();
const Vector mapped_values = mapped_sigma_tau.CopyTo<Vector>();

const auto overshoot = norm_2(trial_values - mapped_values);
const auto stress_scale = std::max(norm_2(trial_values), 1.0e-12);
const auto relative_overshoot = overshoot / stress_scale;

const auto number_of_sub_steps = static_cast<int>(std::ceil(relative_overshoot / mMaxRelativeOvershoot));

return std::clamp(number_of_sub_steps, 1, mMaxNumberOfSubSteps);
}

Geo::SigmaTau InterfaceCoulombLaw::CalculateTrialTractionVector(const Vector& rRelativeDisplacementVector,
double NormalStiffness,
double ShearStiffness) const
Expand All @@ -186,6 +270,7 @@ void InterfaceCoulombLaw::FinalizeMaterialResponseCauchy(Parameters& rConstituti
{
mRelativeDisplacementVectorFinalized = rConstitutiveLawParameters.GetStrainVector();
mTractionVectorFinalized = mTractionVector;
mCalculatedNumberOfSubSteps = 0;
}

Matrix& InterfaceCoulombLaw::CalculateValue(Parameters& rConstitutiveLawParameters,
Expand Down Expand Up @@ -214,6 +299,8 @@ void InterfaceCoulombLaw::save(Serializer& rSerializer) const
rSerializer.save("RelativeDisplacementVectorFinalized", mRelativeDisplacementVectorFinalized);
rSerializer.save("CoulombImpl", mpCoulombImpl);
rSerializer.save("IsModelInitialized", mIsModelInitialized);
rSerializer.save("MaxRelativeOvershoot", mMaxRelativeOvershoot);
rSerializer.save("MaxNumberOfSubSteps", mMaxNumberOfSubSteps);
}

void InterfaceCoulombLaw::load(Serializer& rSerializer)
Expand All @@ -225,6 +312,8 @@ void InterfaceCoulombLaw::load(Serializer& rSerializer)
rSerializer.load("RelativeDisplacementVectorFinalized", mRelativeDisplacementVectorFinalized);
rSerializer.load("CoulombImpl", mpCoulombImpl);
rSerializer.load("IsModelInitialized", mIsModelInitialized);
rSerializer.load("MaxRelativeOvershoot", mMaxRelativeOvershoot);
rSerializer.load("MaxNumberOfSubSteps", mMaxNumberOfSubSteps);
}

} // Namespace Kratos
} // Namespace Kratos
6 changes: 6 additions & 0 deletions applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_law.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceCoulombLaw : public Constit
[[nodiscard]] Geo::SigmaTau CalculateTrialTractionVector(const Vector& rRelativeDisplacementVector,
double NormalStiffness,
double ShearStiffness) const;
std::size_t CalculateAdaptiveNumberOfSubSteps(const Geo::SigmaTau& rTrialTraction,
const Matrix& rElasticMatrix) const;

double mMaxRelativeOvershoot = 1.0;
int mMaxNumberOfSubSteps = 1;
std::size_t mCalculatedNumberOfSubSteps = 0;

friend class Serializer;
void save(Serializer& rSerializer) const override;
Expand Down
Loading
Loading