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.

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

Add sparse matrix construction tests #324

Merged
merged 4 commits into from
Dec 24, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
19 changes: 19 additions & 0 deletions src/solver/compute_b_num_non_zero.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "compute_number_of_non_zeros_constraints.hpp"

namespace openturbine {
[[nodiscard]] static size_t ComputeBNumNonZero(
const Kokkos::View<ConstraintType*>::const_type& constraint_type,
const Kokkos::View<Kokkos::pair<size_t, size_t>*>::const_type& constraint_row_range
) {
auto B_num_non_zero = size_t{0U};
Kokkos::parallel_reduce(
"ComputeNumberOfNonZeros_Constraints", constraint_type.extent(0),
ComputeNumberOfNonZeros_Constraints{constraint_type, constraint_row_range}, B_num_non_zero
);
return B_num_non_zero;
}
} // namespace openturbine
90 changes: 90 additions & 0 deletions src/solver/compute_k_col_inds.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

template <typename RowPtrType, typename IndicesType>
struct ComputeKColIndsFunction {
using RowPtrValueType = typename RowPtrType::value_type;
using IndicesValueType = typename IndicesType::value_type;
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;
Kokkos::View<size_t*>::const_type node_freedom_map_table;
Kokkos::View<size_t*>::const_type num_nodes_per_element;
Kokkos::View<size_t**>::const_type node_state_indices;
typename RowPtrType::const_type K_row_ptrs;
IndicesType K_col_inds;

KOKKOS_FUNCTION
bool ContainsNode(size_t element, size_t node) const {
bool contains_node = false;
for (auto n = 0U; n < num_nodes_per_element(element); ++n) {
contains_node = contains_node || (node_state_indices(element, n) == node);
}
return contains_node;
}

KOKKOS_FUNCTION
RowPtrValueType ComputeColInds(size_t element, size_t node, RowPtrValueType current_dof_index)
const {
for (auto n = 0U; n < num_nodes_per_element(element); ++n) {
const auto node_state_index = node_state_indices(element, n);
if (node_state_index != node) {
const auto target_node_num_dof =
count_active_dofs(node_freedom_allocation_table(node_state_index));
const auto target_node_dof_index = node_freedom_map_table(node_state_index);
for (auto k = 0U; k < target_node_num_dof; ++k, ++current_dof_index) {
const auto col_index = static_cast<IndicesValueType>(target_node_dof_index + k);
K_col_inds(current_dof_index) = col_index;
}
}
}
return current_dof_index;
}

KOKKOS_FUNCTION
void operator()(size_t i) const {
const auto this_node_num_dof = count_active_dofs(node_freedom_allocation_table(i));
const auto this_node_dof_index = node_freedom_map_table(i);

for (auto j = 0U; j < this_node_num_dof; ++j) {
auto current_dof_index = K_row_ptrs(this_node_dof_index + j);

for (auto k = 0U; k < this_node_num_dof; ++k, ++current_dof_index) {
K_col_inds(current_dof_index) = static_cast<int>(this_node_dof_index + k);
}

for (auto e = 0U; e < num_nodes_per_element.extent(0); ++e) {
if (!ContainsNode(e, i)) {
continue;
}
current_dof_index = ComputeColInds(e, i, current_dof_index);
}
}
}
};

template <typename RowPtrType, typename IndicesType>
[[nodiscard]] inline IndicesType ComputeKColInds(
size_t K_num_non_zero,
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table,
const Kokkos::View<size_t*>::const_type& node_freedom_map_table,
const Kokkos::View<size_t*>::const_type& num_nodes_per_element,
const Kokkos::View<size_t**>::const_type& node_state_indices,
const typename RowPtrType::const_type& K_row_ptrs
) {
auto K_col_inds = IndicesType("col_inds", K_num_non_zero);

Kokkos::parallel_for(
"ComputeKColInds", node_freedom_allocation_table.extent(0),
ComputeKColIndsFunction<RowPtrType, IndicesType>{
node_freedom_allocation_table, node_freedom_map_table, num_nodes_per_element,
node_state_indices, K_row_ptrs, K_col_inds
}
);

return K_col_inds;
}
} // namespace openturbine
66 changes: 66 additions & 0 deletions src/solver/compute_k_num_non_zero.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

