Using paralell utilities in DisplacementCriteria and `ResidualCrite…

…ria`
KratosMultiphysics · Nov 17, 2023 · dfc099d · dfc099d
1 parent 873813e
commit dfc099d
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 48 deletions.
diff --git a/kratos/solving_strategies/convergencecriterias/displacement_criteria.h b/kratos/solving_strategies/convergencecriterias/displacement_criteria.h
@@ -19,7 +19,8 @@
 
 // Project includes
 #include "includes/model_part.h"
-#include "includes/define.h"
+#include "utilities/parallel_utilities.h"
+#include "utilities/reduction_utilities.h"
 #include "solving_strategies/convergencecriterias/convergence_criteria.h"
 
 namespace Kratos
@@ -385,18 +386,17 @@ class DisplacementCriteria
      */
     void CalculateReferenceNorm(DofsArrayType& rDofSet)
     {
-        TDataType reference_disp_norm = TDataType();
-        TDataType dof_value;
+        // Auxiliary struct
+        struct TLS {TDataType dof_value{};};
 
-        #pragma omp parallel for reduction(+:reference_disp_norm)
-        for (int i = 0; i < static_cast<int>(rDofSet.size()); i++) {
-            auto it_dof = rDofSet.begin() + i;
-
-            if(it_dof->IsFree()) {
-                dof_value = it_dof->GetSolutionStepValue();
-                reference_disp_norm += dof_value * dof_value;
+        const TDataType reference_disp_norm = block_for_each<SumReduction<TDataType>>(rDofSet, TLS(), [](auto& rDof, TLS& rTLS) {
+            if(rDof.IsFree()) {
+                rTLS.dof_value = rDof.GetSolutionStepValue();
+                return std::pow(rTLS.dof_value, 2);
+            } else {
+                return TDataType();
             }
-        }
+        });
         mReferenceDispNorm = std::sqrt(reference_disp_norm);
     }
 
@@ -417,21 +417,22 @@ class DisplacementCriteria
         TDataType final_correction_norm = TDataType();
         SizeType dof_num = 0;
 
+        // Custom reduction
+        using CustomReduction = CombinedReduction<SumReduction<TDataType>,SumReduction<SizeType>>;
+
+        // Auxiliary struct
+        struct TLS {TDataType dof_value{};IndexType dof_id{}; TDataType variation_dof_value{};};
+
         // Loop over Dofs
-        #pragma omp parallel for reduction(+:final_correction_norm,dof_num)
-        for (int i = 0; i < static_cast<int>(rDofSet.size()); i++) {
-            auto it_dof = rDofSet.begin() + i;
-
-            IndexType dof_id;
-            TDataType variation_dof_value;
-
-            if (it_dof->IsFree()) {
-                dof_id = it_dof->EquationId();
-                variation_dof_value = Dx[dof_id];
-                final_correction_norm += std::pow(variation_dof_value, 2);
-                dof_num++;
+        std::tie(final_correction_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [&Dx](auto& rDof, TLS& rTLS) {
+            if (rDof.IsFree()) {
+                rTLS.dof_id = rDof.EquationId();
+                rTLS.variation_dof_value = Dx[rTLS.dof_id];
+                return std::make_tuple(std::pow(rTLS.variation_dof_value, 2), 1);
+            } else {
+                return std::make_tuple(0.0, 0);
             }
-        }
+        });
 
         rDofNum = dof_num;
         return std::sqrt(final_correction_norm);

diff --git a/kratos/solving_strategies/convergencecriterias/residual_criteria.h b/kratos/solving_strategies/convergencecriterias/residual_criteria.h
@@ -20,6 +20,8 @@
 #include "includes/model_part.h"
 #include "includes/define.h"
 #include "utilities/constraint_utilities.h"
+#include "utilities/parallel_utilities.h"
+#include "utilities/reduction_utilities.h"
 #include "solving_strategies/convergencecriterias/convergence_criteria.h"
 
 namespace Kratos
@@ -360,37 +362,33 @@ 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<int>(rDofSet.size());
+        // Custom reduction
+        using CustomReduction = CombinedReduction<SumReduction<TDataType>,SumReduction<SizeType>>;
+
+        // Auxiliary struct
+        struct TLS {TDataType residual_dof_value{};};
 
         // Loop over Dofs
         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();
-
+            std::tie(residual_solution_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [this, &rb](auto& rDof, TLS& rTLS) {
+                const IndexType dof_id = rDof.EquationId();
                 if (mActiveDofs[dof_id] == 1) {
-                    residual_dof_value = TSparseSpace::GetValue(rb,dof_id);
-                    residual_solution_norm += std::pow(residual_dof_value, 2);
-                    dof_num++;
+                    rTLS.residual_dof_value = TSparseSpace::GetValue(rb, dof_id);
+                    return std::make_tuple(std::pow(rTLS.residual_dof_value, 2), 1);
+                } else {
+                    return std::make_tuple(0.0, 0);
                 }
-            }
+            });
         } 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(rb,dof_id);
-                    residual_solution_norm += std::pow(residual_dof_value, 2);
-                    dof_num++;
+            std::tie(residual_solution_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [&rb](auto& rDof, TLS& rTLS) {
+                if (!rDof.IsFixed()) {
+                    const IndexType dof_id = rDof.EquationId();
+                    rTLS.residual_dof_value = TSparseSpace::GetValue(rb, dof_id);
+                    return std::make_tuple(std::pow(rTLS.residual_dof_value, 2), 1);
+                } else {
+                    return std::make_tuple(0.0, 0);
                 }
-            }
+            });
         }
 
         rDofNum = dof_num;