Split functions within the step folder into .hpp and .cpp files

src/constraints/constraints.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <array>
+#include <numeric>
 #include <vector>
 
 #include <Kokkos_Core.hpp>

src/solver/solver.hpp

@@ -88,7 +88,7 @@ struct Solver {
         const Kokkos::View<Kokkos::pair<size_t, size_t>*>::const_type& constraint_row_range
     )
         : num_system_nodes(node_IDs.extent(0)),
-          num_system_dofs(num_system_nodes * kLieAlgebraComponents),
+          num_system_dofs(num_system_nodes * 6U),
           num_dofs(num_system_dofs + num_constraint_dofs),
           R("R", num_dofs),
           x("x", num_dofs) {

src/step/CMakeLists.txt

@@ -1 +1,21 @@
-target_sources(openturbine_library PRIVATE)
+target_sources(openturbine_library PRIVATE 
+  assemble_constraints_matrix.cpp
+  assemble_constraints_residual.cpp
+  assemble_inertia_matrix.cpp
+  assemble_residual_vector.cpp
+  assemble_stiffness_matrix.cpp
+  assemble_system_matrix.cpp
+  assemble_system_residual.cpp
+  assemble_tangent_operator.cpp
+  calculate_convergence_error.cpp
+  post_step_clean_up.cpp
+  predict_next_state.cpp
+  reset_constraints.cpp
+  solve_system.cpp
+  step.cpp
+  update_constraint_prediction.cpp
+  update_constraint_variables.cpp
+  update_state_prediction.cpp
+  update_system_variables.cpp
+  update_tangent_operator.cpp
+)

src/step/assemble_constraints_matrix.cpp

@@ -0,0 +1,55 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_Profiling_ScopedRegion.hpp>
+
+#include "src/constraints/constraints.hpp"
+#include "src/solver/copy_constraints_to_sparse_matrix.hpp"
+#include "src/solver/copy_sparse_values_to_transpose.hpp"
+#include "src/solver/solver.hpp"
+
+namespace openturbine {
+void AssembleConstraintsMatrix(Solver& solver, const Constraints& constraints) {
+    auto region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Matrix");
+
+    if (constraints.num == 0) {
+        return;
+    }
+
+    {
+        auto const_region = Kokkos::Profiling::ScopedRegion("Constraints Matrix");
+        auto constraint_policy =
+            Kokkos::TeamPolicy<>(static_cast<int>(constraints.num), Kokkos::AUTO());
+
+        Kokkos::parallel_for(
+            "CopyConstraintsToSparseMatrix", constraint_policy,
+            CopyConstraintsToSparseMatrix<Solver::CrsMatrixType>{
+                constraints.row_range, constraints.base_node_col_range,
+                constraints.target_node_col_range, solver.B, constraints.base_gradient_terms,
+                constraints.target_gradient_terms}
+        );
+    }
+
+    {
+        auto trans_region = Kokkos::Profiling::ScopedRegion("Transpose Constraints Matrix");
+        auto B_num_rows = solver.B.numRows();
+        auto constraint_transpose_policy = Kokkos::TeamPolicy<>(B_num_rows, Kokkos::AUTO());
+        auto tmp_row_map = Solver::RowPtrType(
+            Kokkos::view_alloc(Kokkos::WithoutInitializing, "tmp_row_map"),
+            solver.B_t.graph.row_map.extent(0)
+        );
+        Kokkos::deep_copy(tmp_row_map, solver.B_t.graph.row_map);
+        Kokkos::parallel_for(
+            "CopySparseValuesToTranspose", constraint_transpose_policy,
+            CopySparseValuesToTranspose<Solver::CrsMatrixType>{solver.B, tmp_row_map, solver.B_t}
+        );
+        KokkosSparse::sort_crs_matrix(solver.B_t);
+    }
+
+    {
+        auto mult_region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Matrix");
+        KokkosSparse::spgemm_numeric(
+            solver.constraints_spgemm_handle, solver.B, false, solver.T, false,
+            solver.constraints_matrix
+        );
+    }
+}
+}  // namespace openturbine

src/step/assemble_constraints_matrix.hpp

@@ -1,57 +1,8 @@
 #pragma once
 
-#include <Kokkos_Core.hpp>
-#include <Kokkos_Profiling_ScopedRegion.hpp>
-
-#include "src/constraints/constraints.hpp"
-#include "src/solver/copy_constraints_to_sparse_matrix.hpp"
-#include "src/solver/copy_sparse_values_to_transpose.hpp"
-#include "src/solver/solver.hpp"
-
 namespace openturbine {
-inline void AssembleConstraintsMatrix(Solver& solver, Constraints& constraints) {
-    auto region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Matrix");
-
-    if (constraints.num == 0) {
-        return;
-    }
-
-    {
-        auto const_region = Kokkos::Profiling::ScopedRegion("Constraints Matrix");
-        auto constraint_policy =
-            Kokkos::TeamPolicy<>(static_cast<int>(constraints.num), Kokkos::AUTO());
-
-        Kokkos::parallel_for(
-            "CopyConstraintsToSparseMatrix", constraint_policy,
-            CopyConstraintsToSparseMatrix<Solver::CrsMatrixType>{
-                constraints.row_range, constraints.base_node_col_range,
-                constraints.target_node_col_range, solver.B, constraints.base_gradient_terms,
-                constraints.target_gradient_terms}
-        );
-    }
-
-    {
-        auto trans_region = Kokkos::Profiling::ScopedRegion("Transpose Constraints Matrix");
-        auto B_num_rows = solver.B.numRows();
-        auto constraint_transpose_policy = Kokkos::TeamPolicy<>(B_num_rows, Kokkos::AUTO());
-        auto tmp_row_map = Solver::RowPtrType(
-            Kokkos::view_alloc(Kokkos::WithoutInitializing, "tmp_row_map"),
-            solver.B_t.graph.row_map.extent(0)
-        );
-        Kokkos::deep_copy(tmp_row_map, solver.B_t.graph.row_map);
-        Kokkos::parallel_for(
-            "CopySparseValuesToTranspose", constraint_transpose_policy,
-            CopySparseValuesToTranspose<Solver::CrsMatrixType>{solver.B, tmp_row_map, solver.B_t}
-        );
-        KokkosSparse::sort_crs_matrix(solver.B_t);
-    }
+struct Solver;
+struct Constraints;
 
-    {
-        auto mult_region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Matrix");
-        KokkosSparse::spgemm_numeric(
-            solver.constraints_spgemm_handle, solver.B, false, solver.T, false,
-            solver.constraints_matrix
-        );
-    }
-}
+void AssembleConstraintsMatrix(Solver& solver, const Constraints& constraints);
 }  // namespace openturbine

src/step/assemble_constraints_residual.cpp

@@ -0,0 +1,49 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_Profiling_ScopedRegion.hpp>
+
+#include "src/constraints/constraints.hpp"
+#include "src/solver/contribute_constraints_system_residual_to_vector.hpp"
+#include "src/solver/copy_constraints_residual_to_vector.hpp"
+#include "src/solver/solver.hpp"
+
+namespace openturbine {
+
+void AssembleConstraintsResidual(const Solver& solver, const Constraints& constraints) {
+    auto resid_region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Residual");
+
+    if (constraints.num == 0) {
+        return;
+    }
+
+    Kokkos::parallel_for(
+        "ContributeConstraintsSystemResidualToVector", constraints.num,
+        ContributeConstraintsSystemResidualToVector{
+            constraints.target_node_index, solver.R, constraints.system_residual_terms}
+    );
+
+    auto R = Solver::ValuesType(
+        Kokkos::view_alloc(Kokkos::WithoutInitializing, "R_local"), solver.num_system_dofs
+    );
+    Kokkos::deep_copy(
+        R, Kokkos::subview(solver.R, Kokkos::make_pair(size_t{0U}, solver.num_system_dofs))
+    );
+    auto lambda = Solver::ValuesType(
+        Kokkos::view_alloc(Kokkos::WithoutInitializing, "lambda"), constraints.lambda.extent(0)
+    );
+    Kokkos::deep_copy(lambda, constraints.lambda);
+    auto spmv_handle = Solver::SpmvHandle();
+    KokkosSparse::spmv(&spmv_handle, "T", 1., solver.B, lambda, 1., R);
+    Kokkos::deep_copy(
+        Kokkos::subview(solver.R, Kokkos::make_pair(size_t{0U}, solver.num_system_dofs)), R
+    );
+
+    Kokkos::parallel_for(
+        "CopyConstraintsResidualToVector", constraints.num,
+        CopyConstraintsResidualToVector{
+            constraints.row_range,
+            Kokkos::subview(solver.R, Kokkos::make_pair(solver.num_system_dofs, solver.num_dofs)),
+            constraints.residual_terms}
+    );
+}
+
+}  // namespace openturbine

src/step/assemble_constraints_residual.hpp

@@ -1,51 +1,9 @@
 #pragma once
 
-#include <Kokkos_Core.hpp>
-#include <Kokkos_Profiling_ScopedRegion.hpp>
-
-#include "src/constraints/constraints.hpp"
-#include "src/solver/contribute_constraints_system_residual_to_vector.hpp"
-#include "src/solver/copy_constraints_residual_to_vector.hpp"
-#include "src/solver/solver.hpp"
-
 namespace openturbine {
+struct Solver;
+struct Constraints;
 
-inline void AssembleConstraintsResidual(Solver& solver, Constraints& constraints) {
-    auto resid_region = Kokkos::Profiling::ScopedRegion("Assemble Constraints Residual");
-
-    if (constraints.num == 0) {
-        return;
-    }
-
-    Kokkos::parallel_for(
-        "ContributeConstraintsSystemResidualToVector", constraints.num,
-        ContributeConstraintsSystemResidualToVector{
-            constraints.target_node_index, solver.R, constraints.system_residual_terms}
-    );
-
-    auto R = Solver::ValuesType(
-        Kokkos::view_alloc(Kokkos::WithoutInitializing, "R_local"), solver.num_system_dofs
-    );
-    Kokkos::deep_copy(
-        R, Kokkos::subview(solver.R, Kokkos::make_pair(size_t{0U}, solver.num_system_dofs))
-    );
-    auto lambda = Solver::ValuesType(
-        Kokkos::view_alloc(Kokkos::WithoutInitializing, "lambda"), constraints.lambda.extent(0)
-    );
-    Kokkos::deep_copy(lambda, constraints.lambda);
-    auto spmv_handle = Solver::SpmvHandle();
-    KokkosSparse::spmv(&spmv_handle, "T", 1., solver.B, lambda, 1., R);
-    Kokkos::deep_copy(
-        Kokkos::subview(solver.R, Kokkos::make_pair(size_t{0U}, solver.num_system_dofs)), R
-    );
-
-    Kokkos::parallel_for(
-        "CopyConstraintsResidualToVector", constraints.num,
-        CopyConstraintsResidualToVector{
-            constraints.row_range,
-            Kokkos::subview(solver.R, Kokkos::make_pair(solver.num_system_dofs, solver.num_dofs)),
-            constraints.residual_terms}
-    );
-}
+void AssembleConstraintsResidual(const Solver& solver, const Constraints& constraints);
 
 }  // namespace openturbine

src/step/assemble_inertia_matrix.cpp

@@ -0,0 +1,25 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_Profiling_ScopedRegion.hpp>
+
+#include "src/beams/beams.hpp"
+#include "src/system/integrate_inertia_matrix.hpp"
+
+namespace openturbine {
+
+void AssembleInertiaMatrix(const Beams& beams, double beta_prime, double gamma_prime) {
+    auto region = Kokkos::Profiling::ScopedRegion("Assemble Inertia Matrix");
+    auto range_policy = Kokkos::TeamPolicy<>(static_cast<int>(beams.num_elems), Kokkos::AUTO());
+    auto smem = 2 * Kokkos::View<double* [6][6]>::shmem_size(beams.max_elem_qps) +
+                2 * Kokkos::View<double*>::shmem_size(beams.max_elem_qps) +
+                Kokkos::View<double**>::shmem_size(beams.max_elem_nodes, beams.max_elem_qps);
+    range_policy.set_scratch_size(1, Kokkos::PerTeam(smem));
+    Kokkos::parallel_for(
+        "IntegrateInertiaMatrix", range_policy,
+        IntegrateInertiaMatrix{
+            beams.num_nodes_per_element, beams.num_qps_per_element, beams.qp_weight,
+            beams.qp_jacobian, beams.shape_interp, beams.qp_Muu, beams.qp_Guu, beta_prime,
+            gamma_prime, beams.inertia_matrix_terms}
+    );
+}
+
+}  // namespace openturbine

src/step/assemble_inertia_matrix.hpp

@@ -1,27 +1,8 @@
 #pragma once
 
-#include <Kokkos_Core.hpp>
-#include <Kokkos_Profiling_ScopedRegion.hpp>
-
-#include "src/beams/beams.hpp"
-#include "src/system/integrate_inertia_matrix.hpp"
-
 namespace openturbine {
+struct Beams;
 
-inline void AssembleInertiaMatrix(const Beams& beams, double beta_prime, double gamma_prime) {
-    auto region = Kokkos::Profiling::ScopedRegion("Assemble Inertia Matrix");
-    auto range_policy = Kokkos::TeamPolicy<>(static_cast<int>(beams.num_elems), Kokkos::AUTO());
-    auto smem = 2 * Kokkos::View<double* [6][6]>::shmem_size(beams.max_elem_qps) +
-                2 * Kokkos::View<double*>::shmem_size(beams.max_elem_qps) +
-                Kokkos::View<double**>::shmem_size(beams.max_elem_nodes, beams.max_elem_qps);
-    range_policy.set_scratch_size(1, Kokkos::PerTeam(smem));
-    Kokkos::parallel_for(
-        "IntegrateInertiaMatrix", range_policy,
-        IntegrateInertiaMatrix{
-            beams.num_nodes_per_element, beams.num_qps_per_element, beams.qp_weight,
-            beams.qp_jacobian, beams.shape_interp, beams.qp_Muu, beams.qp_Guu, beta_prime,
-            gamma_prime, beams.inertia_matrix_terms}
-    );
-}
+void AssembleInertiaMatrix(const Beams& beams, double beta_prime, double gamma_prime);
 
 }  // namespace openturbine

src/step/assemble_residual_vector.cpp

@@ -0,0 +1,30 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_Profiling_ScopedRegion.hpp>
+
+#include "src/beams/beams.hpp"
+#include "src/system/integrate_residual_vector.hpp"
+
+namespace openturbine {
+
+void AssembleResidualVector(const Beams& beams) {
+    auto region = Kokkos::Profiling::ScopedRegion("Assemble Residual");
+    auto range_policy = Kokkos::TeamPolicy<>(static_cast<int>(beams.num_elems), Kokkos::AUTO());
+    const auto shape_size =
+        Kokkos::View<double**>::shmem_size(beams.max_elem_nodes, beams.max_elem_qps);
+    const auto weight_size = Kokkos::View<double*>::shmem_size(beams.max_elem_qps);
+    const auto node_variable_size = Kokkos::View<double* [6]>::shmem_size(beams.max_elem_nodes);
+    const auto qp_variable_size = Kokkos::View<double* [6]>::shmem_size(beams.max_elem_qps);
+    range_policy.set_scratch_size(
+        1,
+        Kokkos::PerTeam(2 * shape_size + 2 * weight_size + node_variable_size + 4 * qp_variable_size)
+    );
+    Kokkos::parallel_for(
+        "IntegrateResidualVector", range_policy,
+        IntegrateResidualVector{
+            beams.num_nodes_per_element, beams.num_qps_per_element, beams.qp_weight,
+            beams.qp_jacobian, beams.shape_interp, beams.shape_deriv, beams.node_FX, beams.qp_Fc,
+            beams.qp_Fd, beams.qp_Fi, beams.qp_Fg, beams.residual_vector_terms}
+    );
+}
+
+}  // namespace openturbine

src/step/assemble_residual_vector.hpp

@@ -1,32 +1,8 @@
 #pragma once
 
-#include <Kokkos_Core.hpp>
-#include <Kokkos_Profiling_ScopedRegion.hpp>
-
-#include "src/beams/beams.hpp"
-#include "src/system/integrate_residual_vector.hpp"
-
 namespace openturbine {
+struct Beams;
 
-inline void AssembleResidualVector(const Beams& beams) {
-    auto region = Kokkos::Profiling::ScopedRegion("Assemble Residual");
-    auto range_policy = Kokkos::TeamPolicy<>(static_cast<int>(beams.num_elems), Kokkos::AUTO());
-    const auto shape_size =
-        Kokkos::View<double**>::shmem_size(beams.max_elem_nodes, beams.max_elem_qps);
-    const auto weight_size = Kokkos::View<double*>::shmem_size(beams.max_elem_qps);
-    const auto node_variable_size = Kokkos::View<double* [6]>::shmem_size(beams.max_elem_nodes);
-    const auto qp_variable_size = Kokkos::View<double* [6]>::shmem_size(beams.max_elem_qps);
-    range_policy.set_scratch_size(
-        1,
-        Kokkos::PerTeam(2 * shape_size + 2 * weight_size + node_variable_size + 4 * qp_variable_size)
-    );
-    Kokkos::parallel_for(
-        "IntegrateResidualVector", range_policy,
-        IntegrateResidualVector{
-            beams.num_nodes_per_element, beams.num_qps_per_element, beams.qp_weight,
-            beams.qp_jacobian, beams.shape_interp, beams.shape_deriv, beams.node_FX, beams.qp_Fc,
-            beams.qp_Fd, beams.qp_Fi, beams.qp_Fg, beams.residual_vector_terms}
-    );
-}
+void AssembleResidualVector(const Beams& beams);
 
 }  // namespace openturbine