[GeomechanicsApplication] Geo/linear elastic strategy #12479

aronnoordam · 2024-06-25T13:16:38Z

added another strategy and builder and solver meant for linear elastic problems.

Mostly, this strategy exploits the fact that the internal forces dont have to be recalculated every iteration.

stiffness matrix is constant
damping matrix is constant
mass matrix is constant
rhs, only condition contributions are updated
the linear system is pre factorized in the first iteration/ time-step, if possible

Only external forces are updated

aronnoordam · 2024-07-08T07:52:41Z

applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py

dont forget to revert this

# Conflicts: # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.hpp # applications/GeoMechanicsApplication/custom_python/add_custom_processes_to_python.cpp # applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp # applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py

# Conflicts: # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp

rfaasse

Hi @aronnoordam, thanks for processing the first round of feedback. This helps a lot in understanding the code and giving better feedback and in my opinion both the tests and production code improved a lot: less duplication, clearer responsibilities of classes and more specific testing on what the strategy does and (more importantly) does not do/call.

This round, I tried to be more detailed in my comments, due to larger scale design becoming more clear. For the builder and solver, I would like to discuss a bit more with the team about if it makes sense to move more newmark functionality to the scheme and reduce a bit more duplication (which is hard as you mentioned, due to the differences in private non-virtual functions).

Next to that, @avdg81 and @WPK4FEM will also have another look at the PR 👍

rfaasse · 2024-10-10T07:03:52Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+        : ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(
+              rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, false, MoveMeshFlag)
+    {
+        // new constructor


This comment can be removed

rfaasse · 2024-10-10T07:06:09Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+
+        BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);
+
+        // Note that constraints are not applied in this predict, nor is an update performed


Maybe it's good to mention here that the constraints are applied in the builder and solver.

And did you already find what the reason is that the constraints cannot be applied in the same way as for the base strategy?

yes actually, as for constraints, there is a dotproduct with the constraint transformation matrix and the RHS, now in my case, i already transformed the mass and damping matrix, thus when applying a constraint again to the rhs, part of the RHS is transformed twice.

rfaasse · 2024-10-10T07:10:37Z

...chanicsApplication/custom_strategies/schemes/incremental_newmark_linear_elastic_U_scheme.hpp

Since the formulations are different for this new incremental scheme, this should be documented in the README.md with the other schemes.

