Skip to content

[Core] Add parallel try catch to remaining raw omp loops in builder and solver #14319

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rfaasse wants to merge 4 commits into
base: master
Choose a base branch
from
+74 −26
Open

[Core] Add parallel try catch to remaining raw omp loops in builder and solver #14319

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6325b8a
Added the parallel try-catch macros to the remaining element/conditio…
rfaasse Mar 19, 2026
44bf0ff
Added the parallel try-catch macros to the remaining loops
rfaasse Mar 25, 2026
f963785
Added thread number as an argument to the catch macro
rfaasse Mar 30, 2026
e901edc
Made the omp_get_thread_num call conditional using a preprocessor dir…
rfaasse Mar 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 62 additions & 14 deletions kratos/solving_strategies/builder_and_solvers/residualbased_block_builder_and_solver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,11 +223,15 @@ class ResidualBasedBlockBuilderAndSolver
// Assemble all elements
const auto timer = BuiltinTimer();

KRATOS_PREPARE_CATCH_THREAD_EXCEPTION

#pragma omp parallel firstprivate(nelements,nconditions, LHS_Contribution, RHS_Contribution, EquationId )
{
# pragma omp for schedule(guided, 512) nowait
for (int k = 0; k < nelements; k++) {
auto it_elem = el_begin + k;
for (int i = 0; i < nelements; i++) {
KRATOS_TRY

auto it_elem = el_begin + i;

if (it_elem->IsActive()) {
// Calculate elemental contribution
Expand All @@ -237,11 +241,14 @@ class ResidualBasedBlockBuilderAndSolver
Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
}

KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}

#pragma omp for schedule(guided, 512)
for (int k = 0; k < nconditions; k++) {
auto it_cond = cond_begin + k;
for (int i = 0; i < nconditions; i++) {
KRATOS_TRY

auto it_cond = cond_begin + i;

if (it_cond->IsActive()) {
// Calculate elemental contribution
Expand All @@ -250,9 +257,13 @@ class ResidualBasedBlockBuilderAndSolver
// Assemble the elemental contribution
Assemble(A, b, LHS_Contribution, RHS_Contribution, EquationId);
}
KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())

}
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION

KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", this->GetEchoLevel() >= 1) << "Build time: " << timer << std::endl;

KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", this->GetEchoLevel() > 2) << "Finished parallel building" << std::endl;
Expand Down Expand Up @@ -297,11 +308,15 @@ class ResidualBasedBlockBuilderAndSolver
// Assemble all elements
const auto timer = BuiltinTimer();

KRATOS_PREPARE_CATCH_THREAD_EXCEPTION

#pragma omp parallel firstprivate(nelements, nconditions, lhs_contribution, equation_id )
{
# pragma omp for schedule(guided, 512) nowait
for (int k = 0; k < nelements; ++k) {
auto it_elem = it_elem_begin + k;
for (int i = 0; i < nelements; ++i) {
KRATOS_TRY

auto it_elem = it_elem_begin + i;

// Detect if the element is active or not. If the user did not make any choice the element is active by default
if (it_elem->IsActive()) {
Expand All @@ -311,11 +326,15 @@ class ResidualBasedBlockBuilderAndSolver
// Assemble the elemental contribution
AssembleLHS(rA, lhs_contribution, equation_id);
}

KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}

#pragma omp for schedule(guided, 512)
for (int k = 0; k < nconditions; ++k) {
auto it_cond = it_cond_begin + k;
for (int i = 0; i < nconditions; ++i) {
KRATOS_TRY

auto it_cond = it_cond_begin + i;

// Detect if the element is active or not. If the user did not make any choice the element is active by default
if (it_cond->IsActive()) {
Expand All @@ -325,9 +344,13 @@ class ResidualBasedBlockBuilderAndSolver
// Assemble the elemental contribution
AssembleLHS(rA, lhs_contribution, equation_id);
}

KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION

KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", this->GetEchoLevel() >= 1) << "Build time LHS: " << timer << std::endl;

KRATOS_INFO_IF("ResidualBasedBlockBuilderAndSolver", this->GetEchoLevel() > 2) << "Finished parallel building LHS" << std::endl;
Expand Down Expand Up @@ -1195,7 +1218,7 @@ class ResidualBasedBlockBuilderAndSolver
AssembleRHS(b, RHS_Contribution, EquationId);
}

KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}

LHS_Contribution.resize(0, 0, false);
Expand All @@ -1217,7 +1240,7 @@ class ResidualBasedBlockBuilderAndSolver
AssembleRHS(b, RHS_Contribution, EquationId);
}

KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}
}

Expand All @@ -1239,21 +1262,27 @@ class ResidualBasedBlockBuilderAndSolver

std::vector<LockObject> lock_array(indices.size());

KRATOS_PREPARE_CATCH_THREAD_EXCEPTION

#pragma omp parallel
{
Element::EquationIdVectorType slave_ids(3);
Element::EquationIdVectorType master_ids(3);
std::unordered_map<IndexType, std::unordered_set<IndexType>> temp_indices;

#pragma omp for schedule(guided, 512) nowait
for (int i_const = 0; i_const < static_cast<int>(rModelPart.MasterSlaveConstraints().size()); ++i_const) {
auto it_const = it_const_begin + i_const;
for (int i = 0; i < static_cast<int>(rModelPart.MasterSlaveConstraints().size()); ++i) {
KRATOS_TRY

auto it_const = it_const_begin + i;
it_const->EquationIdVector(slave_ids, master_ids, r_current_process_info);

// Slave DoFs
for (auto &id_i : slave_ids) {
temp_indices[id_i].insert(master_ids.begin(), master_ids.end());
}

KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}

// Merging all the temporal indexes
Expand All @@ -1264,6 +1293,8 @@ class ResidualBasedBlockBuilderAndSolver
}
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION

mSlaveIds.clear();
mMasterIds.clear();
for (int i = 0; i < static_cast<int>(indices.size()); ++i) {
Expand Down Expand Up @@ -1337,13 +1368,17 @@ class ResidualBasedBlockBuilderAndSolver
// We clear the set
mInactiveSlaveDofs.clear();

KRATOS_PREPARE_CATCH_THREAD_EXCEPTION

#pragma omp parallel firstprivate(transformation_matrix, constant_vector, slave_equation_ids, master_equation_ids)
{
std::unordered_set<IndexType> auxiliar_inactive_slave_dofs;

#pragma omp for schedule(guided, 512)
for (int i_const = 0; i_const < number_of_constraints; ++i_const) {
auto it_const = rModelPart.MasterSlaveConstraints().begin() + i_const;
for (int i = 0; i < number_of_constraints; ++i) {
KRATOS_TRY

auto it_const = rModelPart.MasterSlaveConstraints().begin() + i;
it_const->EquationIdVector(slave_equation_ids, master_equation_ids, r_current_process_info);

// If the constraint is active
Expand All @@ -1364,6 +1399,8 @@ class ResidualBasedBlockBuilderAndSolver
} else { // Taking into account inactive constraints
auxiliar_inactive_slave_dofs.insert(slave_equation_ids.begin(), slave_equation_ids.end());
}

KRATOS_CATCH_THREAD_EXCEPTION(GetCurrentThreadId())
}

// We merge all the sets in one thread
Expand All @@ -1373,6 +1410,8 @@ class ResidualBasedBlockBuilderAndSolver
}
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION

// Setting the master dofs into the T and C system
for (auto eq_id : mMasterIds) {
mConstantVector[eq_id] = 0.0;
Expand Down Expand Up @@ -1705,6 +1744,15 @@ class ResidualBasedBlockBuilderAndSolver
return pos;
}

int GetCurrentThreadId() const
{
#ifdef KRATOS_SMP_OPENMP
return omp_get_thread_num();
#else
return 0;
#endif
}

///@}
///@name Private Operations
///@{
Expand Down
24 changes: 12 additions & 12 deletions kratos/utilities/parallel_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,16 +41,16 @@
#define KRATOS_PREPARE_CATCH_THREAD_EXCEPTION std::stringstream err_stream;

#ifndef KRATOS_NO_TRY_CATCH
#define KRATOS_CATCH_THREAD_EXCEPTION \
#define KRATOS_CATCH_THREAD_EXCEPTION(ThreadNumber) \
} catch(Exception& e) { \
KRATOS_CRITICAL_SECTION \
err_stream << "Thread #" << i << " caught exception: " << e.what(); \
err_stream << "Thread #" << ThreadNumber << " caught exception: " << e.what(); \
} catch(std::exception& e) { \
KRATOS_CRITICAL_SECTION \
err_stream << "Thread #" << i << " caught exception: " << e.what(); \
err_stream << "Thread #" << ThreadNumber << " caught exception: " << e.what(); \
} catch(...) { \
KRATOS_CRITICAL_SECTION \
err_stream << "Thread #" << i << " caught unknown exception:"; \
err_stream << "Thread #" << ThreadNumber << " caught unknown exception:"; \
}
#else
#define KRATOS_CATCH_THREAD_EXCEPTION {}
Expand Down Expand Up @@ -189,7 +189,7 @@ class BlockPartition
for (auto it = mBlockPartition[i]; it != mBlockPartition[i+1]; ++it) {
f(*it); //note that we pass the value to the function, not the iterator
}
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand All @@ -214,7 +214,7 @@ class BlockPartition
local_reducer.LocalReduce(f(*it));
}
global_reducer.ThreadSafeReduce(local_reducer);
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}

KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand Down Expand Up @@ -244,7 +244,7 @@ class BlockPartition
for (auto it = mBlockPartition[i]; it != mBlockPartition[i+1]; ++it){
f(*it, thread_local_storage); // note that we pass the value to the function, not the iterator
}
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand Down Expand Up @@ -278,7 +278,7 @@ class BlockPartition
local_reducer.LocalReduce(f(*it, thread_local_storage));
}
global_reducer.ThreadSafeReduce(local_reducer);
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand Down Expand Up @@ -525,7 +525,7 @@ class IndexPartition
for (auto k = mBlockPartition[i]; k < mBlockPartition[i+1]; ++k) {
f(k); //note that we pass a reference to the value, not the iterator
}
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
}
Expand All @@ -549,7 +549,7 @@ class IndexPartition
local_reducer.LocalReduce(f(k));
}
global_reducer.ThreadSafeReduce(local_reducer);
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
return global_reducer.GetValue();
Expand Down Expand Up @@ -578,7 +578,7 @@ class IndexPartition
for (auto k = mBlockPartition[i]; k < mBlockPartition[i+1]; ++k) {
f(k, thread_local_storage); //note that we pass a reference to the value, not the iterator
}
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand Down Expand Up @@ -612,7 +612,7 @@ class IndexPartition
local_reducer.LocalReduce(f(k, thread_local_storage));
}
global_reducer.ThreadSafeReduce(local_reducer);
KRATOS_CATCH_THREAD_EXCEPTION
KRATOS_CATCH_THREAD_EXCEPTION(i)
}
}
KRATOS_CHECK_AND_THROW_THREAD_EXCEPTION
Expand Down
Loading