Skip to content
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

[Core] Adding debug check in CheckAndCorrectZeroDiagonalValues #12782

Merged
merged 10 commits into from
Nov 21, 2024
Merged

[Core] Adding debug check in CheckAndCorrectZeroDiagonalValues #12782

Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
b5422c6
[Core] Adding debug check in `CheckAndCorrectZeroDiagonalValues`
loumalouomega Oct 24, 2024
2982990
@riccardorossi suggestion refactor
loumalouomega Oct 24, 2024
81c57e4
Refactor debug check in `ublas_space.h`
loumalouomega Oct 24, 2024
0bfd872
Refactor debug check in ublas_space.h with proper iterators
loumalouomega Oct 25, 2024
0b99dd3
Iterators do not require &
loumalouomega Oct 25, 2024
52805a8
Iterators do not require &
loumalouomega Oct 25, 2024
ea6950b
Minimum a warning in release
loumalouomega Oct 25, 2024
374f064
@matekelemen suggestion
loumalouomega Oct 25, 2024
da287c4
Simplification
loumalouomega Oct 25, 2024
b89c341
Merge master
loumalouomega Nov 21, 2024
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
29 changes: 23 additions & 6 deletions kratos/spaces/ublas_space.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -708,10 +708,16 @@ class UblasSpace
const SCALING_DIAGONAL ScalingDiagonal = SCALING_DIAGONAL::NO_SCALING
)
{
// The system size
const std::size_t system_size = rA.size1();

const auto& Avalues = rA.value_data();
const auto& Arow_indices = rA.index1_data();
// The matrix data
auto& r_Avalues = rA.value_data();
const auto& r_Arow_indices = rA.index1_data();
const auto& r_Acol_indices = rA.index2_data();

// Define the iterators
const auto it_Acol_indices_begin = r_Acol_indices.begin();

// Define zero value tolerance
const double zero_tolerance = std::numeric_limits<double>::epsilon();
Expand All @@ -723,18 +729,29 @@ class UblasSpace
IndexPartition(system_size).for_each([&](std::size_t Index){
bool empty = true;

const std::size_t col_begin = Arow_indices[Index];
const std::size_t col_end = Arow_indices[Index + 1];
const std::size_t col_begin = r_Arow_indices[Index];
const std::size_t col_end = r_Arow_indices[Index + 1];

for (std::size_t j = col_begin; j < col_end; ++j) {
if(std::abs(Avalues[j]) > zero_tolerance) {
if(std::abs(r_Avalues[j]) > zero_tolerance) {
empty = false;
break;
}
}

if(empty) {
rA(Index, Index) = scale_factor;
const auto it_Acol_indices_row_begin = it_Acol_indices_begin + col_begin;
const auto it_Acol_indices_row_end = it_Acol_indices_begin + col_end;

const auto lower = std::lower_bound(it_Acol_indices_row_begin, it_Acol_indices_row_end, Index);
const auto upper = std::upper_bound(it_Acol_indices_row_begin, it_Acol_indices_row_end, Index);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const auto upper = std::upper_bound(it_Acol_indices_row_begin, it_Acol_indices_row_end, Index);
const auto upper = std::upper_bound(lower, it_Acol_indices_row_end, Index);

sry for the pickiness, but if we're doing it, we may as well do it properly.

No need to start searching from the range begin again. lower_bound <= upper_bound always

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i don't think you need the upper (in line with the suggestion of @matekelemen):

if(empty) {
                const auto it_Acol_indices_row_begin = it_Acol_indices_begin + col_begin;
                const auto it_Acol_indices_row_end = it_Acol_indices_begin + col_end;

                const auto lower = std::lower_bound(it_Acol_indices_row_begin, it_Acol_indices_row_end, Index);
                if (lower != it_Acol_indices_row_end && *lower==j) { // Index was found
                    r_Avalues[std::distance(it_Acol_indices_begin, lower)] = scale_factor;
                } else {
            #ifdef KRATOS_DEBUG
                    KRATOS_ERROR << "Diagonal term (" << Index << ", " << Index << ") is not defined in the system matrix" << std::endl;
            #else

per regards the error i am not really sure.

This s a scaling operation. If the row is not there than i would say it makes no sense to scale it ... not that this is an error ...
What is actually wrong in my opinion is to put a nonzero diagonal where there was no diagonal (in a sense equivalent to having a zero diagonal)


if (lower != upper) { // Index was found
r_Avalues[std::distance(it_Acol_indices_begin, lower)] = scale_factor;
} else {
KRATOS_DEBUG_ERROR << "Diagonal term (" << Index << ", " << Index << ") is not defined in the system matrix" << std::endl;
KRATOS_WARNING("UblasSpace") << "Diagonal term (" << Index << ", " << Index << ") is not defined in the system matrix" << std::endl;
}
loumalouomega marked this conversation as resolved.
Show resolved Hide resolved
rb[Index] = 0.0;
}
});
Expand Down
Loading