Skip to content

[GeoMechanicsApplication] Include variable prediction in all newmark schemes #13449

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

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

[GeoMechanicsApplication] Include variable prediction in all newmark schemes #13449

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f795b45
Moved the predict from dynamic to the generalized newmark scheme
rfaasse May 19, 2025
b890fdc
Added prediction of first order scalar variable when dt_variable is n…
rfaasse May 20, 2025
2ccf9eb
Added prediction of first order scalar variable when dt_variable is f…
rfaasse May 20, 2025
13e86d3
Fixed unit test by fixing logic (if -> else if)
rfaasse May 20, 2025
6ccec00
Minor refactoring
rfaasse May 20, 2025
78b5a69
Changed predict of quasistatic UPw scheme, such that it does not pred…
rfaasse Jun 3, 2025
2855780
Merge remote-tracking branch 'origin/master' into geo/11957-include-v…
rfaasse Jun 4, 2025
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
6 changes: 4 additions & 2 deletions applications/GeoMechanicsApplication/custom_constitutive/small_strain_udsm_3D_law.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,14 +33,16 @@ constexpr auto props_size = SizeType{50};

array_1d<double, props_size> MakePropsVector(const Vector& rUMatParameters)
{
KRATOS_DEBUG_ERROR_IF(rUMatParameters.size() > props_size) << "Number of UMAT_PARAMETERS (" << rUMatParameters.size() << ") exceeds the maximum number of " << props_size << "\n";
KRATOS_DEBUG_ERROR_IF(rUMatParameters.size() > props_size)
<< "Number of UMAT_PARAMETERS (" << rUMatParameters.size()
<< ") exceeds the maximum number of " << props_size << "\n";

auto result = array_1d<double, props_size>{props_size, 0.0};
std::copy(rUMatParameters.begin(), rUMatParameters.end(), result.begin());
return result;
}

}
} // namespace

namespace Kratos
{
Expand Down
93 changes: 93 additions & 0 deletions ...ications/GeoMechanicsApplication/custom_strategies/schemes/generalized_newmark_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
using BaseType = Scheme<TSparseSpace, TDenseSpace>;
using TSystemMatrixType = typename BaseType::TSystemMatrixType;
using TSystemVectorType = typename BaseType::TSystemVectorType;
using DofsArrayType = typename BaseType::DofsArrayType;

GeneralizedNewmarkScheme(const std::vector<FirstOrderScalarVariable>& rFirstOrderScalarVariables,
const std::vector<SecondOrderVectorVariable>& rSecondOrderVectorVariables,
Expand Down Expand Up @@ -105,6 +106,17 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
KRATOS_CATCH("")
}

void Predict(ModelPart& rModelPart, DofsArrayType& rDofSet, TSystemMatrixType& A, TSystemVectorType& Dx, TSystemVectorType& b) override
{
KRATOS_TRY

PredictVariables(rModelPart);
// Update (Angular) Acceleration, (Angular) Velocity and DtPressure
this->UpdateVariablesDerivatives(rModelPart);

KRATOS_CATCH("")
}

protected:
inline void SetTimeFactors(ModelPart& rModelPart) override
{
Expand Down Expand Up @@ -181,6 +193,87 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
}

private:
void PredictVariables(const ModelPart& rModelPart)
{
block_for_each(rModelPart.Nodes(), [this](Node& rNode) { PredictVariablesForNode(rNode); });
}

void PredictVariablesForNode(Node& rNode)
{
for (const auto& r_first_order_scalar_variable : this->GetFirstOrderScalarVariables()) {
if (!rNode.SolutionStepsDataHas(r_first_order_scalar_variable.instance)) continue;
PredictFirstOrderScalarVariableForNode(rNode, r_first_order_scalar_variable);
}

for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
PredictSecondOrderVectorVariableForNode(rNode, r_second_order_vector_variable);
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

@WPK4FEM should I move this part (for the SecondOrderVectorVariables) back to the dynamic version of the newmark scheme, since only there we have non-zero accelerations/velocities?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, at least the part for acceleration.
The newmark_quasistatic_damped_U_Pw_scheme has no acceleration, but has velocities and a damping matrix, so would have a Predict without accelerations, but with prescribed velocities.

}

void PredictFirstOrderScalarVariableForNode(Node& rNode, const FirstOrderScalarVariable& rFirstOrderScalarVariable)
{
const double previous_variable =
rNode.FastGetSolutionStepValue(rFirstOrderScalarVariable.instance, 1);
const double previous_first_time_derivative =
rNode.FastGetSolutionStepValue(rFirstOrderScalarVariable.first_time_derivative, 1);
const double current_first_time_derivative =
rNode.FastGetSolutionStepValue(rFirstOrderScalarVariable.first_time_derivative, 0);

if (rNode.IsFixed(rFirstOrderScalarVariable.first_time_derivative)) {
rNode.FastGetSolutionStepValue(rFirstOrderScalarVariable.instance) =
previous_variable + (1 - this->GetTheta()) * this->GetDeltaTime() * previous_first_time_derivative +
current_first_time_derivative * this->GetTheta() * this->GetDeltaTime();
} else if (!rNode.IsFixed(rFirstOrderScalarVariable.instance)) {
rNode.FastGetSolutionStepValue(rFirstOrderScalarVariable.instance) =
previous_variable + this->GetDeltaTime() * previous_first_time_derivative;
}
}

void PredictSecondOrderVectorVariableForNode(Node& rNode, const SecondOrderVectorVariable& rSecondOrderVectorVariable)
{
const std::vector<std::string> components = {"X", "Y", "Z"};

for (const auto& component : components) {
const auto& instance_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVectorVariable.instance.Name(), component);

if (!rNode.HasDofFor(instance_component)) continue;

const auto& first_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVectorVariable.first_time_derivative.Name(), component);
const auto& second_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVectorVariable.second_time_derivative.Name(), component);