struct ComputeKNumNonZero_OffDiagonal {
Kokkos::View<size_t*>::const_type num_nodes_per_element;
Kokkos::View<size_t**>::const_type node_state_indices;
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;

KOKKOS_FUNCTION
void operator()(size_t i, size_t& update) const {
auto num_element_dof = 0UL;
for (auto j = 0U; j < num_nodes_per_element(i); ++j) {
const auto num_node_dof =
count_active_dofs(node_freedom_allocation_table(node_state_indices(i, j)));
num_element_dof += num_node_dof;
}
const auto num_element_non_zero = num_element_dof * num_element_dof;
auto num_diagonal_non_zero = 0UL;
for (auto j = 0U; j < num_nodes_per_element(i); ++j) {
const auto num_node_dof =
count_active_dofs(node_freedom_allocation_table(node_state_indices(i, j)));
num_diagonal_non_zero += num_node_dof * num_node_dof;
}
update += num_element_non_zero - num_diagonal_non_zero;
}
};

struct ComputeKNumNonZero_Diagonal {
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;

KOKKOS_FUNCTION
void operator()(size_t i, size_t& update) const {
const auto num_node_dof = count_active_dofs(node_freedom_allocation_table(i));
const auto num_diagonal_non_zero = num_node_dof * num_node_dof;
update += num_diagonal_non_zero;
}
};

/// Computes the number of non-zero entries in the stiffness matrix K for sparse storage
[[nodiscard]] inline size_t ComputeKNumNonZero(
const Kokkos::View<size_t*>::const_type& num_nodes_per_element,
const Kokkos::View<size_t**>::const_type& node_state_indices,
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table
) {
auto K_num_non_zero_off_diagonal = 0UL;
Kokkos::parallel_reduce(
"ComputeKNumNonZero_OffDiagonal", num_nodes_per_element.extent(0),
ComputeKNumNonZero_OffDiagonal{
num_nodes_per_element, node_state_indices, node_freedom_allocation_table
},
K_num_non_zero_off_diagonal
);
auto K_num_non_zero_diagonal = 0UL;
Kokkos::parallel_reduce(
"ComputeKNumNonZero_Diagonal", node_freedom_allocation_table.extent(0),
ComputeKNumNonZero_Diagonal{node_freedom_allocation_table}, K_num_non_zero_diagonal
);
return K_num_non_zero_off_diagonal + K_num_non_zero_diagonal;
}

} // namespace openturbine
90 changes: 90 additions & 0 deletions src/solver/compute_k_row_ptrs.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

template <typename RowPtrType>
struct ComputeKRowEntries {
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;
Kokkos::View<size_t*>::const_type node_freedom_map_table;
Kokkos::View<size_t*>::const_type num_nodes_per_element;
Kokkos::View<size_t**>::const_type node_state_indices;
RowPtrType K_row_entries;

KOKKOS_FUNCTION
void operator()(size_t i) const {
const auto this_node_num_dof = count_active_dofs(node_freedom_allocation_table(i));
const auto this_node_dof_index = node_freedom_map_table(i);

auto num_entries_in_row = this_node_num_dof;
bool node_found_in_system = false;

// contributions to non-diagonal block from coupled nodes
for (auto e = 0U; e < num_nodes_per_element.extent(0); ++e) {
bool contains_node = false;
auto num_entries_in_element = 0UL;
for (auto j = 0U; j < num_nodes_per_element(e); ++j) {
contains_node = contains_node || (node_state_indices(e, j) == i);
num_entries_in_element +=
count_active_dofs(node_freedom_allocation_table(node_state_indices(e, j)));
}
if (contains_node) {
node_found_in_system = true;
num_entries_in_row += num_entries_in_element - this_node_num_dof;
}
}
if (node_found_in_system) {
for (auto j = 0U; j < this_node_num_dof; ++j) {
K_row_entries(this_node_dof_index + j) = num_entries_in_row;
}
}
}
};

template <typename RowPtrType>
struct ComputeKRowPtrsScanner {
typename RowPtrType::const_type K_row_entries;
RowPtrType K_row_ptrs;

KOKKOS_FUNCTION
void operator()(size_t i, size_t& update, bool is_final) const {
update += K_row_entries(i);
if (is_final) {
K_row_ptrs(i + 1) = update;
}
}
};

