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

rfaasse · 2025-05-20T08:51:32Z

📝 Description
This PR adds prediction to all Newmark schemes, both for vector as well as scalar variables.

…ot fixed (and changed other tests accordingly)

…ixed + formatting

rfaasse · 2025-05-20T09:37:30Z

applications/GeoMechanicsApplication/custom_strategies/schemes/generalized_newmark_scheme.hpp

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


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

markelov208

Hi Richard, I think it was difficult to spot such a thing. I have only a comment on Predict function's description in common_strategies\schemes\README.md. Shall it be updated?

rfaasse · 2025-05-20T12:59:43Z

applications/GeoMechanicsApplication/custom_strategies/schemes/generalized_newmark_scheme.hpp

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


@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?

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.

rfaasse · 2025-05-20T13:01:46Z

Hi Richard, I think it was difficult to spot such a thing. I have only a comment on Predict function's description in common_strategies\schemes\README.md. Shall it be updated?

Good catch, I'll update the README

RiccardoRossi

hi, i hope you don't mind if i chime in in this conversation.

in order for the Predict to work correctly, we need to guarantee consistency between the DISPLACEMENT, VELOCITY and ACCELERATION (i'll focus my the case on standard solid elements...of course this can be expanded to ROTATION, ANGULAR_VELOCITY and ANGULAR_ACCELERATION if needed.

the point is that the standard "Update" method, called after every system solution guarantees that the VELOCITY and ACCELERATION are updated in accordance to the provided DISPLACEMENT.
One way to undersand the "Predict" is to consider it as a "cheap approximation of a Solve", that is, the Predict provides an approximation to the end of step DISPLACEMENT from which VEL and ACC can be computed consistently.

Such approximation is exact where the DISPLACEMENT is Fixed (dirichlet constraint).

The problem comes when one wants to impose VELOCITY or ACCELERATION as a fixed value instead of the DISPLACEMENT. The current method does not support this, "as is".

In my opinion the correct way of handling this is to guarantee that when the VELOCITY or the ACCELERATION is fixed the DISPLACEMENT is also fixed automaticially, even if not done in the input (of course we will need to unfix it in the InitializeSolutionStep when fixing it automatically). The DISPLACEMENT on the automatically fixed dofs should then be provided the value that guarantees that the desired value of ACCELERATION and/or VELOCITY is obtained.

RiccardoRossi

ok, looking into the PR i think my comment should have been shorter:

"artificially fix DISPLACEMENT_" whenever you want to fix "VELOCITY_" or "ACCELERATION_*" and then simply apply call "Update"
(and then unfix the list of dofs you fix this way in the FinalizeSolutionStep)

i cannot say if the formulas you are using for the applied displacement are the correct ones ... but provided that they are, this should work

…ict displacements based on its derivatives (only water pressure)

rfaasse · 2025-06-04T07:02:59Z

ok, looking into the PR i think my comment should have been shorter:

"artificially fix DISPLACEMENT__" whenever you want to fix "VELOCITY__" or "ACCELERATION_*" and then simply apply call "Update" (and then unfix the list of dofs you fix this way in the FinalizeSolutionStep)

i cannot say if the formulas you are using for the applied displacement are the correct ones ... but provided that they are, this should work

Thanks for your comment, indeed at this point we need to fix the displacement dof manually when we want to fix the velocity or acceleration. This was also what was missing in our prescribed derivatives tests that were failing before when changing the order of Predict and InitializeSolutionStep (related to #13350).

We now improved that test by manually fixing displacements (see #13478 if you're interested), which is why we don't see any issues anymore in our test suite with switching the order of Predict and InitializeSolutionStep.

At some point indeed it would be nice to be able to automatically fix displacements when prescribing its (second) derivative, but for now this is good enough for us

…ariable-prediction-in-all-newmark-schemes

rfaasse added 5 commits May 19, 2025 09:58

Moved the predict from dynamic to the generalized newmark scheme

f795b45

Added prediction of first order scalar variable when dt_variable is n…

b890fdc

…ot fixed (and changed other tests accordingly)

Added prediction of first order scalar variable when dt_variable is f…

2ccf9eb

…ixed + formatting

Fixed unit test by fixing logic (if -> else if)

13e86d3

Minor refactoring

6ccec00

rfaasse requested a review from a team as a code owner May 20, 2025 08:51

rfaasse linked an issue May 20, 2025 that may be closed by this pull request

[GeoMechanicsApplication] Include variable prediction in all Newmark schemes #11957

Closed

rfaasse commented May 20, 2025

View reviewed changes

rfaasse requested a review from markelov208 May 20, 2025 09:38

markelov208 previously approved these changes May 20, 2025

View reviewed changes

rfaasse commented May 20, 2025

View reviewed changes

RiccardoRossi reviewed May 24, 2025

View reviewed changes

Changed predict of quasistatic UPw scheme, such that it does not pred…

78b5a69

…ict displacements based on its derivatives (only water pressure)

rfaasse dismissed markelov208’s stale review via 78b5a69 June 3, 2025 13:08

Merge remote-tracking branch 'origin/master' into geo/11957-include-v…

2855780

…ariable-prediction-in-all-newmark-schemes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

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

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

rfaasse commented May 20, 2025

Uh oh!

rfaasse May 20, 2025

Uh oh!

markelov208 left a comment

Uh oh!

rfaasse May 20, 2025

Uh oh!

WPK4FEM May 20, 2025

Uh oh!

rfaasse commented May 20, 2025

Uh oh!

RiccardoRossi left a comment

Uh oh!

RiccardoRossi left a comment

Uh oh!

rfaasse commented Jun 4, 2025 •

edited

Loading

Uh oh!

Uh oh!

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

Are you sure you want to change the base?

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

Conversation

rfaasse commented May 20, 2025

Uh oh!

rfaasse May 20, 2025

Choose a reason for hiding this comment

Uh oh!

markelov208 left a comment

Choose a reason for hiding this comment

Uh oh!

rfaasse May 20, 2025

Choose a reason for hiding this comment

Uh oh!

WPK4FEM May 20, 2025

Choose a reason for hiding this comment

Uh oh!

rfaasse commented May 20, 2025

Uh oh!

RiccardoRossi left a comment

Choose a reason for hiding this comment

Uh oh!

RiccardoRossi left a comment

Choose a reason for hiding this comment

Uh oh!

rfaasse commented Jun 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

rfaasse commented Jun 4, 2025 •

edited

Loading