const double previous_variable = rNode.FastGetSolutionStepValue(instance_component, 1);
const double current_first_time_derivative =
rNode.FastGetSolutionStepValue(first_time_derivative_component, 0);
const double previous_first_time_derivative =
rNode.FastGetSolutionStepValue(first_time_derivative_component, 1);
const double current_second_time_derivative =
rNode.FastGetSolutionStepValue(second_time_derivative_component, 0);
const double previous_second_time_derivative =
rNode.FastGetSolutionStepValue(second_time_derivative_component, 1);
if (rNode.IsFixed(second_time_derivative_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable + this->GetDeltaTime() * previous_first_time_derivative +
this->GetDeltaTime() * this->GetDeltaTime() *
((0.5 - this->GetBeta()) * previous_second_time_derivative +
this->GetBeta() * current_second_time_derivative);
} else if (rNode.IsFixed(first_time_derivative_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable +
this->GetDeltaTime() * ((this->GetBeta() / this->GetGamma()) *
(current_first_time_derivative - previous_first_time_derivative) +
previous_first_time_derivative);
} else if (!rNode.IsFixed(instance_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable + this->GetDeltaTime() * previous_first_time_derivative +
0.5 * this->GetDeltaTime() * this->GetDeltaTime() * previous_second_time_derivative;
}
}
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Most of this just moved from the dynamic scheme, the only new part is the PredictFirstOrderScalarVariableForNode function


std::optional<double> mBeta;
std::optional<double> mGamma;
std::optional<double> mTheta;
Expand Down
69 changes: 0 additions & 69 deletions ...cations/GeoMechanicsApplication/custom_strategies/schemes/newmark_dynamic_U_Pw_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@ class NewmarkDynamicUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace, TD
KRATOS_CLASS_POINTER_DEFINITION(NewmarkDynamicUPwScheme);

using BaseType = Scheme<TSparseSpace, TDenseSpace>;
using DofsArrayType = typename BaseType::DofsArrayType;
using TSystemMatrixType = typename BaseType::TSystemMatrixType;
using TSystemVectorType = typename BaseType::TSystemVectorType;
using LocalSystemVectorType = typename BaseType::LocalSystemVectorType;
Expand All @@ -52,74 +51,6 @@ class NewmarkDynamicUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace, TD
mVelocityVector.resize(num_threads);
}

void Predict(ModelPart& rModelPart, DofsArrayType& rDofSet, TSystemMatrixType& A, TSystemVectorType& Dx, TSystemVectorType& b) override
{
KRATOS_TRY

PredictVariables(rModelPart);
// Update (Angular) Acceleration, (Angular) Velocity and DtPressure
this->UpdateVariablesDerivatives(rModelPart);

KRATOS_CATCH("")
}

void PredictVariables(const ModelPart& rModelPart)
{
block_for_each(rModelPart.Nodes(), [this](Node& rNode) { PredictVariablesForNode(rNode); });
}

void PredictVariablesForNode(Node& rNode)
{
for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
PredictVariableForNode(rNode, r_second_order_vector_variable);
}
}

