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 all 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
94 changes: 93 additions & 1 deletion ...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 @@ -180,7 +192,87 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
}
}

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

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

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

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

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&, TSystemMatrixType&, TSystemVectorType&, TSystemVectorType&) 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
7 changes: 7 additions & 0 deletions ...ons/GeoMechanicsApplication/custom_strategies/schemes/newmark_quasistatic_U_Pw_scheme.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,13 @@ class NewmarkQuasistaticUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace
KRATOS_CATCH("")
}

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

} // namespace Kratos
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