rfaasse · 2024-10-10T07:12:45Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+        }
+
+        if (is_converged) {
+            // here only the derivatives are updated


I think it's good to add here that only the derivatives are updated in the specific scheme we're using. In general schemes will update both the derivatives and the instance variable

rfaasse · 2024-10-10T07:24:40Z

...rategies/builder_and_solvers/residualbased_block_builder_and_solver_linear_elastic_dynamic.h

+    void AddDynamicsToLhs(TSystemMatrixType& rA, const ModelPart& rModelPart)
+    {
+        const double delta_time = rModelPart.GetProcessInfo()[DELTA_TIME];
+
+        double*       a_values = rA.value_data().begin();
+        const double* m_values = mMassMatrix.value_data().begin();
+        const double* c_values = mDampingMatrix.value_data().begin();
+
+        // add mass and damping contribution to LHS sparse matrix
+        // mass contribution: 1.0 / (mBeta * delta_time * delta_time) * M
+        // damping contribution: mGamma / (mBeta * delta_time) * C
+        for (std::size_t i = 0; i < rA.size1(); i++) {
+            const std::size_t col_begin = rA.index1_data()[i];
+            const std::size_t col_end   = rA.index1_data()[i + 1];
+
+            for (std::size_t j = col_begin; j < col_end; ++j) {
+                a_values[j] += (1.0 / (mBeta * delta_time * delta_time)) * m_values[j];
+                a_values[j] += (mGamma / (mBeta * delta_time)) * c_values[j];
+            }
+        }
+    }


In my opinion the code is much more readable (we don't need the comments, because the code itself just shows the formula), so could we just try it? If the performance drops significantly, we can revert it, but it's worth checking in case it is similar.

rfaasse · 2024-10-10T13:56:55Z

...pplication/tests/test_1d_wave_prop_drained_soil_linear_elastic_solver/ProjectParameters.json

+        "rayleigh_m":                         0.02,
+        "rayleigh_k":                         6e-6,
+        "strategy_type":                      "newton_raphson_linear_elastic",
+        "convergence_criterion":              "displacement_criterion",


Indeed, but there are still functions copied from the base and moved around, so I think it's a good idea to use residual criterion for one of the test cases (shouldn't change the results anyway, just a json change), just to make sure. But I leave it up to you 👍

rfaasse · 2024-10-11T07:18:04Z

applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/geo_custom_element.cpp

Could you also rename the file-name (.cpp and .h) to geo_mock_element to align with the class name?

rfaasse · 2024-10-11T07:19:54Z

applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/geo_custom_element.cpp

+//
+//  Main authors:    Aron Noordam
+//
+#include "custom_utilities/dof_utilities.h"


Could we just include dof.h instead of dof_utilities? or are we using something from the utils?

rfaasse · 2024-10-11T07:25:57Z

...ions/GeoMechanicsApplication/tests/cpp_tests/test_newton_raphson_strategy_linear_elastic.cpp

+    std::vector<double> expected_displacement_x = {0.00673, 0.0505, 0.189, 0.485, 0.961, 1.58,
+                                                   2.23,    2.76,   3.0,   2.85,  2.28,  1.40};
+    std::vector<double> expected_displacement_y = {0.364, 1.35, 2.68, 4.00, 4.95, 5.34,
+                                                   5.13,  4.48, 3.64, 2.90, 2.44, 2.31};
+
+    double delta_time                     = 0.28;
+    bool   use_iterations                 = false;
+    bool   calculate_initial_acceleration = true;
+    bool   use_direct_solver              = true;


These can all be const (same holds for the other tests)

rfaasse · 2024-10-11T07:31:23Z

...ions/GeoMechanicsApplication/tests/cpp_tests/test_newton_raphson_strategy_linear_elastic.cpp

+    KRATOS_EXPECT_EQ(geo_custom_element.GetCountCalculateLeftHandSideCalled(), 1);
+    KRATOS_EXPECT_EQ(geo_custom_element.GetCountCalculateMassMatrixCalled(), 1);
+    KRATOS_EXPECT_EQ(geo_custom_element.GetCalculateDampingMatrixCalled(), 1);
+    KRATOS_EXPECT_EQ(geo_custom_element.GetCountCalculateRightHandSideCalled(), 0);
+
+    KRATOS_EXPECT_EQ(geo_custom_condition.GetCountCalculateLeftHandSideCalled(), 1);
+    KRATOS_EXPECT_EQ(geo_custom_condition.GetCountCalculateMassMatrixCalled(), 1);
+    KRATOS_EXPECT_EQ(geo_custom_condition.GetCalculateDampingMatrixCalled(), 1);
+
+    // rhs for conditions is called each solution step and during initialisation
+    if (UseIterations) {
+        // Two iterations are performed per solution step, thus rhs for conditions is called twice each solution step step and during initialisation
+        KRATOS_EXPECT_EQ(geo_custom_condition.GetCountCalculateRightHandSideCalled(), n_steps * 2 + 1);
+    } else {
+        KRATOS_EXPECT_EQ(geo_custom_condition.GetCountCalculateRightHandSideCalled(), n_steps + 1);
+    }


Thanks for adding these checks, it also functions as documentation (it's really clear what is called how many times from this test)

avdg81

This is a big piece of work that, as I understood from you, yields a significant performance improvement, which is great. As you can see, I have several suggestions. For me the most important ones are aimed at making the new scheme's initialization and finalization behavior symmetric. Also, I have a suggestion to remove the duplication of the nonlinear iteration, which nicely ties in with my previous suggestion. I hope my suggestions will contribute to making your code clearer and more maintainable. Feel free to reach out when my comments are unclear or when you believe that some of them are not feasible. Thanks.

avdg81 · 2024-10-11T07:47:45Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+        typename TLinearSolver::Pointer            pNewLinearSolver,
+        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
+        typename TBuilderAndSolverType::Pointer    pNewBuilderAndSolver,
+        Parameters&                                rParameters,


In addition, also pNewLinearSolver is no longer (?) used. So I would suggest to remove it as well (including the corresponding Doxygen description).

avdg81 · 2024-10-11T09:21:15Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+     */
+    void Initialize() override
+    {
+        KRATOS_TRY;


Please remove the empty statement:

Suggested change

KRATOS_TRY;

KRATOS_TRY

avdg81 · 2024-10-11T09:21:42Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+            // initialize the system matrices and the initial second derivative
+            this->InititalizeSystemAndState();
+
+        KRATOS_CATCH("");


Please remove the empty statement:

Suggested change

KRATOS_CATCH("");

KRATOS_CATCH("")

avdg81 · 2024-10-11T09:28:23Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+    {
+        KRATOS_TRY;
+
+        BaseType::Initialize();


The implementation of the base class's Initialize suggests that it may be called more than once, since there is a flag (mInitializeWasPerformed) that makes sure the actual initialization is performed only once. Would we need a similar guarding mechanism here?

avdg81 · 2024-10-11T09:29:11Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+
+        BaseType::Initialize();
+
+        // Note that FindNeighbourElementsOfConditionsProcess and DeactivateConditionsOnInactiveElements are required to be perfomed before initializing the System and State


A typo:

Suggested change

// Note that FindNeighbourElementsOfConditionsProcess and DeactivateConditionsOnInactiveElements are required to be perfomed before initializing the System and State

// Note that FindNeighbourElementsOfConditionsProcess and DeactivateConditionsOnInactiveElements are required to be performed before initializing the System and State

avdg81 · 2024-10-11T09:51:20Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+
+        DofsArrayType& r_dof_set = BaseType::GetBuilderAndSolver()->GetDofSet();
+
+        BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);


Since rA, rDx, and rb are only used once, I'd suggest to eliminate this local variables:

Suggested change

BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);

BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, *BaseType::mpA, *BaseType::mpDx, *BaseType::mpb);

You could take this one step further and eliminate r_dof_set in the same way.

avdg81 · 2024-10-21T09:10:16Z

...chanicsApplication/custom_strategies/schemes/incremental_newmark_linear_elastic_U_scheme.hpp

+        KRATOS_CATCH("")
+    }
+
+    void FinalizeNonLinIteration(ModelPart& rModelPart, TSystemMatrixType&, TSystemVectorType&, TSystemVectorType&) override


