Merge pull request #12375 from KratosMultiphysics/core/nonconvereged_…

…solutions_for_nr_strategy

[Core] Adding a method to fetch the non-converged solutions from the NewtonRaphson strategy

kratos/python/add_strategies_to_python.cpp

@@ -661,6 +661,8 @@ namespace Kratos:: Python
             .def("GetUseOldStiffnessInFirstIterationFlag", &ResidualBasedNewtonRaphsonStrategyType::GetUseOldStiffnessInFirstIterationFlag)
             .def("SetReformDofSetAtEachStepFlag", &ResidualBasedNewtonRaphsonStrategyType::SetReformDofSetAtEachStepFlag)
             .def("GetReformDofSetAtEachStepFlag", &ResidualBasedNewtonRaphsonStrategyType::GetReformDofSetAtEachStepFlag)
+            .def("GetNonconvergedSolutions", &ResidualBasedNewtonRaphsonStrategyType::GetNonconvergedSolutions)
+            .def("SetUpNonconvergedSolutionsFlag", &ResidualBasedNewtonRaphsonStrategyType::SetUpNonconvergedSolutionsFlag)
             ;
 
         // ARC-LENGTH

kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h

@@ -879,6 +879,41 @@ class ResidualBasedNewtonRaphsonStrategy
         KRATOS_CATCH("");
     }
 
+    /**
+    * @brief Gets the current solution vector
+    * @param rDofSet The set of degrees of freedom (DOFs)
+    * @param this_solution The vector that will be filled with the current solution values
+    * @details This method retrieves the current solution values for the provided DOF set.
+    * The provided solution vector will be resized to match the size of the DOF set if necessary,
+    * and will be filled with the solution values corresponding to each DOF. Each value is accessed
+    * using the equation ID associated with each DOF.
+    */
+    void GetCurrentSolution(DofsArrayType& rDofSet, Vector& this_solution) {
+        this_solution.resize(rDofSet.size());
+        block_for_each(rDofSet, [&](const auto& r_dof) {
+            this_solution[r_dof.EquationId()] = r_dof.GetSolutionStepValue();
+        });
+    }
+
+    /**
+     * @brief Gets the matrix of non-converged solutions and the DOF set
+     * @return A tuple containing the matrix of non-converged solutions and the DOF set
+     * @details This method returns a tuple where the first element is a matrix of non-converged solutions and the second element is the corresponding DOF set.
+     */
+    std::tuple<Matrix, DofsArrayType> GetNonconvergedSolutions(){
+        return std::make_tuple(mNonconvergedSolutionsMatrix, GetBuilderAndSolver()->GetDofSet());
+    }
+
+    /**
+    * @brief Sets the state for storing non-converged solutions.
+    * @param state The state to set for storing non-converged solutions (true to enable, false to disable).
+    * @details This method enables or disables the storage of non-converged solutions at each iteration
+    * by setting the appropriate flag. When the flag is set to true, non-converged solutions will be stored.
+    */
+    void SetUpNonconvergedSolutionsFlag(bool state) {
+        mStoreNonconvergedSolutionsFlag = state;
+    }
+
     /**
      * @brief Solves the current step. This function returns true if a solution has been found, false otherwise.
      */
@@ -889,6 +924,13 @@ class ResidualBasedNewtonRaphsonStrategy
         typename TSchemeType::Pointer p_scheme = GetScheme();
         typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver();
         auto& r_dof_set = p_builder_and_solver->GetDofSet();
+        std::vector<Vector> NonconvergedSolutions;
+
+        if (mStoreNonconvergedSolutionsFlag) {
+            Vector initial;
+            GetCurrentSolution(r_dof_set,initial);
+            NonconvergedSolutions.push_back(initial);
+        }
 
         TSystemMatrixType& rA  = *mpA;
         TSystemVectorType& rDx = *mpDx;
@@ -929,6 +971,12 @@ class ResidualBasedNewtonRaphsonStrategy
         p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb);
         mpConvergenceCriteria->FinalizeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
 