void PredictVariableForNode(Node& rNode, const SecondOrderVectorVariable& rSecondOrderVariables)
{
const std::vector<std::string> components = {"X", "Y", "Z"};

for (const auto& component : components) {
const auto& instance_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVariables.instance.Name(), component);

if (!rNode.HasDofFor(instance_component)) continue;

const auto& first_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVariables.first_time_derivative.Name(), component);
const auto& second_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
rSecondOrderVariables.second_time_derivative.Name(), component);

const double previous_variable = rNode.FastGetSolutionStepValue(instance_component, 1);
const double current_first_time_derivative =
rNode.FastGetSolutionStepValue(first_time_derivative_component, 0);
const double previous_first_time_derivative =
rNode.FastGetSolutionStepValue(first_time_derivative_component, 1);
const double current_second_time_derivative =
rNode.FastGetSolutionStepValue(second_time_derivative_component, 0);
const double previous_second_time_derivative =
rNode.FastGetSolutionStepValue(second_time_derivative_component, 1);
if (rNode.IsFixed(second_time_derivative_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable + this->GetDeltaTime() * previous_first_time_derivative +
this->GetDeltaTime() * this->GetDeltaTime() *
((0.5 - this->GetBeta()) * previous_second_time_derivative +
this->GetBeta() * current_second_time_derivative);
} else if (rNode.IsFixed(first_time_derivative_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable +
this->GetDeltaTime() * ((this->GetBeta() / this->GetGamma()) *
(current_first_time_derivative - previous_first_time_derivative) +
previous_first_time_derivative);
} else if (!rNode.IsFixed(instance_component)) {
rNode.FastGetSolutionStepValue(instance_component) =
previous_variable + this->GetDeltaTime() * previous_first_time_derivative +
0.5 * this->GetDeltaTime() * this->GetDeltaTime() * previous_second_time_derivative;
}
}
}