I would suggest to make the initialization and finalization of a nonlinear iteration symmetric. At present, it isn't, because the member function that initializes a nonlinear iteration initializes both the elements and the conditions. (See GeoMechanicsTimeIntegrationScheme<...>::InitializeNonLinIteration). My suggestion would be to add an override for InitializeNonLinIteration that only initializes the conditions. Then in your new strategy, you can initialize the elements prior to carrying out the first nonlinear iteration.

avdg81 · 2024-10-21T09:20:33Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+            block_for_each(r_model_part.Conditions(), [&r_current_process_info](Condition& r_condition) {
+                if (r_condition.IsActive()) {
+                    r_condition.FinalizeNonLinearIteration(r_current_process_info);
+                }
+            });


If I'm not mistaken, this is exactly what your new scheme does...? So I think you can replace it by:

Suggested change

block_for_each(r_model_part.Conditions(), [&r_current_process_info](Condition& r_condition) {

if (r_condition.IsActive()) {

r_condition.FinalizeNonLinearIteration(r_current_process_info);

}

});

p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb);

avdg81 · 2024-10-21T09:26:38Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+    bool PerformIterationCycle(TSystemMatrixType&   rA,
+                               TSystemVectorType&   rDx,
+                               TSystemVectorType&   rb,
+                               TSystemVectorType&   rDxTot,
+                               std::vector<Vector>& rNonconvergedSolutions,
+                               unsigned int&        rIterationNumber)


I would suggest to extract a member function that does a single nonlinear iteration, and then reuse that for the first iteration as well as for any subsequent iterations. In terms of your implementation, it would roughly look like this:

r_model_part.GetProcessInfo()[NL_ITERATION_NUMBER] = rIterationNumber; p_scheme->InitializeNonLinIteration(r_model_part, rA, rDx, rb); BaseType::mpConvergenceCriteria->InitializeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb); is_converged = BaseType::mpConvergenceCriteria->PreCriteria(r_model_part, r_dof_set, rA, rDx, rb); // call the linear system solver to find the correction mDx for the // it is not called if there is no system to solve if (SparseSpaceType::Size(rDx) != 0) { TSparseSpace::SetToZero(rDx); TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildRHSAndSolve(p_scheme, r_model_part, rA, rDx, rb); } else { KRATOS_WARNING("NO DOFS") << "ATTENTION: no free DOFs!! " << std::endl; } // Debugging info BaseType::EchoInfo(rIterationNumber); // Updating the results stored in the database this->UpdateSolutionStepValue(rDx, rDxTot); p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb); BaseType::mpConvergenceCriteria->FinalizeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb); if (BaseType::mStoreNonconvergedSolutionsFlag) { Vector ith; BaseType::GetCurrentSolution(r_dof_set, ith); rNonconvergedSolutions.push_back(ith); } if (is_converged) { if (BaseType::mpConvergenceCriteria->GetActualizeRHSflag()) { TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb); } is_converged = BaseType::mpConvergenceCriteria->PostCriteria(r_model_part, r_dof_set, rA, rDx, rb); }

This will take away quite a bit of duplication.

avdg81 · 2024-10-21T09:29:40Z

...ustom_strategies/strategies/residualbased_newton_raphson_strategy_linear_elastic_dynamic.hpp

+        // only initialize conditions, not that this cannot be put in the scheme, as the scheme is required to InitializeNonLinearIteration both elements and conditions
+        const auto& r_current_process_info = r_model_part.GetProcessInfo();
+        block_for_each(r_model_part.Conditions(), [&r_current_process_info](Condition& r_condition) {
+            if (r_condition.IsActive()) {
+                r_condition.InitializeNonLinearIteration(r_current_process_info);
+            }
+        });


I would suggest to do the initialization of the elements here, and leave the initialization of the conditions to new your new scheme (see one of my other suggestions). That would make you new scheme symmetric in terms of initialization and finalization. But it also highlights that the handling of elements is special and needs to be done only once rather than in each nonlinear iteration.

…vent calculation of internal forces

…calculation

…uasi static scheme

# Conflicts: # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.hpp

# Conflicts: # applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp # applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.h # applications/GeoMechanicsApplication/custom_strategies/schemes/backward_euler_quasistatic_U_Pw_scheme.hpp # applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_inremental_linear_elastic_3D_law.cpp # applications/GeoMechanicsApplication/tests/cpp_tests/custom_strategies/test_backward_euler_quasistatic_U_Pw_scheme.cpp

aronnoordam added 6 commits June 21, 2024 17:27

added linear elastic dynamic strategy

9c45b7d

improved linear elastic dynamic strategy

0a862d3

added clarification

1f1e3e0

added calculation of initial acceleration

507aeb9

connected solver in python

346e221

reverted unintended commit

c79fe65

aronnoordam commented Jul 8, 2024

View reviewed changes

applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py Outdated

Copy link

Member Author

aronnoordam Jul 8, 2024

Choose a reason for hiding this comment

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

dont forget to revert this

aronnoordam added 16 commits August 2, 2024 12:02

added find_neighbour_elements and solved issue in validation

19f2f70

cleanup

d15b379

added benchmark test

dd32a7f

added benchmark test

74aead2

added dynamic test for multistage

fee5a20

clang tidy custom processes and strategies

b5c7729

clang tidy solver linear elastic

23473f7

clang tidy test solver linear elastic

aaf13ac

clang tidy test solver linear elastic

17f6f8e

removed forcing of using direct solver

8979322

corrected test after merge

e4a1e45

resolved warnings and added test for iterative solver

dd6ab80

removed unused variable

4fce2a9

resolved errors in CI

a678257

resolved errors in CI

4a172e1

aronnoordam self-assigned this Sep 13, 2024

aronnoordam added 6 commits September 13, 2024 09:53