/// Computes the row pointers for the sparse stiffness matrix K in CSR format
template <typename RowPtrType>
[[nodiscard]] inline RowPtrType ComputeKRowPtrs(
size_t K_num_rows,
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table,
const Kokkos::View<size_t*>::const_type& node_freedom_map_table,
const Kokkos::View<size_t*>::const_type& num_nodes_per_element,
const Kokkos::View<size_t**>::const_type& node_state_indices
) {
auto K_row_ptrs = RowPtrType("row_ptrs", K_num_rows + 1);
const auto K_row_entries = RowPtrType("row_ptrs", K_num_rows);

Kokkos::parallel_for(
"ComputeKRowEntries", node_freedom_allocation_table.extent(0),
ComputeKRowEntries<RowPtrType>{
node_freedom_allocation_table, node_freedom_map_table, num_nodes_per_element,
node_state_indices, K_row_entries
}
);

auto result = 0UL;
Kokkos::parallel_scan(
"ComputeKRowPtrs", K_row_entries.extent(0),
ComputeKRowPtrsScanner<RowPtrType>{K_row_entries, K_row_ptrs}, result
);

return K_row_ptrs;
}

} // namespace openturbine
30 changes: 30 additions & 0 deletions src/solver/compute_num_system_dofs.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

struct ComputeNumSystemDofsReducer {
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;

KOKKOS_FUNCTION
void operator()(size_t i, size_t& update) const {
update += count_active_dofs(node_freedom_allocation_table(i));
}
};

/// Computes the total number of active degrees of freedom in the system
[[nodiscard]] inline size_t ComputeNumSystemDofs(
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table
) {
auto total_system_dofs = 0UL;
Kokkos::parallel_reduce(
"ComputeNumSystemDofs", node_freedom_allocation_table.extent(0),
ComputeNumSystemDofsReducer{node_freedom_allocation_table}, total_system_dofs
);
return total_system_dofs;
}

} // namespace openturbine
51 changes: 51 additions & 0 deletions src/solver/compute_t_col_inds.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

template <typename RowPtrType, typename IndicesType>
struct ComputeTColIndsFunction {
using IndicesValuesType = typename IndicesType::value_type;
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;
Kokkos::View<size_t*>::const_type node_freedom_map_table;
typename RowPtrType::const_type T_row_ptrs;
IndicesType T_col_inds;

KOKKOS_FUNCTION
void operator()(size_t i) const {
const auto this_node_num_dof = count_active_dofs(node_freedom_allocation_table(i));
const auto this_node_dof_index = node_freedom_map_table(i);

for (auto j = 0U; j < this_node_num_dof; ++j) {
for (auto k = 0UL, current_dof_index = T_row_ptrs(this_node_dof_index + j);
k < this_node_num_dof; ++k, ++current_dof_index) {
const auto col_index = static_cast<IndicesValuesType>(this_node_dof_index + k);
T_col_inds(current_dof_index) = col_index;
}
}
}
};

template <typename RowPtrType, typename IndicesType>
[[nodiscard]] inline IndicesType ComputeTColInds(
size_t T_num_non_zero,
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table,
const Kokkos::View<size_t*>::const_type& node_freedom_map_table,
const typename RowPtrType::const_type& T_row_ptrs
) {
auto T_col_inds = IndicesType("T_indices", T_num_non_zero);

Kokkos::parallel_for(
"ComputeTColInds", node_freedom_allocation_table.extent(0),
ComputeTColIndsFunction<RowPtrType, IndicesType>{
node_freedom_allocation_table, node_freedom_map_table, T_row_ptrs, T_col_inds
}
);

return T_col_inds;
}

} // namespace openturbine
31 changes: 31 additions & 0 deletions src/solver/compute_t_num_non_zero.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
#pragma once

#include <Kokkos_Core.hpp>

#include "src/dof_management/freedom_signature.hpp"

namespace openturbine {

struct ComputeTNumNonZeroReducer {
Kokkos::View<FreedomSignature*>::const_type node_freedom_allocation_table;

KOKKOS_FUNCTION
void operator()(size_t i, size_t& update) const {
const auto num_node_dof = count_active_dofs(node_freedom_allocation_table(i));
const auto num_diagonal_non_zero = num_node_dof * num_node_dof;
update += num_diagonal_non_zero;
}
};

[[nodiscard]] inline size_t ComputeTNumNonZero(
const Kokkos::View<FreedomSignature*>::const_type& node_freedom_allocation_table
) {
auto T_num_non_zero = 0UL;
Kokkos::parallel_reduce(
"ComputeTNumNonZero", node_freedom_allocation_table.extent(0),
ComputeTNumNonZeroReducer{node_freedom_allocation_table}, T_num_non_zero
);

return T_num_non_zero;
}
} // namespace openturbine
Loading
Loading