From f14d733a2e771bf874d2b6fe1221a66f4d10213e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Wed, 6 Mar 2019 10:08:19 +0100 Subject: [PATCH 01/14] [Core] Residual criteria and block B&S correction --- .../convergencecriterias/residual_criteria.h | 55 +++++++++++++------ .../residualbased_newton_raphson_strategy.h | 47 ++++++++-------- 2 files changed, 59 insertions(+), 43 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index 9f9c72aa42f..b3d588ea6a4 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -89,11 +89,35 @@ class ResidualCriteria ///@name Life Cycle ///@{ + //* Constructor. + explicit ResidualCriteria(Kratos::Parameters Settings) + : BaseType() + { + if (Settings.Has("residual_absolute_tolerance")) { + mAlwaysConvergedNorm = Settings["residual_absolute_tolerance"].GetDouble(); + } else if (Settings.Has("absolute_tolerance")) { + mAlwaysConvergedNorm = Settings["absolute_tolerance"].GetDouble(); + } else { + KRATOS_WARNING("ResidualCriteria") << "residual_absolute_tolerance or absolute_tolerance nor defined on settings. Using default 1.0e-9" << std::endl; + mAlwaysConvergedNorm = 1.0e-9; + } + if (Settings.Has("residual_relative_tolerance")) { + mRatioTolerance = Settings["residual_relative_tolerance"].GetDouble(); + } else if (Settings.Has("relative_tolerance")) { + mRatioTolerance = Settings["relative_tolerance"].GetDouble(); + } else { + KRATOS_WARNING("ResidualCriteria") << "residual_relative_tolerance or relative_tolerance nor defined on settings. Using default 1.0e-4" << std::endl; + mRatioTolerance = 1.0e-4; + } + } + //* Constructor. explicit ResidualCriteria( TDataType NewRatioTolerance, TDataType AlwaysConvergedNorm) - : ConvergenceCriteria< TSparseSpace, TDenseSpace >() + : BaseType(), + mRatioTolerance(NewRatioTolerance), + mAlwaysConvergedNorm(AlwaysConvergedNorm) { mRatioTolerance = NewRatioTolerance; mAlwaysConvergedNorm = AlwaysConvergedNorm; @@ -103,7 +127,6 @@ class ResidualCriteria //* Copy constructor. explicit ResidualCriteria( ResidualCriteria const& rOther ) :BaseType(rOther) - ,mInitialResidualIsSet(rOther.mInitialResidualIsSet) ,mRatioTolerance(rOther.mRatioTolerance) ,mInitialResidualNorm(rOther.mInitialResidualNorm) ,mCurrentResidualNorm(rOther.mCurrentResidualNorm) @@ -141,14 +164,9 @@ class ResidualCriteria if (size_b != 0) { //if we are solving for something SizeType size_residual; - if (mInitialResidualIsSet == false) { - CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, b); - mInitialResidualIsSet = true; - } - - TDataType ratio = 0.0; CalculateResidualNorm(mCurrentResidualNorm, size_residual, rDofSet, b); + TDataType ratio = 0.0; if(mInitialResidualNorm < std::numeric_limits::epsilon()) { ratio = 0.0; } else { @@ -187,19 +205,23 @@ class ResidualCriteria * @brief This function initializes the solution step * @param rModelPart Reference to the ModelPart containing the problem. * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) - * @param A System matrix (unused) - * @param Dx Vector of results (variations on nodal variables) - * @param b RHS vector (residual + reactions) + * @param rA System matrix (unused) + * @param rDx Vector of results (variations on nodal variables) + * @param rb RHS vector (residual + reactions) */ void InitializeSolutionStep( ModelPart& rModelPart, DofsArrayType& rDofSet, - const TSystemMatrixType& A, - const TSystemVectorType& Dx, - const TSystemVectorType& b + const TSystemMatrixType& rA, + const TSystemVectorType& rDx, + const TSystemVectorType& rb ) override { - mInitialResidualIsSet = false; + BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb); + SizeType size_residual; + CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, rb); + + KRATOS_WATCH(mInitialResidualNorm) } /** @@ -349,9 +371,6 @@ class ResidualCriteria ///@name Member Variables ///@{ - - bool mInitialResidualIsSet; /// This "flag" is set in order to set that the initial residual is already computed - TDataType mRatioTolerance; /// The ratio threshold for the norm of the residual TDataType mInitialResidualNorm; /// The reference norm of the residual diff --git a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h index 3bc6b4a0df3..bc540cd7485 100644 --- a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h +++ b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h @@ -557,13 +557,13 @@ class ResidualBasedNewtonRaphsonStrategy { KRATOS_TRY; - if (mSolutionStepIsInitialized == false) - { - //pointers needed in the solution + if (!mSolutionStepIsInitialized) { + // Pointers needed in the solution typename TSchemeType::Pointer p_scheme = GetScheme(); typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver(); + ModelPart& r_model_part = BaseType::GetModelPart(); - const int rank = BaseType::GetModelPart().GetCommunicator().MyPID(); + const int rank = r_model_part.GetCommunicator().MyPID(); //set up the system, operation performed just once unless it is required //to reform the dof set at each iteration @@ -573,20 +573,20 @@ class ResidualBasedNewtonRaphsonStrategy { //setting up the list of the DOFs to be solved BuiltinTimer setup_dofs_time; - p_builder_and_solver->SetUpDofSet(p_scheme, BaseType::GetModelPart()); + p_builder_and_solver->SetUpDofSet(p_scheme, r_model_part); KRATOS_INFO_IF("Setup Dofs Time", BaseType::GetEchoLevel() > 0 && rank == 0) << setup_dofs_time.ElapsedSeconds() << std::endl; //shaping correctly the system BuiltinTimer setup_system_time; - p_builder_and_solver->SetUpSystem(BaseType::GetModelPart()); + p_builder_and_solver->SetUpSystem(r_model_part); KRATOS_INFO_IF("Setup System Time", BaseType::GetEchoLevel() > 0 && rank == 0) << setup_system_time.ElapsedSeconds() << std::endl; //setting up the Vectors involved to the correct size BuiltinTimer system_matrix_resize_time; p_builder_and_solver->ResizeAndInitializeVectors(p_scheme, mpA, mpDx, mpb, - BaseType::GetModelPart()); + r_model_part); KRATOS_INFO_IF("System Matrix Resize Time", BaseType::GetEchoLevel() > 0 && rank == 0) << system_matrix_resize_time.ElapsedSeconds() << std::endl; } @@ -598,11 +598,17 @@ class ResidualBasedNewtonRaphsonStrategy TSystemVectorType& rDx = *mpDx; TSystemVectorType& rb = *mpb; - //initial operations ... things that are constant over the Solution Step - p_builder_and_solver->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb); + // Initial operations ... things that are constant over the Solution Step + p_builder_and_solver->InitializeSolutionStep(r_model_part, rA, rDx, rb); - //initial operations ... things that are constant over the Solution Step - p_scheme->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb); + // Initial operations ... things that are constant over the Solution Step + p_scheme->InitializeSolutionStep(r_model_part, rA, rDx, rb); + + // Initialisation of the convergence criteria + TSparseSpace::SetToZero(rb); + p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb); + + mpConvergenceCriteria->InitializeSolutionStep(r_model_part, p_builder_and_solver->GetDofSet(), rA, rDx, rb); mSolutionStepIsInitialized = true; } @@ -667,23 +673,19 @@ class ResidualBasedNewtonRaphsonStrategy //initializing the parameters of the Newton-Raphson cycle unsigned int iteration_number = 1; BaseType::GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] = iteration_number; - // BaseType::GetModelPart().GetProcessInfo().SetNonLinearIterationNumber(iteration_number); bool is_converged = false; bool residual_is_updated = false; p_scheme->InitializeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb); is_converged = mpConvergenceCriteria->PreCriteria(BaseType::GetModelPart(), p_builder_and_solver->GetDofSet(), rA, rDx, rb); - //function to perform the building and the solving phase. - if (BaseType::mRebuildLevel > 0 || BaseType::mStiffnessMatrixIsBuilt == false) - { + // Function to perform the building and the solving phase. + if (BaseType::mRebuildLevel > 0 || BaseType::mStiffnessMatrixIsBuilt == false) { TSparseSpace::SetToZero(rA); TSparseSpace::SetToZero(rDx); TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildAndSolve(p_scheme, BaseType::GetModelPart(), rA, rDx, rb); - } - else - { + } else { TSparseSpace::SetToZero(rDx); //Dx=0.00; TSparseSpace::SetToZero(rb); @@ -698,13 +700,8 @@ class ResidualBasedNewtonRaphsonStrategy p_scheme->FinalizeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb); - if (is_converged == true) - { - //initialisation of the convergence criteria - mpConvergenceCriteria->InitializeSolutionStep(BaseType::GetModelPart(), p_builder_and_solver->GetDofSet(), rA, rDx, rb); - - if (mpConvergenceCriteria->GetActualizeRHSflag() == true) - { + if (is_converged) { + if (mpConvergenceCriteria->GetActualizeRHSflag()) { TSparseSpace::SetToZero(rb); p_builder_and_solver->BuildRHS(p_scheme, BaseType::GetModelPart(), rb); From 3225be4f712a1fd232dd84e58d4a53dddff330e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Wed, 6 Mar 2019 10:09:32 +0100 Subject: [PATCH 02/14] conditionally computing RHS and adding correction for MPCs --- .../convergencecriterias/residual_criteria.h | 39 ++++++++++++++++--- .../residualbased_newton_raphson_strategy.h | 7 +++- 2 files changed, 38 insertions(+), 8 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index b3d588ea6a4..b32ba75961f 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -109,6 +109,8 @@ class ResidualCriteria KRATOS_WARNING("ResidualCriteria") << "residual_relative_tolerance or relative_tolerance nor defined on settings. Using default 1.0e-4" << std::endl; mRatioTolerance = 1.0e-4; } + + this->mActualizeRHSIsNeeded = true; } //* Constructor. @@ -119,9 +121,7 @@ class ResidualCriteria mRatioTolerance(NewRatioTolerance), mAlwaysConvergedNorm(AlwaysConvergedNorm) { - mRatioTolerance = NewRatioTolerance; - mAlwaysConvergedNorm = AlwaysConvergedNorm; - mInitialResidualIsSet = false; + this->mActualizeRHSIsNeeded = true; } //* Copy constructor. @@ -133,6 +133,7 @@ class ResidualCriteria ,mAlwaysConvergedNorm(rOther.mAlwaysConvergedNorm) ,mReferenceDispNorm(rOther.mReferenceDispNorm) { + this->mActualizeRHSIsNeeded = true; } //* Destructor. @@ -218,6 +219,31 @@ class ResidualCriteria ) override { BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb); + + mActiveDofs.resize(rDofSet.size()); + + #pragma omp parallel for + for(int i=0; i(mActiveDofs.size()); ++i) + mActiveDofs[i] = true; + + #pragma omp parallel for + for (int i = 0; i < static_cast(rDofSet.size()); i++) { + auto it_dof = rDofSet.begin() + i; + if(it_dof->IsFixed()) + mActiveDofs[it_dof->EquationId()] = false; + } + + for(auto& mpc : rModelPart.MasterSlaveConstraints()) + { + auto& slave_dofs = mpc.GetSlaveDofsVector(); + auto& master_dofs = mpc.GetMasterDofsVector(); + + for(auto& dof : slave_dofs) + mActiveDofs[dof->EquationId()] = false; + for(auto& dof : master_dofs) + mActiveDofs[dof->EquationId()] = false; + } + SizeType size_residual; CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, rb); @@ -329,11 +355,10 @@ class ResidualCriteria for (int i = 0; i < static_cast(rDofSet.size()); i++) { auto it_dof = rDofSet.begin() + i; - IndexType dof_id; + IndexType dof_id = it_dof->EquationId(); TDataType residual_dof_value; - if (it_dof->IsFree()) { - dof_id = it_dof->EquationId(); + if (mActiveDofs[dof_id]) { residual_dof_value = TSparseSpace::GetValue(b,dof_id); residual_solution_norm += std::pow(residual_dof_value, 2); dof_num++; @@ -381,6 +406,8 @@ class ResidualCriteria TDataType mReferenceDispNorm; /// The norm at the beginning of the iterations + std::vector mActiveDofs; + ///@} ///@name Private Operators diff --git a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h index bc540cd7485..071f51b37e7 100644 --- a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h +++ b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h @@ -605,8 +605,11 @@ class ResidualBasedNewtonRaphsonStrategy p_scheme->InitializeSolutionStep(r_model_part, rA, rDx, rb); // Initialisation of the convergence criteria - TSparseSpace::SetToZero(rb); - p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb); + if (mpConvergenceCriteria->GetActualizeRHSflag() == true) + { + TSparseSpace::SetToZero(rb); + p_builder_and_solver->BuildRHS(p_scheme, r_model_part, rb); + } mpConvergenceCriteria->InitializeSolutionStep(r_model_part, p_builder_and_solver->GetDofSet(), rA, rDx, rb); From d8275b373567234491d4f6baf859bafd639056fb Mon Sep 17 00:00:00 2001 From: riccardo Date: Thu, 28 Feb 2019 18:27:06 +0100 Subject: [PATCH 03/14] missing flag --- .../solving_strategies/convergencecriterias/and_criteria.h | 6 ++++++ .../solving_strategies/convergencecriterias/or_criteria.h | 7 +++++++ 2 files changed, 13 insertions(+) diff --git a/kratos/solving_strategies/convergencecriterias/and_criteria.h b/kratos/solving_strategies/convergencecriterias/and_criteria.h index 3fb884a85f8..5b0e13aa77e 100644 --- a/kratos/solving_strategies/convergencecriterias/and_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/and_criteria.h @@ -96,6 +96,9 @@ class And_Criteria mpFirstCriterion(pFirstCriterion), mpSecondCriterion(pSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** @@ -107,6 +110,9 @@ class And_Criteria mpFirstCriterion(rOther.mpFirstCriterion), mpSecondCriterion(rOther.mpSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** Destructor. diff --git a/kratos/solving_strategies/convergencecriterias/or_criteria.h b/kratos/solving_strategies/convergencecriterias/or_criteria.h index fe001d49657..d9d64839d9c 100755 --- a/kratos/solving_strategies/convergencecriterias/or_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/or_criteria.h @@ -95,6 +95,10 @@ class Or_Criteria mpFirstCriterion(pFirstCriterion), mpSecondCriterion(pSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + } /** Copy constructor. @@ -104,6 +108,9 @@ class Or_Criteria mpFirstCriterion(rOther.mpFirstCriterion), mpSecondCriterion(rOther.mpSecondCriterion) { + this->mActualizeRHSIsNeeded = false; + if(mpFirstCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; + if(mpSecondCriterion->GetActualizeRHSflag()) this->mActualizeRHSIsNeeded = true; } /** Destructor. From fe5b282b6361792d941bbd94d62a035a92ad2093 Mon Sep 17 00:00:00 2001 From: Philipp Bucher Date: Thu, 28 Feb 2019 20:37:32 +0100 Subject: [PATCH 04/14] some minor improvements --- .../convergencecriterias/residual_criteria.h | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index b32ba75961f..740a200b5c2 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -223,31 +223,29 @@ class ResidualCriteria mActiveDofs.resize(rDofSet.size()); #pragma omp parallel for - for(int i=0; i(mActiveDofs.size()); ++i) + for(int i=0; i(mActiveDofs.size()); ++i) { mActiveDofs[i] = true; + } #pragma omp parallel for - for (int i = 0; i < static_cast(rDofSet.size()); i++) { + for (int i=0; i(rDofSet.size()); ++i) { auto it_dof = rDofSet.begin() + i; - if(it_dof->IsFixed()) + if (it_dof->IsFixed()) { mActiveDofs[it_dof->EquationId()] = false; + } } - for(auto& mpc : rModelPart.MasterSlaveConstraints()) - { - auto& slave_dofs = mpc.GetSlaveDofsVector(); - auto& master_dofs = mpc.GetMasterDofsVector(); - - for(auto& dof : slave_dofs) - mActiveDofs[dof->EquationId()] = false; - for(auto& dof : master_dofs) - mActiveDofs[dof->EquationId()] = false; + for (const auto& r_mpc : rModelPart.MasterSlaveConstraints()) { + for (const auto& r_dof : r_mpc.GetMasterDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } + for (const auto& r_dof : r_mpc.GetSlaveDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } } SizeType size_residual; CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, rb); - - KRATOS_WATCH(mInitialResidualNorm) } /** From 3c783ed3d05408b733f1514580dd79682446bb8d Mon Sep 17 00:00:00 2001 From: Philipp Bucher Date: Thu, 28 Feb 2019 20:44:45 +0100 Subject: [PATCH 05/14] adding const --- .../convergencecriterias/residual_criteria.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index 740a200b5c2..0a90f923fcd 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -229,7 +229,7 @@ class ResidualCriteria #pragma omp parallel for for (int i=0; i(rDofSet.size()); ++i) { - auto it_dof = rDofSet.begin() + i; + const auto it_dof = rDofSet.begin() + i; if (it_dof->IsFixed()) { mActiveDofs[it_dof->EquationId()] = false; } @@ -353,7 +353,7 @@ class ResidualCriteria for (int i = 0; i < static_cast(rDofSet.size()); i++) { auto it_dof = rDofSet.begin() + i; - IndexType dof_id = it_dof->EquationId(); + const IndexType dof_id = it_dof->EquationId(); TDataType residual_dof_value; if (mActiveDofs[dof_id]) { From ef6a693c26e2473727a5eb57b865cf3157a08df4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Mon, 4 Mar 2019 22:54:00 +0100 Subject: [PATCH 06/14] Computing only when needed --- .../convergencecriterias/residual_criteria.h | 76 ++++++++++++------- 1 file changed, 49 insertions(+), 27 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index 0a90f923fcd..e37dcba5074 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -165,7 +165,7 @@ class ResidualCriteria if (size_b != 0) { //if we are solving for something SizeType size_residual; - CalculateResidualNorm(mCurrentResidualNorm, size_residual, rDofSet, b); + CalculateResidualNorm(rModelPart, mCurrentResidualNorm, size_residual, rDofSet, b); TDataType ratio = 0.0; if(mInitialResidualNorm < std::numeric_limits::epsilon()) { @@ -220,32 +220,35 @@ class ResidualCriteria { BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb); - mActiveDofs.resize(rDofSet.size()); + // Filling mActiveDofs when MPC exist + if (rModelPart.NumberOfMasterSlaveConstraints() > 0) { + mActiveDofs.resize(rDofSet.size()); - #pragma omp parallel for - for(int i=0; i(mActiveDofs.size()); ++i) { - mActiveDofs[i] = true; - } - - #pragma omp parallel for - for (int i=0; i(rDofSet.size()); ++i) { - const auto it_dof = rDofSet.begin() + i; - if (it_dof->IsFixed()) { - mActiveDofs[it_dof->EquationId()] = false; + #pragma omp parallel for + for(int i=0; i(mActiveDofs.size()); ++i) { + mActiveDofs[i] = true; } - } - - for (const auto& r_mpc : rModelPart.MasterSlaveConstraints()) { - for (const auto& r_dof : r_mpc.GetMasterDofsVector()) { - mActiveDofs[r_dof->EquationId()] = false; + + #pragma omp parallel for + for (int i=0; i(rDofSet.size()); ++i) { + const auto it_dof = rDofSet.begin() + i; + if (it_dof->IsFixed()) { + mActiveDofs[it_dof->EquationId()] = false; + } } - for (const auto& r_dof : r_mpc.GetSlaveDofsVector()) { - mActiveDofs[r_dof->EquationId()] = false; + + for (const auto& r_mpc : rModelPart.MasterSlaveConstraints()) { + for (const auto& r_dof : r_mpc.GetMasterDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } + for (const auto& r_dof : r_mpc.GetSlaveDofsVector()) { + mActiveDofs[r_dof->EquationId()] = false; + } } } SizeType size_residual; - CalculateResidualNorm(mInitialResidualNorm, size_residual, rDofSet, rb); + CalculateResidualNorm(rModelPart, mInitialResidualNorm, size_residual, rDofSet, rb); } /** @@ -332,12 +335,14 @@ class ResidualCriteria /** * @brief This method computes the norm of the residual * @details It checks if the dof is fixed + * @param rModelPart Reference to the ModelPart containing the problem. * @param rResidualSolutionNorm The norm of the residual * @param rDofNum The number of DoFs * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) * @param b RHS vector (residual + reactions) */ virtual void CalculateResidualNorm( + ModelPart& rModelPart, TDataType& rResidualSolutionNorm, SizeType& rDofNum, DofsArrayType& rDofSet, @@ -348,15 +353,32 @@ class ResidualCriteria TDataType residual_solution_norm = TDataType(); SizeType dof_num = 0; + // Auxiliar values + TDataType residual_dof_value; + const auto it_dof_begin = rDofSet.begin(); + const int number_of_dof = static_cast(rDofSet.size()); + // Loop over Dofs - #pragma omp parallel for reduction(+:residual_solution_norm,dof_num) - for (int i = 0; i < static_cast(rDofSet.size()); i++) { - auto it_dof = rDofSet.begin() + i; + if (rModelPart.NumberOfMasterSlaveConstraints() > 0) { + #pragma omp parallel for reduction(+:residual_solution_norm, dof_num, residual_dof_value) + for (int i = 0; i < number_of_dof; i++) { + auto it_dof = it_dof_begin + i; + + const IndexType dof_id = it_dof->EquationId(); + + if (mActiveDofs[dof_id]) { + residual_dof_value = TSparseSpace::GetValue(b,dof_id); + residual_solution_norm += std::pow(residual_dof_value, 2); + dof_num++; + } + } + } else { + #pragma omp parallel for reduction(+:residual_solution_norm, dof_num, residual_dof_value) + for (int i = 0; i < number_of_dof; i++) { + auto it_dof = it_dof_begin + i; - const IndexType dof_id = it_dof->EquationId(); - TDataType residual_dof_value; + const IndexType dof_id = it_dof->EquationId(); - if (mActiveDofs[dof_id]) { residual_dof_value = TSparseSpace::GetValue(b,dof_id); residual_solution_norm += std::pow(residual_dof_value, 2); dof_num++; @@ -404,7 +426,7 @@ class ResidualCriteria TDataType mReferenceDispNorm; /// The norm at the beginning of the iterations - std::vector mActiveDofs; + std::vector mActiveDofs; /// This vector contains the dofs that are active ///@} From bb7c720ebf6eca4f421acfbe63038a7f65490c4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Mon, 4 Mar 2019 22:58:29 +0100 Subject: [PATCH 07/14] Update Trilinos version --- .../convergencecriterias/trilinos_residual_criteria.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h b/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h index aace8dc2fdd..657afb44751 100644 --- a/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h +++ b/applications/trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h @@ -82,12 +82,14 @@ class TrilinosResidualCriteria : public ResidualCriteria< TSparseSpace, TDenseSp /** * @brief This method computes the norm of the residual * @details It checks if the dof is fixed + * @param rModelPart Reference to the ModelPart containing the problem. * @param rResidualSolutionNorm The norm of the residual * @param rDofNum The number of DoFs * @param rDofSet Reference to the container of the problem's degrees of freedom (stored by the BuilderAndSolver) * @param b RHS vector (residual + reactions) */ void CalculateResidualNorm( + ModelPart& rModelPart, TDataType& rResidualSolutionNorm, typename BaseType::SizeType& rDofNum, typename BaseType::DofsArrayType& rDofSet, From 3f15728a14c938cb030729cd86722ca034e25409 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Tue, 5 Mar 2019 09:19:49 +0100 Subject: [PATCH 08/14] Compilation error due to firstprivate --- .../convergencecriterias/residual_criteria.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index e37dcba5074..d2a0108f593 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -354,13 +354,13 @@ class ResidualCriteria SizeType dof_num = 0; // Auxiliar values - TDataType residual_dof_value; + TDataType residual_dof_value = 0.0; const auto it_dof_begin = rDofSet.begin(); const int number_of_dof = static_cast(rDofSet.size()); // Loop over Dofs if (rModelPart.NumberOfMasterSlaveConstraints() > 0) { - #pragma omp parallel for reduction(+:residual_solution_norm, dof_num, residual_dof_value) + #pragma omp parallel for firstprivate(residual_dof_value) reduction(+:residual_solution_norm, dof_num) for (int i = 0; i < number_of_dof; i++) { auto it_dof = it_dof_begin + i; @@ -373,7 +373,7 @@ class ResidualCriteria } } } else { - #pragma omp parallel for reduction(+:residual_solution_norm, dof_num, residual_dof_value) + #pragma omp parallel for firstprivate(residual_dof_value) reduction(+:residual_solution_norm, dof_num) for (int i = 0; i < number_of_dof; i++) { auto it_dof = it_dof_begin + i; @@ -428,7 +428,6 @@ class ResidualCriteria std::vector mActiveDofs; /// This vector contains the dofs that are active - ///@} ///@name Private Operators ///@{ From 867c697036e4824aad84effd1ed77282d826fe74 Mon Sep 17 00:00:00 2001 From: riccardo Date: Tue, 5 Mar 2019 09:20:05 +0100 Subject: [PATCH 09/14] resetting the value of b to leave it untouched after convergence criteria initialization --- .../strategies/residualbased_newton_raphson_strategy.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h index 071f51b37e7..34abd6306a3 100644 --- a/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h +++ b/kratos/solving_strategies/strategies/residualbased_newton_raphson_strategy.h @@ -613,6 +613,9 @@ class ResidualBasedNewtonRaphsonStrategy mpConvergenceCriteria->InitializeSolutionStep(r_model_part, p_builder_and_solver->GetDofSet(), rA, rDx, rb); + if (mpConvergenceCriteria->GetActualizeRHSflag() == true) + TSparseSpace::SetToZero(rb); + mSolutionStepIsInitialized = true; } From 653ecce936bab21ef15b1b0cec4ad0d4ffcf85cb Mon Sep 17 00:00:00 2001 From: riccardo Date: Tue, 5 Mar 2019 11:39:00 +0100 Subject: [PATCH 10/14] fix to work correctly in the case of elimination builder and solver --- .../convergencecriterias/residual_criteria.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index d2a0108f593..1ea371b7d69 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -377,10 +377,12 @@ class ResidualCriteria for (int i = 0; i < number_of_dof; i++) { auto it_dof = it_dof_begin + i; - const IndexType dof_id = it_dof->EquationId(); - - residual_dof_value = TSparseSpace::GetValue(b,dof_id); - residual_solution_norm += std::pow(residual_dof_value, 2); + if(it_dof->IsFixed() == false) + { + const IndexType dof_id = it_dof->EquationId(); + residual_dof_value = TSparseSpace::GetValue(b,dof_id); + residual_solution_norm += std::pow(residual_dof_value, 2); + } dof_num++; } } From 68260f95a8f08e029ed314c4d0b653bed3fc0c92 Mon Sep 17 00:00:00 2001 From: Philipp Bucher Date: Tue, 5 Mar 2019 11:46:59 +0100 Subject: [PATCH 11/14] minor --- .../convergencecriterias/residual_criteria.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index 1ea371b7d69..a256b266966 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -377,8 +377,7 @@ class ResidualCriteria for (int i = 0; i < number_of_dof; i++) { auto it_dof = it_dof_begin + i; - if(it_dof->IsFixed() == false) - { + if (!it_dof->IsFixed()) { const IndexType dof_id = it_dof->EquationId(); residual_dof_value = TSparseSpace::GetValue(b,dof_id); residual_solution_norm += std::pow(residual_dof_value, 2); From 00245724e9bb49ab7d7ec4f9f7093e4524f4854b Mon Sep 17 00:00:00 2001 From: Philipp Bucher Date: Tue, 5 Mar 2019 13:03:02 +0100 Subject: [PATCH 12/14] correct dof-number --- .../solving_strategies/convergencecriterias/residual_criteria.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h index a256b266966..196febf0dd9 100644 --- a/kratos/solving_strategies/convergencecriterias/residual_criteria.h +++ b/kratos/solving_strategies/convergencecriterias/residual_criteria.h @@ -381,8 +381,8 @@ class ResidualCriteria const IndexType dof_id = it_dof->EquationId(); residual_dof_value = TSparseSpace::GetValue(b,dof_id); residual_solution_norm += std::pow(residual_dof_value, 2); + dof_num++; } - dof_num++; } } From 7d7c91cffd9113dc7d0181201f747809533a6825 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Tue, 5 Mar 2019 17:21:29 +0100 Subject: [PATCH 13/14] Adjusting tolerance --- .../tests/test_patch_test_truss.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py index 31b00d38ede..29bca3a1080 100644 --- a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py +++ b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py @@ -221,7 +221,7 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): -70650929.0390236,-71205408.69085957,-71758918.27250087,-72311464.28340018, -72863053.1484657,-73413691.21926463,-73963384.77520159,-74512140.02467461, -75059963.10620539,] - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep],6) + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep], delta=1.0e-3) ##node2 node_temp = mp.Nodes[2] @@ -231,7 +231,7 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): #pointLoad self.assertAlmostEqual(load_temp,Force_i) #reaction_x - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1),6) + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1), delta=1.0e-3) #displacement_y EA = 210e9*0.01 L = sqrt(4+1) @@ -282,7 +282,7 @@ def _check_results_cable(self,mp,Force_X): self.assertAlmostEqual(disp_u_2, 0.022296019142446475,6) - self.assertAlmostEqual(r_u_1, -Force_X,6) + self.assertAlmostEqual(r_u_1, -Force_X,delta=1.0e-3) self.assertAlmostEqual(r_u_3, 0.00 ,4) def _check_results_dynamic_explicit(self,mp,time_i,time_step,linear_flag): From e788a52cc6957d6fe0605bd1f58c4a9eaa80e824 Mon Sep 17 00:00:00 2001 From: KlausBSautter Date: Wed, 6 Mar 2019 08:39:28 +0100 Subject: [PATCH 14/14] fixed truss tests --- .../tests/test_patch_test_truss.py | 39 +++---------------- 1 file changed, 5 insertions(+), 34 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py index 29bca3a1080..0d5a78fa2ef 100644 --- a/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py +++ b/applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py @@ -190,38 +190,9 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): reaction_y_node1 = Force_i*(-1) self.assertAlmostEqual(reac_temp[1],reaction_y_node1,6) #reaction_x - reaction_x_node1 = [741464.9276510741,1485888.977636112,2233316.9009164227, - 2983794.615716884,3737369.2509119534,4494089.191516033,5254004.126397622, - 6017165.0983619075,6783624.556737246,7553436.412628534,8326656.0969987875, - 9103340.621764094,9883548.644087134,10667340.53408764,11454778.446184492, - 12245926.394319195,13040850.331311818,13839618.232645638,14642300.18497437, - 15448968.47968232,16259697.711865114,17074564.885111626,17893649.522510014, - 18717033.784337185,19544802.592936598,20377043.765315644,21213848.154068694, - 22055309.79726978,22901526.07803619,23752597.894554257,24608629.841397412, - 25469730.40309285,26336012.160943132,27207592.014241144,28084591.41712631, - 28967136.632448073,29855359.004160944,30749395.24993747,31649387.775854316, - 32555485.015236698,33467841.79396308,34386619.72480394,35311987.633672595, - 36244122.02100001,37183207.5618443,38129437.64879564,39083014.9822317, - 40044152.21308944,41013072.64397469,41990010.99523344,42975214.24348978, - 43968942.54123137,44971470.22722923,45983086.93901924,47004098.840344116, - 48034829.97843962,49075623.78836585,50126844.76436112,51188880.321473464, - 52262142.87466787,53347072.16731208,54444137.88663925,55553842.61067612, - 56676725.13951358,57813364.2740804,58964383.118212394,60130453.99549516, - 61312304.09186774,62510721.95949238,63726565.048343465,64960768.47141639, - 66214355.26005781,67488448.43149039,68784285.27627705,70103234.38657254, - 71446816.09699601,72816727.21376368,74214871.18642414,75643395.26295514, - 77104736.71277626,78601680.98043492,80137435.76669273,81715726.72014935, - 83340922.98755075,85018204.87282823,86753792.28015186,88555263.27634439, - 90432010.47500862,92395916.01811945,94462388.67855237,96652033.38549507, - 98993500.26523624,101528726.98334643,104323616.16359182,107493197.43582612, - 111276440.23647068,116390127.39236663,-62782528.388332605,-63351316.30823133, - -63919034.598836,-64485690.945303164,-65051292.93836311,-65615848.075949684, - -66179363.76479947,-66741847.3220098,-67303305.9765681,-67863746.8708426, - -68423177.06204486,-68981603.5236578,-69539033.1468354,-70095472.74176757, - -70650929.0390236,-71205408.69085957,-71758918.27250087,-72311464.28340018, - -72863053.1484657,-73413691.21926463,-73963384.77520159,-74512140.02467461, - -75059963.10620539,] - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep], delta=1.0e-3) + reaction_x_node1 = [741464.9276515746, 1485888.977636112, + 2233316.9009164227, 2983794.615716549, 3737369.2509122863] + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep], 6) ##node2 node_temp = mp.Nodes[2] @@ -231,7 +202,7 @@ def _check_results_nonlinear(self,mp,timestep,Force_i): #pointLoad self.assertAlmostEqual(load_temp,Force_i) #reaction_x - self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1), delta=1.0e-3) + self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1), 6) #displacement_y EA = 210e9*0.01 L = sqrt(4+1) @@ -282,7 +253,7 @@ def _check_results_cable(self,mp,Force_X): self.assertAlmostEqual(disp_u_2, 0.022296019142446475,6) - self.assertAlmostEqual(r_u_1, -Force_X,delta=1.0e-3) + self.assertAlmostEqual(r_u_1, -100000000.00093779,6) self.assertAlmostEqual(r_u_3, 0.00 ,4) def _check_results_dynamic_explicit(self,mp,time_i,time_step,linear_flag):