resolved errors in CI

a5accd0

cleaned test_dynamics

236def6

cleaned test_dynamics

94ea951

sonar cloud suggestions

93e2dcb

Merge branch 'master' into geo/linear_elastic_strategy

810d3c1

# Conflicts: # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp

resolved double declaration

de657ba

aronnoordam added 2 commits October 10, 2024 12:00

resolved sonarcloud issues

6cd438d

removed unused variable

b4f9be7

rfaasse reviewed Oct 11, 2024

View reviewed changes

avdg81 reviewed Oct 21, 2024

View reviewed changes

aronnoordam added 26 commits October 24, 2024 15:06

solved bug for multistage and resolved comments

cbeb8c4

rhs of elements are now called once, with 0 zero displacements to pre…

7e391af

…vent calculation of internal forces

added test for gravity load on soil column and fixed bug in init acc …

70bb353

…calculation

added test for gravity load on soil column and fixed bug in init acc …

a71e1ca

…calculation

corrected c++ test

7aef93c

internal forces calculation can now be skipped in lysmer

3f2b0c3

moved initialisation of dxTot

1335743

added test to check if solvesolutionstep can be called multiple times

2a3dee6

removed calculation of velocity and acceleration for backward euler q…

cb32891

…uasi static scheme

commented meaningless tests

fe4a6c1

added checks if rotation dofs are added

9593e2b

added validate and assign defaults for absorbing boundary

80d599b

added skip mark for settlement test

38c788d

Merge branch 'master' into geo/linear_elastic_strategy

91d79cd

# Conflicts: # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.cpp # applications/GeoMechanicsApplication/custom_conditions/U_Pw_normal_lysmer_absorbing_condition.hpp

added function to reset incremental linear elastic law

9941a9c

style

efc79de

corrected max_iterations call

3c6aefc

added incremental newmark formulations to the README.md

6103040

cleanup

b0bec9f

corrected incremental newmark equilibrium equation

3ce7ea6

cleanup

4981c3e

Merge branch 'master' into geo/linear_elastic_strategy

1838cfa

cleanup

3062f9f

added try, catch statements

ea03fcd

solved index error

a13b1ec

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[GeomechanicsApplication] Geo/linear elastic strategy #12479

[GeomechanicsApplication] Geo/linear elastic strategy #12479

aronnoordam commented Jun 25, 2024 •

edited

Loading

aronnoordam Jul 8, 2024

rfaasse left a comment

rfaasse Oct 10, 2024

rfaasse Oct 10, 2024

rfaasse Oct 10, 2024

aronnoordam Oct 11, 2024

rfaasse Oct 10, 2024

rfaasse Oct 10, 2024

rfaasse Oct 10, 2024

rfaasse Oct 10, 2024

rfaasse Oct 11, 2024

rfaasse Oct 11, 2024

aronnoordam Oct 11, 2024

rfaasse Oct 11, 2024

rfaasse Oct 11, 2024

avdg81 left a comment

avdg81 Oct 11, 2024

avdg81 Oct 11, 2024

avdg81 Oct 11, 2024

avdg81 Oct 11, 2024

avdg81 Oct 11, 2024

avdg81 Oct 11, 2024

avdg81 Oct 21, 2024

avdg81 Oct 21, 2024

avdg81 Oct 21, 2024

avdg81 Oct 21, 2024


		BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);

		// Note that constraints are not applied in this predict, nor is an update performed


		BaseType::Initialize();

		// Note that FindNeighbourElementsOfConditionsProcess and DeactivateConditionsOnInactiveElements are required to be perfomed before initializing the System and State


		DofsArrayType& r_dof_set = BaseType::GetBuilderAndSolver()->GetDofSet();

		BaseType::GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);

[GeomechanicsApplication] Geo/linear elastic strategy #12479

Are you sure you want to change the base?

[GeomechanicsApplication] Geo/linear elastic strategy #12479

Conversation

aronnoordam commented Jun 25, 2024 • edited Loading

Choose a reason for hiding this comment

rfaasse left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

avdg81 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aronnoordam commented Jun 25, 2024 •

edited

Loading