Skip to content

Commit

Permalink
Merge pull request #4350 from KratosMultiphysics/release/residual-cit…
Browse files Browse the repository at this point in the history 
…eria-NR-fix

[Core][Release] Residual criteria and block B&S correction
  • Loading branch information
@roigcarlo
roigcarlo authored Mar 6, 2019
2 parents 646f2f2 + e788a52 commit 64271eb
Show file tree
Hide file tree
Showing 6 changed files with 147 additions and 92 deletions.
39 changes: 5 additions & 34 deletions applications/StructuralMechanicsApplication/tests/test_patch_test_truss.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],6)
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]
Expand All @@ -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),6)
self.assertAlmostEqual(reac_temp[0],reaction_x_node1[timestep]*(-1), 6)
#displacement_y
EA = 210e9*0.01
L = sqrt(4+1)
Expand Down Expand Up @@ -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,6)
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):
Expand Down
2 changes: 2 additions & 0 deletions .../trilinos_application/custom_strategies/convergencecriterias/trilinos_residual_criteria.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
6 changes: 6 additions & 0 deletions kratos/solving_strategies/convergencecriterias/and_criteria.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

/**
Expand All @@ -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.
Expand Down
7 changes: 7 additions & 0 deletions kratos/solving_strategies/convergencecriterias/or_criteria.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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.
Expand Down
132 changes: 99 additions & 33 deletions kratos/solving_strategies/convergencecriterias/residual_criteria.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,27 +89,51 @@ 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;
}

this->mActualizeRHSIsNeeded = true;
}

//* Constructor.
explicit ResidualCriteria(
TDataType NewRatioTolerance,
TDataType AlwaysConvergedNorm)
: ConvergenceCriteria< TSparseSpace, TDenseSpace >()
: BaseType(),
mRatioTolerance(NewRatioTolerance),
mAlwaysConvergedNorm(AlwaysConvergedNorm)
{
mRatioTolerance = NewRatioTolerance;
mAlwaysConvergedNorm = AlwaysConvergedNorm;
mInitialResidualIsSet = false;
this->mActualizeRHSIsNeeded = true;
}

//* Copy constructor.
explicit ResidualCriteria( ResidualCriteria const& rOther )
:BaseType(rOther)
,mInitialResidualIsSet(rOther.mInitialResidualIsSet)
,mRatioTolerance(rOther.mRatioTolerance)
,mInitialResidualNorm(rOther.mInitialResidualNorm)
,mCurrentResidualNorm(rOther.mCurrentResidualNorm)
,mAlwaysConvergedNorm(rOther.mAlwaysConvergedNorm)
,mReferenceDispNorm(rOther.mReferenceDispNorm)
{
this->mActualizeRHSIsNeeded = true;
}

//* Destructor.
Expand Down Expand Up @@ -141,14 +165,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;
}
CalculateResidualNorm(rModelPart, mCurrentResidualNorm, size_residual, rDofSet, b);

TDataType ratio = 0.0;
CalculateResidualNorm(mCurrentResidualNorm, size_residual, rDofSet, b);

if(mInitialResidualNorm < std::numeric_limits<TDataType>::epsilon()) {
ratio = 0.0;
} else {
Expand Down Expand Up @@ -187,19 +206,49 @@ 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);

// Filling mActiveDofs when MPC exist
if (rModelPart.NumberOfMasterSlaveConstraints() > 0) {
mActiveDofs.resize(rDofSet.size());

#pragma omp parallel for
for(int i=0; i<static_cast<int>(mActiveDofs.size()); ++i) {
mActiveDofs[i] = true;
}

#pragma omp parallel for
for (int i=0; i<static_cast<int>(rDofSet.size()); ++i) {
const auto it_dof = rDofSet.begin() + i;
if (it_dof->IsFixed()) {
mActiveDofs[it_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(rModelPart, mInitialResidualNorm, size_residual, rDofSet, rb);
}

/**
Expand Down Expand Up @@ -286,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,
Expand All @@ -302,19 +353,36 @@ class ResidualCriteria
TDataType residual_solution_norm = TDataType();
SizeType dof_num = 0;

// Auxiliar values
TDataType residual_dof_value = 0.0;
const auto it_dof_begin = rDofSet.begin();
const int number_of_dof = static_cast<int>(rDofSet.size());

// Loop over Dofs
#pragma omp parallel for reduction(+:residual_solution_norm,dof_num)
for (int i = 0; i < static_cast<int>(rDofSet.size()); i++) {
auto it_dof = rDofSet.begin() + i;

IndexType dof_id;
TDataType residual_dof_value;

if (it_dof->IsFree()) {
dof_id = it_dof->EquationId();
residual_dof_value = TSparseSpace::GetValue(b,dof_id);
residual_solution_norm += std::pow(residual_dof_value, 2);
dof_num++;
if (rModelPart.NumberOfMasterSlaveConstraints() > 0) {
#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;

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 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;

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

Expand Down Expand Up @@ -349,9 +417,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
Expand All @@ -362,6 +427,7 @@ class ResidualCriteria

TDataType mReferenceDispNorm; /// The norm at the beginning of the iterations

std::vector<bool> mActiveDofs; /// This vector contains the dofs that are active

///@}
///@name Private Operators
Expand Down
Loading

0 comments on commit 64271eb

Please sign in to comment.