From 3ecd8ee8c9fe789f5cd0acb5116639bf23ffcc32 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Mon, 28 Oct 2024 17:08:40 +0100 Subject: [PATCH] [GeoMechanicsApplication] Clean up strategies (#12792) Changes include: - Removed a redundant overridden member. - Removed unused function parameters (all related to linear solvers) from the constructors of the strategies. The base class of all our strategies (`GeoMechanicsNewtonRaphsonStrategy`) used to receive a pointer to a linear solver, but that parameter was never used. Consequently, there is no point in supplying one. In fact, requesting one but not using it is very confusing. The constructors of the derived classes simply forwarded the linear solver pointers to the base class constructor. Also there, the corresponding parameters have been removed from the constructors. - Moved a few member functions and data members that were defined in class `GeoMechanicsNewtonRaphsonStrategy` to one of the derived classes (`GeoMechanicsRammArcLengthStrategy`), since only there they were actually being used. As a result, the moved members could be made `private` rather than `protected`. - Replaced a local implementation of the L2 norm calculation by calling Boost's `norm_2` function. There is one thing we need to keep in mind here: the local implementation used to be done in parallel (although it was using deprecated functionality to achieve that). In other words, there might be a performance impact here. The changes made also fix three code smells found by SonarQube. --- .../add_custom_strategies_to_python.cpp | 6 +- ...ewton_raphson_erosion_process_strategy.hpp | 12 +--- .../geo_mechanics_newton_raphson_strategy.hpp | 72 +------------------ ...geo_mechanics_ramm_arc_length_strategy.hpp | 56 ++++++++++----- .../solving_strategy_factory.hpp | 2 +- .../custom_workflows/dgeoflow.cpp | 4 +- .../python_scripts/geomechanics_solver.py | 3 - 7 files changed, 50 insertions(+), 105 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp index a985ff46abcf..2e1e5183b4b4 100644 --- a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp +++ b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp @@ -110,18 +110,18 @@ void AddCustomStrategiesToPython(pybind11::module& m) py::class_( m, "GeoMechanicsNewtonRaphsonStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsNewtonRaphsonErosionProcessStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsRammArcLengthStrategy") - .def(py::init()) .def("UpdateLoads", &GeoMechanicsRammArcLengthStrategyType::UpdateLoads); diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp index 2431333c34f5..aca6fca4a05d 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp @@ -54,7 +54,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy GeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart& model_part, typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -63,16 +62,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { rank = model_part.GetCommunicator().MyPID(); mPipingIterations = rParameters["max_piping_iterations"].GetInt(); diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp index f8360def0513..c9017fe3965d 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp @@ -52,16 +52,14 @@ class GeoMechanicsNewtonRaphsonStrategy * @brief Default constructor * @param rModelPart The model part of the problem * @param pScheme The integration scheme - * @param pNewLinearSolver The linear solver employed * @param pNewConvergenceCriteria The convergence criteria employed * @param MaxIterations The maximum number of iterations * @param CalculateReactions The flag for the reaction calculation * @param ReformDofSetAtEachStep The flag that allows to compute the modification of the DOF * @param MoveMeshFlag The flag that allows to move the mesh */ - GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -70,15 +68,7 @@ class GeoMechanicsNewtonRaphsonStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : ResidualBasedNewtonRaphsonStrategy( - model_part, - pScheme, - /*pNewLinearSolver,*/ - pNewConvergenceCriteria, - pNewBuilderAndSolver, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { // only include validation with c++11 since raw_literals do not exist in c++03 Parameters default_parameters(R"( @@ -100,62 +90,6 @@ class GeoMechanicsNewtonRaphsonStrategy // Validate against defaults -- this also ensures no type mismatch rParameters.ValidateAndAssignDefaults(default_parameters); - - // Set Load SubModelParts and Variable names - if (rParameters["loads_sub_model_part_list"].size() > 0) { - mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); - mVariableNames.resize(rParameters["loads_variable_list"].size()); - - if (mSubModelPartList.size() != mVariableNames.size()) - KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" - << std::endl; - - for (unsigned int i = 0; i < mVariableNames.size(); ++i) { - mSubModelPartList[i] = - &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); - mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); - } - } - } - -protected: - /// Member Variables - std::vector mSubModelPartList; /// List of every SubModelPart associated to an external load - std::vector mVariableNames; /// Name of the nodal variable associated to every SubModelPart - - int Check() override - { - KRATOS_TRY - - return MotherType::Check(); - - KRATOS_CATCH("") - } - - double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) - { - double ReferenceDofsNorm = 0.0; - - int NumThreads = ParallelUtilities::GetNumThreads(); - OpenMPUtils::PartitionVector DofSetPartition; - OpenMPUtils::DivideInPartitions(rDofSet.size(), NumThreads, DofSetPartition); - -#pragma omp parallel reduction(+ : ReferenceDofsNorm) - { - int k = OpenMPUtils::ThisThread(); - - typename DofsArrayType::iterator DofsBegin = rDofSet.begin() + DofSetPartition[k]; - typename DofsArrayType::iterator DofsEnd = rDofSet.begin() + DofSetPartition[k + 1]; - - for (typename DofsArrayType::iterator itDof = DofsBegin; itDof != DofsEnd; ++itDof) { - if (itDof->IsFree()) { - const double& temp = itDof->GetSolutionStepValue(); - ReferenceDofsNorm += temp * temp; - } - } - } - - return sqrt(ReferenceDofsNorm); } }; // Class GeoMechanicsNewtonRaphsonStrategy diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp index 2ecd12520b26..e13287c0bcaf 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp @@ -19,6 +19,9 @@ // Application includes #include "geo_mechanics_application_variables.h" +#include +#include + namespace Kratos { @@ -50,12 +53,9 @@ class GeoMechanicsRammArcLengthStrategy using GrandMotherType::mpDx; // Delta x of iteration i using GrandMotherType::mpScheme; using GrandMotherType::mReformDofSetAtEachStep; - using MotherType::mSubModelPartList; - using MotherType::mVariableNames; - GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -64,17 +64,24 @@ class GeoMechanicsRammArcLengthStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { + // Set Load SubModelParts and Variable names + if (rParameters["loads_sub_model_part_list"].size() > 0) { + mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); + mVariableNames.resize(rParameters["loads_variable_list"].size()); + + if (mSubModelPartList.size() != mVariableNames.size()) + KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" + << std::endl; + + for (unsigned int i = 0; i < mVariableNames.size(); ++i) { + mSubModelPartList[i] = + &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); + mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); + } + } + mDesiredIterations = rParameters["desired_iterations"].GetInt(); mMaxRadiusFactor = rParameters["max_radius_factor"].GetDouble(); mMinRadiusFactor = rParameters["min_radius_factor"].GetDouble(); @@ -306,7 +313,7 @@ class GeoMechanicsRammArcLengthStrategy else if (mRadius < mMinRadiusFactor * mRadius_0) mRadius = mMinRadiusFactor * mRadius_0; // Update Norm of x - mNormxEquilibrium = this->CalculateReferenceDofsNorm(rDofSet); + mNormxEquilibrium = CalculateReferenceDofsNorm(rDofSet); } else { std::cout << "************ NO CONVERGENCE: restoring equilibrium path ************" << std::endl; @@ -534,6 +541,23 @@ class GeoMechanicsRammArcLengthStrategy // Save the applied Lambda factor mLambda_old = mLambda; } + +private: + static double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) + { + auto is_free_dof = [](const auto& rDof) { return rDof.IsFree(); }; + auto free_dofs = rDofSet | boost::adaptors::filtered(is_free_dof); + + auto free_dof_values = Vector{boost::size(free_dofs)}; + auto get_dof_value = [](const auto& rDof) { return rDof.GetSolutionStepValue(); }; + std::transform(std::begin(free_dofs), std::end(free_dofs), free_dof_values.begin(), get_dof_value); + + return norm_2(free_dof_values); + } + + std::vector mSubModelPartList; // List of every SubModelPart associated to an external load + std::vector mVariableNames; // Name of the nodal variable associated to every SubModelPart + }; // Class GeoMechanicsRammArcLengthStrategy } // namespace Kratos \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp index 03296529d279..6e21ce30c7c4 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp @@ -69,7 +69,7 @@ class SolvingStrategyFactory ParametersUtilities::CopyOptionalParameters(rSolverSettings, strategy_entries); auto result = std::make_unique>( - rModelPart, scheme, solver, criteria, builder_and_solver, strategy_parameters, + rModelPart, scheme, criteria, builder_and_solver, strategy_parameters, max_iterations, calculate_reactions, reform_dof_set_at_each_step, move_mesh_flag); result->SetEchoLevel(echo_level); return result; diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp index dbc7c4350839..e0abc0639d35 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp @@ -138,8 +138,8 @@ KratosExecute::GeoMechanicsNewtonRaphsonErosionProcessStrategyType::Pointer Krat auto pSolvingStrategy = Kratos::make_unique>( - rModelPart, p_scheme, p_solver, p_criteria, p_builder_and_solver, p_parameters, - MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); + rModelPart, p_scheme, p_criteria, p_builder_and_solver, p_parameters, MaxIterations, + CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); pSolvingStrategy->Check(); return pSolvingStrategy; diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py index c5580cb9f167..4a0d43e24856 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py @@ -489,7 +489,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -504,7 +503,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("max_piping_iterations", self.settings["max_piping_iterations"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonErosionProcessStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -545,7 +543,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsRammArcLengthStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params,