feat(deterministic-parallel-assembly): step 1 — tbb-element-evaluation

CMakeLists.txt

@@ -29,6 +29,7 @@ add_library(fesa_core STATIC
     src/fesa/analysis/linear_static_analysis.cpp
     src/fesa/analysis/run_solver.cpp
     src/fesa/assembly/contribution.cpp
+    src/fesa/assembly/parallel_assembler.cpp
     src/fesa/assembly/serial_assembler.cpp
     src/fesa/constraints/essential_bc.cpp
     src/fesa/core/version.cpp
@@ -55,7 +56,7 @@ target_include_directories(fesa_core
 
 target_compile_features(fesa_core PUBLIC cxx_std_20)
 target_compile_options(fesa_core PRIVATE /W4 /permissive- /EHsc)
-target_link_libraries(fesa_core PRIVATE MKL::MKL HDF5::HDF5)
+target_link_libraries(fesa_core PRIVATE MKL::MKL TBB::tbb HDF5::HDF5)
 
 add_executable(fesa
     src/fesa/cli/main.cpp

include/fesa/assembly/assembler.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <cstddef>
+
+#include <fesa/assembly/equation_system.hpp>
+#include <fesa/fem/dof_manager.hpp>
+#include <fesa/model/domain.hpp>
+
+namespace fesa {
+
+struct AssemblyOptions final {
+    std::size_t max_threads;
+    std::size_t grain_size;
+};
+
+[[nodiscard]] EquationSystem assemble_parallel(
+    const Domain& domain,
+    const DofManager& dofs,
+    AssemblyOptions options);
+
+}  // namespace fesa

src/fesa/assembly/parallel_assembler.cpp

@@ -0,0 +1,208 @@
+#include <fesa/assembly/assembler.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <iterator>
+#include <limits>
+#include <numeric>
+#include <optional>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <fesa/assembly/contribution.hpp>
+#include <fesa/elements/beam/beam3d2.hpp>
+
+#include <oneapi/tbb/blocked_range.h>
+#include <oneapi/tbb/parallel_for.h>
+#include <oneapi/tbb/task_arena.h>
+
+namespace fesa {
+namespace {
+
+struct ElementEvaluation final {
+    std::vector<MatrixContribution> contributions;
+    std::optional<std::string> error;
+};
+
+auto origin_key(const EntityOrigin& origin) {
+    return std::tie(
+        origin.instance_name,
+        origin.local_label,
+        origin.part_name);
+}
+
+std::string kernel_error_message(
+    const BeamElement& element,
+    const BeamKernelResult& result) {
+    std::string message =
+        "Beam element " + std::to_string(element.origin.local_label) +
+        " kernel failed";
+    for (const Diagnostic& diagnostic : result.diagnostics) {
+        message += ": " + diagnostic.code + " - " + diagnostic.message;
+    }
+    return message;
+}
+
+std::vector<std::size_t> canonical_element_order(const Domain& domain) {
+    std::vector<std::size_t> order(domain.beam_elements().size());
+    std::iota(order.begin(), order.end(), std::size_t{0});
+    std::ranges::sort(
+        order,
+        [&domain](const std::size_t left, const std::size_t right) {
+            return origin_key(domain.beam_elements()[left].origin) <
+                   origin_key(domain.beam_elements()[right].origin);
+        });
+    return order;
+}
+
+std::vector<ElementId> canonical_contribution_element_ids(
+    const std::vector<std::size_t>& order) {
+    // Domain ElementIds are not ordered by input identity. These tie-break
+    // IDs encode the existing serial assembler's element-origin order.
+    std::vector<ElementId> ids(order.size(), ElementId{0});
+    for (std::size_t rank = 0; rank < order.size(); ++rank) {
+        ids[order[rank]] = ElementId{static_cast<std::int64_t>(rank)};
+    }
+    return ids;
+}
+
+ElementEvaluation evaluate_element(
+    const Domain& domain,
+    const DofManager& dofs,
+    const BeamElement& element,
+    const ElementId canonical_id) {
+    const BeamKernelResult result = compute_beam3d2({
+        {
+            domain.node(element.nodes[0]).position,
+            domain.node(element.nodes[1]).position,
+        },
+        domain.material(element.material),
+        domain.section(element.section),
+    });
+    if (!result.contribution.has_value()) {
+        return {{}, kernel_error_message(element, result)};
+    }
+
+    ElementEvaluation evaluation;
+    evaluation.contributions.reserve(78);
+    const std::array<std::size_t, 12> full_dofs =
+        dofs.element_full_dofs(element);
+    for (std::size_t local_row = 0; local_row < full_dofs.size();
+         ++local_row) {
+        for (std::size_t local_column = local_row;
+             local_column < full_dofs.size();
+             ++local_column) {
+            evaluation.contributions.push_back({
+                std::min(
+                    full_dofs[local_row],
+                    full_dofs[local_column]),
+                std::max(
+                    full_dofs[local_row],
+                    full_dofs[local_column]),
+                canonical_id,
+                static_cast<std::uint16_t>(
+                    local_row * full_dofs.size() + local_column),
+                result.contribution
+                    ->global_stiffness[local_row][local_column],
+            });
+        }
+    }
+    return evaluation;
+}
+
+std::vector<double> assemble_force(
+    const Domain& domain,
+    const DofManager& dofs) {
+    std::vector<double> force(dofs.full_dof_count(), 0.0);
+    for (const NodalLoad& load : domain.step().nodal_loads) {
+        for (std::size_t component = 0; component < load.values.size();
+             ++component) {
+            const auto dof = static_cast<NodeDof>(component);
+            force[dofs.full_dof({load.node, dof})] +=
+                load.values[component];
+        }
+    }
+    return force;
+}
+
+void validate_options(const AssemblyOptions options) {
+    if (options.max_threads == 0) {
+        throw std::invalid_argument{
+            "Assembly max_threads must be greater than zero."};
+    }
+    if (options.max_threads >
+        static_cast<std::size_t>(std::numeric_limits<int>::max())) {
+        throw std::invalid_argument{
+            "Assembly max_threads exceeds the TBB task arena range."};
+    }
+    if (options.grain_size == 0) {
+        throw std::invalid_argument{
+            "Assembly grain_size must be greater than zero."};
+    }
+}
+
+}  // namespace
+
+EquationSystem assemble_parallel(
+    const Domain& domain,
+    const DofManager& dofs,
+    const AssemblyOptions options) {
+    validate_options(options);
+
+    const std::vector<std::size_t> element_order =
+        canonical_element_order(domain);
+    const std::vector<ElementId> element_ids =
+        canonical_contribution_element_ids(element_order);
+    std::vector<ElementEvaluation> evaluations(
+        domain.beam_elements().size());
+
+    oneapi::tbb::task_arena arena{
+        static_cast<int>(options.max_threads)};
+    arena.execute([&] {
+        oneapi::tbb::parallel_for(
+            oneapi::tbb::blocked_range<std::size_t>{
+                0,
+                evaluations.size(),
+                options.grain_size,
+            },
+            [&](const oneapi::tbb::blocked_range<std::size_t>& range) {
+                for (std::size_t index = range.begin();
+                     index != range.end();
+                     ++index) {
+                    ElementEvaluation local = evaluate_element(
+                        domain,
+                        dofs,
+                        domain.beam_elements()[index],
+                        element_ids[index]);
+                    evaluations[index] = std::move(local);
+                }
+            });
+    });
+
+    std::vector<MatrixContribution> contributions;
+    contributions.reserve(domain.beam_elements().size() * 78);
+    for (const std::size_t element_index : element_order) {
+        ElementEvaluation& evaluation = evaluations[element_index];
+        if (evaluation.error.has_value()) {
+            throw std::runtime_error{std::move(*evaluation.error)};
+        }
+        contributions.insert(
+            contributions.end(),
+            std::make_move_iterator(evaluation.contributions.begin()),
+            std::make_move_iterator(evaluation.contributions.end()));
+    }
+
+    const std::vector<MatrixContribution> canonical =
+        canonicalize_contributions(contributions);
+    return {
+        merge_contributions(dofs.full_dof_count(), canonical),
+        assemble_force(domain, dofs),
+    };
+}
+
+}  // namespace fesa

tests/CMakeLists.txt

@@ -409,6 +409,43 @@ add_test(
         --gtest_filter=DeterministicMerge.*
 )
 
+add_executable(fesa_parallel_assembly_tests
+    unit/assembly/parallel_assembler_test.cpp
+)
+
+target_compile_features(fesa_parallel_assembly_tests PRIVATE cxx_std_20)
+target_compile_options(
+    fesa_parallel_assembly_tests
+    PRIVATE
+        /W4
+        /permissive-
+        /EHsc
+)
+
+target_link_libraries(fesa_parallel_assembly_tests
+    PRIVATE
+        fesa_core
+        GTest::gtest_main
+)
+
+add_test(
+    NAME ParallelAssembly
+    COMMAND "$<TARGET_FILE:fesa_parallel_assembly_tests>"
+        --gtest_filter=ParallelAssembly.*
+)
+
+add_test(
+    NAME TbbElementEvaluation
+    COMMAND "$<TARGET_FILE:fesa_parallel_assembly_tests>"
+        --gtest_filter=TbbElementEvaluation.*
+)
+
+set_property(
+    TEST ParallelAssembly TbbElementEvaluation
+    PROPERTY ENVIRONMENT_MODIFICATION
+        ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS}
+)
+
 add_executable(fesa_constraint_tests
     unit/constraints/essential_bc_test.cpp
 )

tests/unit/assembly/parallel_assembler_test.cpp

@@ -0,0 +1,189 @@
+#include <fesa/assembly/assembler.hpp>
+#include <fesa/assembly/serial_assembler.hpp>
+
+#include <array>
+#include <bit>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <fesa/model/domain_builder.hpp>
+
+#include <gtest/gtest.h>
+
+namespace {
+
+fesa::Domain finish_domain(
+    fesa::DomainBuilder builder,
+    fesa::StepDefinition step) {
+    builder.set_step(std::move(step));
+    auto result = std::move(builder).build();
+    if (!result.domain.has_value()) {
+        throw std::runtime_error{"Test Domain failed validation."};
+    }
+    return std::move(*result.domain);
+}
+
+fesa::Domain build_beam_domain(const bool branched) {
+    fesa::DomainBuilder builder;
+    const std::array<fesa::Vec3, 6> positions{{
+        {0.0, 0.0, 0.0},
+        {1.0, 0.0, 0.0},
+        {2.0, 0.0, 0.0},
+        {3.0, 0.0, 0.0},
+        {1.0, 0.0, 1.0},
+        {1.0, 1.0, 1.0},
+    }};
+    const std::size_t node_count = branched ? positions.size() : 4;
+    for (std::size_t index = 0; index < node_count; ++index) {
+        builder.add_node({
+            fesa::NodeId{static_cast<std::int64_t>(index)},
+            fesa::EntityOrigin{
+                "BeamPart",
+                "Beam-1",
+                static_cast<std::int64_t>(index + 1),
+            },
+            positions[index],
+        });
+    }
+    builder.add_material({
+        fesa::MaterialId{0},
+        "Steel",
+        210.0e9,
+        0.3,
+    });
+    builder.add_section({
+        fesa::SectionId{0},
+        "General",
+        0.02,
+        3.0e-5,
+        4.0e-5,
+        2.0e-5,
+        0.015,
+        0.016,
+        fesa::ShearPropertySource::input,
+        fesa::Vec3{0.0, 1.0, 0.0},
+        {},
+    });
+
+    const std::array<fesa::BeamElement, 5> elements{{
+        {
+            fesa::ElementId{40},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 300},
+            {fesa::NodeId{2}, fesa::NodeId{3}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+        {
+            fesa::ElementId{10},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 100},
+            {fesa::NodeId{0}, fesa::NodeId{1}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+        {
+            fesa::ElementId{30},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 200},
+            {fesa::NodeId{1}, fesa::NodeId{2}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+        {
+            fesa::ElementId{0},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 500},
+            {fesa::NodeId{1}, fesa::NodeId{5}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+        {
+            fesa::ElementId{20},
+            fesa::EntityOrigin{"BeamPart", "Beam-1", 400},
+            {fesa::NodeId{1}, fesa::NodeId{4}},
+            fesa::MaterialId{0},
+            fesa::SectionId{0},
+        },
+    }};
+    const std::size_t element_count = branched ? elements.size() : 3;
+    for (std::size_t index = 0; index < element_count; ++index) {
+        builder.add_beam_element(elements[index]);
+    }
+
+    return finish_domain(
+        std::move(builder),
+        {
+            "Load",
+            {},
+            {
+                {fesa::NodeId{3}, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}},
+                {fesa::NodeId{1}, {1.0e16, 1.0, 0.0, 0.0, 0.0, 0.0}},
+                {fesa::NodeId{1}, {-1.0e16, 2.0, 0.0, 0.0, 0.0, 0.0}},
+            },
+        });
+}
+
+std::vector<std::uint64_t> bits(const std::vector<double>& values) {
+    std::vector<std::uint64_t> result;
+    result.reserve(values.size());
+    for (const double value : values) {
+        result.push_back(std::bit_cast<std::uint64_t>(value));
+    }
+    return result;
+}
+
+void expect_bitwise_equal(
+    const fesa::EquationSystem& expected,
+    const fesa::EquationSystem& actual) {
+    EXPECT_EQ(actual.stiffness.order, expected.stiffness.order);
+    EXPECT_EQ(
+        actual.stiffness.row_offsets,
+        expected.stiffness.row_offsets);
+    EXPECT_EQ(
+        actual.stiffness.column_indices,
+        expected.stiffness.column_indices);
+    EXPECT_EQ(bits(actual.stiffness.values), bits(expected.stiffness.values));
+    EXPECT_EQ(bits(actual.force), bits(expected.force));
+}
+
+TEST(ParallelAssembly, MatchesSerialForFixedBeamChain) {
+    const fesa::Domain domain = build_beam_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+    const fesa::EquationSystem serial = fesa::assemble_serial(domain, dofs);
+
+    expect_bitwise_equal(
+        serial,
+        fesa::assemble_parallel(domain, dofs, {1, 1}));
+    expect_bitwise_equal(
+        serial,
+        fesa::assemble_parallel(domain, dofs, {4, 2}));
+}
+
+TEST(ParallelAssembly, MatchesSerialForFixedBranchedDomain) {
+    const fesa::Domain domain = build_beam_domain(true);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+    const fesa::EquationSystem serial = fesa::assemble_serial(domain, dofs);
+
+    expect_bitwise_equal(
+        serial,
+        fesa::assemble_parallel(domain, dofs, {1, 2}));
+    expect_bitwise_equal(
+        serial,
+        fesa::assemble_parallel(domain, dofs, {3, 1}));
+}
+
+TEST(TbbElementEvaluation, RejectsZeroExecutionLimits) {
+    const fesa::Domain domain = build_beam_domain(false);
+    const fesa::DofManager dofs = fesa::DofManager::build(domain);
+
+    EXPECT_THROW(
+        static_cast<void>(
+            fesa::assemble_parallel(domain, dofs, {0, 1})),
+        std::invalid_argument);
+    EXPECT_THROW(
+        static_cast<void>(
+            fesa::assemble_parallel(domain, dofs, {1, 0})),
+        std::invalid_argument);
+}
+
+}  // namespace