+        if (mStoreNonconvergedSolutionsFlag) {
+            Vector first;
+            GetCurrentSolution(r_dof_set,first);
+            NonconvergedSolutions.push_back(first);
+        }
+
         if (is_converged) {
             if (mpConvergenceCriteria->GetActualizeRHSflag()) {
                 TSparseSpace::SetToZero(rb);
@@ -996,6 +1044,12 @@ class ResidualBasedNewtonRaphsonStrategy
             p_scheme->FinalizeNonLinIteration(r_model_part, rA, rDx, rb);
             mpConvergenceCriteria->FinalizeNonLinearIteration(r_model_part, r_dof_set, rA, rDx, rb);
 
+            if (mStoreNonconvergedSolutionsFlag == true){
+                Vector ith;
+                GetCurrentSolution(r_dof_set,ith);
+                NonconvergedSolutions.push_back(ith);
+            }
+
             residual_is_updated = false;
 
             if (is_converged == true)
@@ -1039,6 +1093,15 @@ class ResidualBasedNewtonRaphsonStrategy
         if (mCalculateReactionsFlag == true)
             p_builder_and_solver->CalculateReactions(p_scheme, r_model_part, rA, rDx, rb);
 
+        if (mStoreNonconvergedSolutionsFlag) {
+            mNonconvergedSolutionsMatrix = Matrix( r_dof_set.size(), NonconvergedSolutions.size() );
+            for (std::size_t i = 0; i < NonconvergedSolutions.size(); ++i) {
+                block_for_each(r_dof_set, [&](const auto& r_dof) {
+                    mNonconvergedSolutionsMatrix(r_dof.EquationId(), i) = NonconvergedSolutions[i](r_dof.EquationId());
+                });
+            }
+        }
+
         return is_converged;
     }
 
@@ -1265,6 +1328,20 @@ class ResidualBasedNewtonRaphsonStrategy
 
     bool mKeepSystemConstantDuringIterations; // Flag to allow keeping system matrix constant during iterations
 
+    /**
+     * @brief This matrix stores the non-converged solutions
+     * @details The matrix is structured such that each column represents the solution vector at a specific non-converged iteration.
+     */
+    Matrix mNonconvergedSolutionsMatrix;
+
+    /**
+     * @brief Flag indicating whether to store non-converged solutions
+     * @details Only when set to true (by calling the SetUpNonconvergedSolutionsGathering method) will the non-converged solutions at each iteration be stored.
+     */
+    bool mStoreNonconvergedSolutionsFlag = false;
+
+
+
     ///@}
     ///@name Private Operators
     ///@{

kratos/tests/cpp_tests/strategies/strategies/test_strategies.cpp

@@ -227,6 +227,64 @@ namespace Kratos
             }
         }
 
+
+        /**
+         * Checks if the Nonconvereged solutions of the Newton Rapshon strategy performs correctly
+         */
+        KRATOS_TEST_CASE_IN_SUITE(NonconvergedSolutionsNRStrategy, KratosCoreFastSuite)
+        {
+            Model current_model;
+
+            constexpr double tolerance_residual_criteria = 1e-15;  // with this tolerance, I expect the strategy to reach the maximum number of iterations, i.e. 10
+
+
+            ModelPart& model_part = current_model.CreateModelPart("Main");
+
+            SchemeType::Pointer pscheme = SchemeType::Pointer( new ResidualBasedIncrementalUpdateStaticSchemeType() );
+            LinearSolverType::Pointer psolver = LinearSolverType::Pointer( new SkylineLUFactorizationSolverType() );
+            ConvergenceCriteriaType::Pointer pcriteria = ConvergenceCriteriaType::Pointer( new ResidualCriteriaType(tolerance_residual_criteria, tolerance_residual_criteria) );
+            BuilderAndSolverType::Pointer pbuildandsolve = BuilderAndSolverType::Pointer( new ResidualBasedBlockBuilderAndSolverType(psolver) );
+
+            ResidualBasedNewtonRaphsonStrategyType::Pointer pstrategy = ResidualBasedNewtonRaphsonStrategyType::Pointer( new ResidualBasedNewtonRaphsonStrategyType(model_part, pscheme, pcriteria, pbuildandsolve, 10, true));
+
+            DofsArrayType Doftemp = BasicTestStrategyDisplacement(model_part, ResidualType::NON_LINEAR);
+
+            NodeType::Pointer pnode = model_part.pGetNode(1);
+
+            double time = 0.0;
+            const double delta_time = 1.0e-4;
+            const unsigned int number_iterations = 1;
+
+            pstrategy-> SetUpNonconvergedSolutionsFlag(true);
+
+            for (unsigned int iter = 0; iter < number_iterations; ++iter) {
+                time += delta_time;
+
+                model_part.CloneTimeStep(time);
+
+                array_1d<double, 3> init_vector;
+                init_vector[0] = 0.5;
+                init_vector[1] = 0.5;
+                init_vector[2] = 0.5;
+                pnode->FastGetSolutionStepValue(DISPLACEMENT) = init_vector;
+
+                pcriteria->SetEchoLevel(0);
+                pstrategy->SetEchoLevel(0);
+                pstrategy->Solve();
+
+                auto [NonconveregedSolutionsMatrix,Dofs]= pstrategy->GetNonconvergedSolutions();
+
+                unsigned int expected_rows = Dofs.size();
+                unsigned int expected_cols = model_part.GetProcessInfo()[NL_ITERATION_NUMBER] + 1; //+1 because zeroth is included
+
+                KRATOS_CHECK_EQUAL(NonconveregedSolutionsMatrix.size1(), expected_rows);
+                KRATOS_CHECK_EQUAL(NonconveregedSolutionsMatrix.size2(), expected_cols);
+
+
+            }
+        }
+
+
         KRATOS_TEST_CASE_IN_SUITE(LineSearchStrategy, KratosCoreFastSuite)
         {
             Model current_model;