void CalculateSystemContributions(Condition& rCurrentCondition,
LocalSystemMatrixType& rLHS_Contribution,
LocalSystemVectorType& rRHS_Contribution,
Expand Down
2 changes: 1 addition & 1 deletion ...hanicsApplication/tests/cpp_tests/custom_strategies/test_generalized_newmark_T_scheme.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ KRATOS_TEST_CASE_IN_SUITE(NewmarkTSchemeUpdate_SetsDtTemperature, KratosGeoMecha
scheme.Predict(model_part, dof_set, A, Dx, b);

// This is the expected value as calculated by the UpdateVariablesDerivatives
KRATOS_EXPECT_DOUBLE_EQ(node.FastGetSolutionStepValue(DT_TEMPERATURE, 0), 7.0 / 3.0);
KRATOS_EXPECT_DOUBLE_EQ(node.FastGetSolutionStepValue(DT_TEMPERATURE, 0), 3.0);
}

KRATOS_TEST_CASE_IN_SUITE(InitializeNewmarkTScheme_SetsTimeFactors, KratosGeoMechanicsFastSuiteWithoutKernel)
Expand Down
2 changes: 1 addition & 1 deletion ...ions/GeoMechanicsApplication/tests/cpp_tests/custom_strategies/test_newmark_Pw_scheme.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ KRATOS_TEST_CASE_IN_SUITE(NewmarkPwSchemeUpdate_SetsDtPressure, KratosGeoMechani
scheme.Predict(model_part, dof_set, A, Dx, b);

// This is the expected value as calculated by the UpdateVariablesDerivatives
KRATOS_EXPECT_DOUBLE_EQ(node.FastGetSolutionStepValue(DT_WATER_PRESSURE, 0), 7.0 / 3.0);
KRATOS_EXPECT_DOUBLE_EQ(node.FastGetSolutionStepValue(DT_WATER_PRESSURE, 0), 3.0);
}

KRATOS_TEST_CASE_IN_SUITE(InitializeNewmarkPwScheme_SetsTimeFactors, KratosGeoMechanicsFastSuiteWithoutKernel)
Expand Down
3 changes: 2 additions & 1 deletion ...chanicsApplication/tests/cpp_tests/custom_strategies/test_newmark_dynamic_U_Pw_scheme.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,14 +204,15 @@ KRATOS_TEST_CASE_IN_SUITE(NewmarkDynamicUPwSchemePredict_UpdatesVariablesDerivat
Vector Dx;
Vector b;

tester.GetModelPart().GetNode(0).FastGetSolutionStepValue(DT_WATER_PRESSURE, 1) = 3.0;
tester.mScheme.InitializeSolutionStep(tester.GetModelPart(), A, Dx, b); // This is needed to set the time factors

tester.mScheme.Predict(tester.GetModelPart(), dof_set, A, Dx, b);

// These expected numbers result from the calculations in UpdateVariablesDerivatives
const auto expected_acceleration = Kratos::array_1d<double, 3>{13, 14, 15};
const auto expected_velocity = Kratos::array_1d<double, 3>{62, 67, 72};
constexpr auto expected_dt_water_pressure = 1.0 / 3.0;
constexpr auto expected_dt_water_pressure = 3.0;

const auto actual_acceleration =
tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(ACCELERATION, 0);
Expand Down
44 changes: 43 additions & 1 deletion ...icsApplication/tests/cpp_tests/custom_strategies/test_newmark_quasistatic_U_Pw_scheme.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,47 @@ KRATOS_TEST_CASE_IN_SUITE(InitializeNewmarkUPwScheme_SetsTimeFactors, KratosGeoM
expected_velocity_coefficient);
}

KRATOS_TEST_CASE_IN_SUITE(NewmarkUPwSchemePredict_PredictsWaterPressure, KratosGeoMechanicsFastSuiteWithoutKernel)
{
NewmarkQuasistaticUPwSchemeTester tester;

ModelPart::DofsArrayType dof_set;
CompressedMatrix A;
Vector Dx;
Vector b;
tester.GetModelPart().GetProcessInfo()[DELTA_TIME] = 4.0;
tester.GetModelPart().GetNode(0).FastGetSolutionStepValue(DT_WATER_PRESSURE, 1) = 3.0;
tester.mScheme.InitializeSolutionStep(tester.GetModelPart(), A, Dx, b); // This is needed to set the time factors

tester.mScheme.Predict(tester.GetModelPart(), dof_set, A, Dx, b);
constexpr auto expected_water_pressure = 13.0;

KRATOS_EXPECT_DOUBLE_EQ(tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(WATER_PRESSURE, 0),
expected_water_pressure);
}

KRATOS_TEST_CASE_IN_SUITE(NewmarkUPwSchemePredict_PredictsWaterPressureWhenDerivativeIsFixed,
KratosGeoMechanicsFastSuiteWithoutKernel)
{
NewmarkQuasistaticUPwSchemeTester tester;

ModelPart::DofsArrayType dof_set;
CompressedMatrix A;
Vector Dx;
Vector b;
tester.GetModelPart().GetProcessInfo()[DELTA_TIME] = 4.0;
tester.GetModelPart().GetNode(0).FastGetSolutionStepValue(DT_WATER_PRESSURE, 1) = 3.0;
tester.GetModelPart().GetNode(0).FastGetSolutionStepValue(DT_WATER_PRESSURE, 0) = 7.0;
tester.GetModelPart().GetNode(0).Fix(DT_WATER_PRESSURE);
tester.mScheme.InitializeSolutionStep(tester.GetModelPart(), A, Dx, b); // This is needed to set the time factors

tester.mScheme.Predict(tester.GetModelPart(), dof_set, A, Dx, b);
constexpr auto expected_water_pressure = 25.0;

KRATOS_EXPECT_DOUBLE_EQ(tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(WATER_PRESSURE, 0),
expected_water_pressure);
}

KRATOS_TEST_CASE_IN_SUITE(NewmarkUPwSchemePredict_UpdatesVariablesDerivatives, KratosGeoMechanicsFastSuiteWithoutKernel)
{
NewmarkQuasistaticUPwSchemeTester tester;
Expand All @@ -102,6 +143,7 @@ KRATOS_TEST_CASE_IN_SUITE(NewmarkUPwSchemePredict_UpdatesVariablesDerivatives, K
CompressedMatrix A;
Vector Dx;
Vector b;
tester.GetModelPart().GetNode(0).FastGetSolutionStepValue(DT_WATER_PRESSURE, 1) = 3.0;

tester.mScheme.InitializeSolutionStep(tester.GetModelPart(), A, Dx, b); // This is needed to set the time factors

Expand All @@ -110,7 +152,7 @@ KRATOS_TEST_CASE_IN_SUITE(NewmarkUPwSchemePredict_UpdatesVariablesDerivatives, K
// These expected numbers result from the calculations in UpdateVariablesDerivatives
const auto expected_acceleration = Kratos::array_1d<double, 3>{0.0, 0.0, 0.0};
const auto expected_velocity = Kratos::array_1d<double, 3>{0.0, 0.0, 0.0};
constexpr auto expected_dt_water_pressure = 1.0 / 3.0;
constexpr auto expected_dt_water_pressure = 3.0;

const auto actual_acceleration =
tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(ACCELERATION, 0);
Expand Down
Loading