diff --git a/include/fesa/analysis/thermal_analysis_state.h b/include/fesa/analysis/thermal_analysis_state.h new file mode 100644 index 0000000..1114de7 --- /dev/null +++ b/include/fesa/analysis/thermal_analysis_state.h @@ -0,0 +1,95 @@ +#ifndef FESA_ANALYSIS_THERMAL_ANALYSIS_STATE_H_ +#define FESA_ANALYSIS_THERMAL_ANALYSIS_STATE_H_ + +#include +#include +#include + +#include "fesa/analysis/thermal_time_grid.h" +#include "fesa/assembly/thermal_assembler.h" +#include "fesa/results/thermal_result_records.h" + +namespace fesa { + +/// @brief Owns the immutable thermal result inventory and accepted history. +/// @note No Domain/model/assembly references survive Create; state may outlive +/// them. Numerical equilibrium and energy checks belong to the procedure before +/// Commit. +class ThermalAnalysisState final { + public: + /// @brief Validates the inventory/schedule and commits raw initial + /// temperature. + /// @note Initial frame has no recovery payload and is never projected onto + /// BCs. + static Result Create(const AnalysisModel& model, + const DofManager& dofs, + const ThermalAssembly& assembly); + /// @brief Builds a complete typed candidate without changing accepted + /// history. + /// @param temperature Full node-major temperature, including prescribed + /// values. + /// @param reaction Constrained full residual in ReactionDofs() order, even if + /// empty. + /// @param flux Global element fluxes in ElementOrder(), point/trace order. + Result MakeCompletedFrame( + std::uint32_t index, Vector temperature, Vector reaction, + std::vector flux) const; + /// @brief Atomically appends the next finite frame after exact + /// mapping/schedule checks. + /// @note Rejection preserves every existing frame; caller-owned candidates + /// stay independent. + Status CommitAcceptedFrame(const ThermalAcceptedFrame& frame); + /// @brief Returns accepted frames, beginning with raw frame zero. + const std::vector& AcceptedFrames() const noexcept; + /// @brief Returns the most recently committed state; Create always supplies + /// frame zero. + const ThermalAcceptedFrame& CommittedState() const noexcept; + /// @brief Reports whether every scheduled increment has been committed. + bool IsComplete() const noexcept; + /// @brief Returns the exact time-grid schedule used to validate accepted + /// frames. + const ThermalTimeGrid& TimeGrid() const noexcept; + /// @brief Returns original nominal dt for metadata, independent of actual + /// last delta. + double NominalTimeIncrement() const noexcept; + /// @brief Returns the exact input end time for metadata. + double TimePeriod() const noexcept; + /// @brief Returns active Domain node IDs in temperature row order. + const std::vector& NodeOrder() const noexcept; + /// @brief Returns immutable source identities corresponding to NodeOrder(). + const std::vector& NodeSourceIds() const noexcept; + /// @brief Returns Domain element IDs in flux/frame row order. + const std::vector& ElementOrder() const noexcept; + /// @brief Returns immutable source identities corresponding to + /// ElementOrder(). + const std::vector& ElementSourceIds() const noexcept; + /// @brief Returns the N by 5 free=0/prescribed=1 mask. + const std::vector>& ConstraintMask() + const noexcept; + /// @brief Returns prescribed node/source-DOF identities in node/station + /// order. + const std::vector& ReactionDofs() const noexcept; + /// @brief Returns E by 4 fixed local frames, stored once for the entire + /// history. + const std::vector& FluxLocalFrames() const noexcept; + + private: + /// @brief Takes an exact schedule; Create then fills the validated inventory. + ThermalAnalysisState(ThermalTimeGrid grid, double nominal, double period); + /// @brief Checks frame metadata, mapping, payload shape and finite values. + Status ValidateFrame(const ThermalAcceptedFrame& frame) const; + ThermalTimeGrid grid_; + double nominal_; + double period_; + std::vector node_order_; + std::vector node_source_ids_; + std::vector element_order_; + std::vector element_source_ids_; + std::vector> constraint_mask_; + std::vector reaction_dofs_; + std::vector local_frames_; + std::vector frames_; +}; + +} // namespace fesa +#endif // FESA_ANALYSIS_THERMAL_ANALYSIS_STATE_H_ diff --git a/include/fesa/assembly/thermal_assembler.h b/include/fesa/assembly/thermal_assembler.h new file mode 100644 index 0000000..2075987 --- /dev/null +++ b/include/fesa/assembly/thermal_assembler.h @@ -0,0 +1,87 @@ +#ifndef FESA_ASSEMBLY_THERMAL_ASSEMBLER_H_ +#define FESA_ASSEMBLY_THERMAL_ASSEMBLER_H_ + +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/elements/thermal_shell.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/results/thermal_result_records.h" + +namespace fesa { +class ThermalAssembly; + +/// @brief Builds full K,C,Q from checked DS4 kernels in deterministic source +/// order. +/// @note Domain must outlive the returned assembly; neither model nor dofs is +/// retained. Workers write index-owned buffers only. Failures are selected in +/// element order. +Result AssembleThermal(const AnalysisModel& model, + const DofManager& dofs, + const ParallelFor& parallel_for); + +/// @brief Owns full-space operators and local kernels for a fixed thermal +/// Domain. +/// @note The original Domain must remain alive and unmoved for this object's +/// lifetime. +class ThermalAssembly final { + public: + /// @brief Returns full conductivity [power/temperature] in canonical CSR + /// order. + const SparseMatrix& Conductivity() const noexcept; + /// @brief Returns full consistent capacity [energy/temperature]. + const SparseMatrix& Capacity() const noexcept; + /// @brief Returns signed full nodal input [power], including prescribed + /// targets. + const Vector& HeatInput() const noexcept; + /// @brief Bounds local interval error plus ordered summation roundoff + /// entrywise. + /// @note Excludes quadrature truncation and physical model error. + const SparseMatrix& ConductivityError() const noexcept; + /// @copydoc ConductivityError + const SparseMatrix& CapacityError() const noexcept; + /// @brief Bounds source-row input summation roundoff entrywise, in power + /// units. + const Vector& HeatInputError() const noexcept; + /// @brief Returns original Domain element IDs in local-kernel/flux row order. + const std::vector& ElementOrder() const noexcept; + /// @brief Returns checked kernels in stable active element order. + const std::vector& Kernels() const noexcept; + /// @brief Returns node-major twenty-entry full DOF scatters in kernel order. + const std::vector>& Scatters() const noexcept; + /// @brief Checks source Domain identity and the complete supplied equation + /// map. + Status ValidateFor(const AnalysisModel& model, const DofManager& dofs) const; + /// @brief Recovers ordered element-global fluxes from a finite full + /// temperature. + Result> RecoverFlux( + const Vector& temperature, const ParallelFor& parallel_for) const; + + private: + friend Result AssembleThermal(const AnalysisModel&, + const DofManager&, + const ParallelFor&); + /// @brief Takes completed immutable operators after validation and reduction. + ThermalAssembly(const Domain& domain, SparseMatrix conductivity, + SparseMatrix capacity, Vector input, + SparseMatrix conductivity_error, SparseMatrix capacity_error, + Vector input_error, std::vector element_order, + std::vector kernels, + std::vector> scatters); + const Domain* domain_; + SparseMatrix conductivity_; + SparseMatrix capacity_; + Vector input_; + SparseMatrix conductivity_error_; + SparseMatrix capacity_error_; + Vector input_error_; + std::vector element_order_; + std::vector kernels_; + std::vector> scatters_; +}; + +} // namespace fesa +#endif // FESA_ASSEMBLY_THERMAL_ASSEMBLER_H_ diff --git a/include/fesa/results/thermal_result_records.h b/include/fesa/results/thermal_result_records.h new file mode 100644 index 0000000..dcc5119 --- /dev/null +++ b/include/fesa/results/thermal_result_records.h @@ -0,0 +1,53 @@ +#ifndef FESA_RESULTS_THERMAL_RESULT_RECORDS_H_ +#define FESA_RESULTS_THERMAL_RESULT_RECORDS_H_ + +#include +#include +#include +#include + +#include "fesa/core/source_identity.h" +#include "fesa/math/vector.h" +#include "fesa/math/vector3.h" + +namespace fesa { + +/// @brief Identifies a prescribed temperature by Domain node and source DOF +/// 11/15. +struct ThermalReactionDof { + EntityIndex node; + int source_dof; +}; + +/// @brief Holds four CCW surface points times six ordered global flux traces. +/// @note Row=6*point+trace; traces are -1+, -.5, 0-, 0+, .5, 1-. +using ThermalElementFlux = std::array; +/// @brief Holds R columns e1,e2,n at each of four CCW surface points. +using ThermalElementFrames = std::array, 4>; + +/// @brief Holds recovery present only at completed increments, including P=0. +struct ThermalCompletedPayload { + std::vector reaction_dofs; + Vector reaction_heat_rate{0U}; // Power, positive into the model. + std::vector element_order; + std::vector global_heat_flux; // Power/area. +}; + +/// @brief Carries an initial or completed thermal frame with explicit row +/// mapping. +/// @note Source identities are resolved through the state's immutable +/// inventory. +struct ThermalAcceptedFrame { + std::uint64_t frame_index = 0U; + std::uint64_t increment_index = 0U; + double time = 0; + double time_increment = 0; + bool is_initial = true; + std::vector node_order; + Vector temperature{ + 0U}; // Node-major, source DOFs 11..15; signed temperature. + std::optional completed; +}; + +} // namespace fesa +#endif // FESA_RESULTS_THERMAL_RESULT_RECORDS_H_ diff --git a/src/fesa/analysis/thermal_analysis_state.cpp b/src/fesa/analysis/thermal_analysis_state.cpp new file mode 100644 index 0000000..1e2e034 --- /dev/null +++ b/src/fesa/analysis/thermal_analysis_state.cpp @@ -0,0 +1,229 @@ +#include "fesa/analysis/thermal_analysis_state.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +/// @brief Rejects a candidate without touching accepted state. +Status Failure(const std::string& message) { + return Status::Failure(FailureCategory::kSolver, {{Severity::kError, + "invalid-thermal-state", + {}, + "HEAT TRANSFER", + "", + message}}); +} + +/// @brief Checks every stored scalar without imposing a temperature sign bound. +bool Finite(const Vector& vector) { + for (std::size_t i = 0; i < vector.Size(); ++i) + if (!std::isfinite(vector[i])) return false; + return true; +} + +/// @brief Compares prescribed row identities without conflating node and row +/// IDs. +bool SameReactionDofs(const std::vector& left, + const std::vector& right) { + if (left.size() != right.size()) return false; + for (std::size_t i = 0; i < left.size(); ++i) + if (left[i].node != right[i].node || + left[i].source_dof != right[i].source_dof) + return false; + return true; +} + +} // namespace + +Result ThermalAnalysisState::Create( + const AnalysisModel& model, const DofManager& dofs, + const ThermalAssembly& assembly) { + const Status validation = assembly.ValidateFor(model, dofs); + if (!validation.IsOk()) + return Result::Failure(validation); + const auto& step = *model.ThermalStep(); + const auto grid = ThermalTimeGrid::Create( + step.time_period, step.nominal_increment, step.maximum_increments); + if (!grid.HasValue()) + return Result::Failure(grid.GetStatus()); + ThermalAnalysisState state{grid.Value(), step.nominal_increment, + step.time_period}; + const auto& domain = model.GetDomain(); + state.node_order_ = model.ActiveThermalNodes(); + state.element_order_ = assembly.ElementOrder(); + state.node_source_ids_.reserve(state.node_order_.size()); + state.element_source_ids_.reserve(state.element_order_.size()); + state.constraint_mask_.resize(state.node_order_.size()); + for (const auto node : state.node_order_) + state.node_source_ids_.push_back(domain.Nodes()[node].source_id); + for (std::size_t index = 0; index < state.element_order_.size(); ++index) { + state.element_source_ids_.push_back(assembly.Kernels()[index].SourceId()); + const auto frame = assembly.Kernels()[index].Recover({}); + if (!frame.HasValue()) + return Result::Failure(frame.GetStatus()); + state.local_frames_.push_back(frame.Value().frames); + } + for (const auto full : dofs.ConstrainedDofs()) { + const auto row = full / 5; + const auto station = full % 5; + if (row >= state.node_order_.size() || (station != 0 && station != 4)) + return Result::Failure(Failure( + "Only active source DOFs 11/15 may carry thermal reactions.")); + state.constraint_mask_[row][station] = 1; + state.reaction_dofs_.push_back( + {state.node_order_[row], static_cast(station) + 11}); + } + ThermalAcceptedFrame initial; + initial.node_order = state.node_order_; + initial.temperature = Vector{dofs.FullDofCount()}; + for (const auto& declaration : domain.ThermalInitialConditions()) { + for (int station = 11; station <= 15; ++station) { + const auto full = dofs.FullThermalDof( + declaration.node_index, static_cast(station)); + initial.temperature[full] = + declaration.station_values[static_cast(station - 11)]; + } + } + const Status valid_initial = state.ValidateFrame(initial); + if (!valid_initial.IsOk()) + return Result::Failure(valid_initial); + // No reserve(count+1): INC may be UINT32_MAX; history grows only upon commit. + state.frames_.push_back(std::move(initial)); + return Result::Success(std::move(state)); +} + +ThermalAnalysisState::ThermalAnalysisState(ThermalTimeGrid grid, double nominal, + double period) + : grid_{grid}, nominal_{nominal}, period_{period} {} + +Status ThermalAnalysisState::ValidateFrame( + const ThermalAcceptedFrame& frame) const { + if (frame.frame_index != frame.increment_index || + frame.frame_index > grid_.Count() || frame.node_order != node_order_ || + frame.temperature.Size() != node_order_.size() * 5 || + !Finite(frame.temperature)) + return Failure( + "Frame indices, temperature rows or finite values do not match the " + "inventory."); + if (frame.frame_index == 0) { + if (!frame.is_initial || frame.time != 0 || frame.time_increment != 0 || + frame.completed) + return Failure( + "Initial frame must have time/delta zero and temperature only."); + return Status::Ok(); + } + const auto scheduled = + grid_.Frame(static_cast(frame.frame_index)); + if (!scheduled.HasValue()) return scheduled.GetStatus(); + if (frame.is_initial || !frame.completed || + frame.time != scheduled.Value().time || + frame.time_increment != scheduled.Value().delta) + return Failure( + "Completed frame must preserve its exact scheduled time, actual delta " + "and recovery payload."); + const auto& payload = *frame.completed; + if (!SameReactionDofs(payload.reaction_dofs, reaction_dofs_) || + payload.reaction_heat_rate.Size() != reaction_dofs_.size() || + !Finite(payload.reaction_heat_rate) || + payload.element_order != element_order_ || + payload.global_heat_flux.size() != element_order_.size()) + return Failure( + "Completed recovery rows do not match the reaction/element inventory."); + for (const auto& element : payload.global_heat_flux) + for (const auto& flux : element) + if (!flux.IsFinite()) + return Failure("Completed global heat flux is not finite."); + return Status::Ok(); +} + +Result ThermalAnalysisState::MakeCompletedFrame( + std::uint32_t index, Vector temperature, Vector reaction, + std::vector flux) const { + const auto time = grid_.Frame(index); + if (!time.HasValue()) + return Result::Failure(time.GetStatus()); + ThermalAcceptedFrame frame; + frame.frame_index = index; + frame.increment_index = index; + frame.time = time.Value().time; + frame.time_increment = time.Value().delta; + frame.is_initial = false; + frame.node_order = node_order_; + frame.temperature = std::move(temperature); + frame.completed = ThermalCompletedPayload{reaction_dofs_, std::move(reaction), + element_order_, std::move(flux)}; + const Status validation = ValidateFrame(frame); + if (!validation.IsOk()) + return Result::Failure(validation); + return Result::Success(std::move(frame)); +} + +Status ThermalAnalysisState::CommitAcceptedFrame( + const ThermalAcceptedFrame& frame) { + if (frame.frame_index != static_cast(frames_.size()) || + IsComplete()) + return Failure("Only the next scheduled completed frame can be committed."); + const Status validation = ValidateFrame(frame); + if (!validation.IsOk()) return validation; + // std::vector copy insertion provides the strong guarantee for this copyable + // value record: neither validation failure nor an allocation exception alters + // history. + frames_.push_back(frame); + return Status::Ok(); +} +const std::vector& ThermalAnalysisState::AcceptedFrames() + const noexcept { + return frames_; +} +const ThermalAcceptedFrame& ThermalAnalysisState::CommittedState() + const noexcept { + return frames_.back(); +} +bool ThermalAnalysisState::IsComplete() const noexcept { + return static_cast(frames_.size()) == + std::uint64_t{grid_.Count()} + 1; +} +const ThermalTimeGrid& ThermalAnalysisState::TimeGrid() const noexcept { + return grid_; +} +double ThermalAnalysisState::NominalTimeIncrement() const noexcept { + return nominal_; +} +double ThermalAnalysisState::TimePeriod() const noexcept { return period_; } +const std::vector& ThermalAnalysisState::NodeOrder() + const noexcept { + return node_order_; +} +const std::vector& ThermalAnalysisState::NodeSourceIds() + const noexcept { + return node_source_ids_; +} +const std::vector& ThermalAnalysisState::ElementOrder() + const noexcept { + return element_order_; +} +const std::vector& ThermalAnalysisState::ElementSourceIds() + const noexcept { + return element_source_ids_; +} +const std::vector>& +ThermalAnalysisState::ConstraintMask() const noexcept { + return constraint_mask_; +} +const std::vector& ThermalAnalysisState::ReactionDofs() + const noexcept { + return reaction_dofs_; +} +const std::vector& ThermalAnalysisState::FluxLocalFrames() + const noexcept { + return local_frames_; +} + +} // namespace fesa diff --git a/src/fesa/assembly/thermal_assembler.cpp b/src/fesa/assembly/thermal_assembler.cpp new file mode 100644 index 0000000..f67e755 --- /dev/null +++ b/src/fesa/assembly/thermal_assembler.cpp @@ -0,0 +1,374 @@ +#include "fesa/assembly/thermal_assembler.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/elements/thermal_element_layout.h" + +namespace fesa { +namespace { + +/// @brief Preserves source provenance for a failed assembly candidate. +Status Failure(const std::string& message, const SourceLocation& location = {}, + const std::string& identity = "") { + return Status::Failure(FailureCategory::kModel, + {{Severity::kError, "invalid-thermal-assembly", + location, "DS4", identity, message}}); +} + +/// @brief Compares the supplied map with the owner's mapping of this exact +/// model. +Status ValidateDofs(const AnalysisModel& model, const DofManager& dofs) { + if (!model.IsThermal() || model.ThermalStep() == nullptr || + dofs.Profile() != DofProfile::kThermal) { + return Failure("A thermal model and thermal DOF profile are required."); + } + const Status invariant = dofs.ValidateInvariants(); + if (!invariant.IsOk()) return invariant; + const auto expected = DofManager::CreateThermal(model); + if (!expected.HasValue()) return expected.GetStatus(); + const auto& reference = expected.Value(); + if (dofs.FullDofCount() != reference.FullDofCount() || + dofs.FreeDofs() != reference.FreeDofs() || + dofs.ConstrainedDofs() != reference.ConstrainedDofs() || + dofs.GetSparsePattern().row_offsets != + reference.GetSparsePattern().row_offsets || + dofs.GetSparsePattern().column_indices != + reference.GetSparsePattern().column_indices || + dofs.PrescribedValues().Size() != reference.PrescribedValues().Size()) { + return Failure( + "Supplied DOF partition or CSR pattern does not match the model."); + } + for (std::size_t i = 0; i < reference.PrescribedValues().Size(); ++i) { + if (dofs.PrescribedValues()[i] != reference.PrescribedValues()[i]) + return Failure("Supplied prescribed values do not match the model."); + } + for (const EntityIndex node : model.ActiveThermalNodes()) { + for (int station = 11; station <= 15; ++station) { + try { + if (dofs.FullThermalDof(node, + static_cast(station)) != + reference.FullThermalDof(node, + static_cast(station))) + return Failure("Supplied active-node rows do not match the model."); + } catch (const std::out_of_range&) { + return Failure("Supplied active-node rows omit a thermal node."); + } + } + } + return Status::Ok(); +} + +/// @brief Maintains an outward enclosure of a fixed-order scalar reduction. +struct SumInterval { + double lower = 0; + double upper = 0; + /// @brief Adds one value with its pre-existing absolute error enclosure. + void Add(double value, double error) { + const double infinity = std::numeric_limits::infinity(); + const double low = + error == 0 ? value : std::nextafter(value - error, -infinity); + const double high = + error == 0 ? value : std::nextafter(value + error, infinity); + lower = std::nextafter(lower + low, -infinity); + upper = std::nextafter(upper + high, infinity); + } + /// @brief Bounds deviation from the separately reduced binary64 value. + double Error(double value) const { + if (lower == value && upper == value) return 0; + return std::nextafter( + std::max(std::abs(value - lower), std::abs(upper - value)), + std::numeric_limits::infinity()); + } +}; + +struct LocalContribution { + std::optional kernel; + std::array scatter{}; + std::vector conductivity; + std::vector capacity; + Status status = Status::Ok(); +}; + +/// @brief Encloses element errors and the same scalar reduction as FromCoo. +Result AssembleError(const SparseMatrix& matrix, + const std::vector& local, + const SparsePattern& pattern, + bool capacity) { + std::vector sums(matrix.Values().size()); + for (const auto& contribution : local) { + const Matrix& error = capacity ? contribution.kernel->CapacityError() + : contribution.kernel->ConductivityError(); + const auto& values = + capacity ? contribution.capacity : contribution.conductivity; + for (const auto& value : values) { + const auto first = + pattern.column_indices.begin() + pattern.row_offsets[value.row]; + const auto last = + pattern.column_indices.begin() + pattern.row_offsets[value.row + 1]; + const auto found = std::lower_bound(first, last, value.column); + if (found == last || *found != value.column) + return Result::Failure(Failure( + "Error contribution does not match the owner CSR pattern.")); + const auto position = + static_cast(found - pattern.column_indices.begin()); + sums[position].Add(value.value, + error(value.local_order / 20, value.local_order % 20)); + } + } + std::vector errors; + errors.reserve(sums.size()); + for (std::size_t row = 0; row < matrix.Rows(); ++row) { + for (std::size_t p = pattern.row_offsets[row]; + p < pattern.row_offsets[row + 1]; ++p) { + const double error = sums[p].Error(matrix.Values()[p]); + if (!std::isfinite(error)) + return Result::Failure( + Failure("Assembly roundoff enclosure is not finite.")); + errors.push_back({row, pattern.column_indices[p], error, 0, p}); + } + } + return SparseMatrix::FromCoo(matrix.Rows(), matrix.Columns(), + std::move(errors), pattern); +} + +} // namespace + +Result AssembleThermal(const AnalysisModel& model, + const DofManager& dofs, + const ParallelFor& parallel_for) { + const Status validation = ValidateDofs(model, dofs); + if (!validation.IsOk()) return Result::Failure(validation); + const auto& domain = model.GetDomain(); + const auto& order = model.ActiveThermalElements(); + if (order.size() > (std::numeric_limits::max)() / 400) + return Result::Failure( + Failure("Thermal COO storage is not addressable.")); + std::vector local(order.size()); + parallel_for.Execute(order.size(), [&](std::size_t index) { + auto& output = local[index]; + if (order[index] >= domain.Elements().Size() || + index >= domain.ThermalElements().Size()) { + output.status = Failure("Active element index is outside the Domain."); + return; + } + const auto& element = domain.ThermalElements()[index]; + const std::string identity = element.SourceId().instance_name + ":" + + element.SourceId().source_label_text; + if (&domain.Elements()[order[index]] != &element || + element.ThermalSectionIndex() >= domain.ThermalSections().size()) { + output.status = + Failure("Active DS4 definition order/section is inconsistent.", + element.Location(), identity); + return; + } + const auto& section = + domain.ThermalSections()[element.ThermalSectionIndex()]; + if (section.material_index >= domain.ThermalMaterials().Size() || + element.MaterialIndex() >= domain.Materials().Size()) { + output.status = Failure("DS4 thermal material assignment is invalid.", + element.Location(), identity); + return; + } + const auto& material = domain.ThermalMaterials()[section.material_index]; + if (&domain.Materials()[element.MaterialIndex()] != &material) { + output.status = + Failure("Thermal-local and global material mappings disagree.", + element.Location(), identity); + return; + } + ThermalElementDofLayout layout{element.SourceId(), {}}; + std::array coordinates{}; + for (std::size_t i = 0; i < 4; ++i) { + layout.node_indices[i] = element.NodeIndices()[i]; + coordinates[i] = + Vector3{domain.Nodes()[layout.node_indices[i]].coordinates}; + } + const auto scatter = dofs.ThermalElementScatter(layout); + if (!scatter.HasValue()) { + output.status = scatter.GetStatus(); + return; + } + auto kernel = ThermalShell::Create( + element.SourceId(), layout.node_indices, coordinates, + material.Conductivity(), material.Density(), material.SpecificHeat(), + section.thickness, element.Location()); + if (!kernel.HasValue()) { + output.status = kernel.GetStatus(); + return; + } + if (kernel.Value().NodeIndices() != layout.node_indices || + scatter.Value().size() != 20) { + output.status = + Failure("Kernel and thermal scatter inventories disagree.", + element.Location(), identity); + return; + } + std::copy(scatter.Value().begin(), scatter.Value().end(), + output.scatter.begin()); + output.kernel = std::move(kernel.Value()); + output.conductivity.reserve(400); + output.capacity.reserve(400); + for (std::size_t i = 0; i < 20; ++i) + for (std::size_t j = 0; j < 20; ++j) { + output.conductivity.push_back({output.scatter[i], output.scatter[j], + output.kernel->Conductivity()(i, j), + index, 20 * i + j}); + output.capacity.push_back({output.scatter[i], output.scatter[j], + output.kernel->Capacity()(i, j), index, + 20 * i + j}); + } + }); + std::vector conductivity; + std::vector capacity; + conductivity.reserve(400 * local.size()); + capacity.reserve(400 * local.size()); + for (const auto& contribution : local) { + if (!contribution.status.IsOk()) + return Result::Failure(contribution.status); + conductivity.insert(conductivity.end(), contribution.conductivity.begin(), + contribution.conductivity.end()); + capacity.insert(capacity.end(), contribution.capacity.begin(), + contribution.capacity.end()); + } + auto k = + SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(conductivity), dofs.GetSparsePattern()); + if (!k.HasValue()) return Result::Failure(k.GetStatus()); + auto c = SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(capacity), dofs.GetSparsePattern()); + if (!c.HasValue()) return Result::Failure(c.GetStatus()); + auto ek = AssembleError(k.Value(), local, dofs.GetSparsePattern(), false); + if (!ek.HasValue()) return Result::Failure(ek.GetStatus()); + auto ec = AssembleError(c.Value(), local, dofs.GetSparsePattern(), true); + if (!ec.HasValue()) return Result::Failure(ec.GetStatus()); + Vector input{dofs.FullDofCount()}; + Vector input_error{dofs.FullDofCount()}; + std::vector input_intervals(dofs.FullDofCount()); + const auto& loads = model.ThermalStep()->loads; + for (std::size_t index = 0; index < loads.size(); ++index) { + const auto& load = loads[index]; + if (load.source_order != index || !std::isfinite(load.power) || + (load.dof != 11 && load.dof != 15)) + return Result::Failure(Failure( + "Thermal loads must retain finite source-row order and DOF 11/15.", + load.location)); + std::size_t full = 0; + try { + full = dofs.FullThermalDof(load.node_index, + static_cast(load.dof)); + } catch (const std::out_of_range&) { + return Result::Failure( + Failure("Thermal load target is inactive.", load.location)); + } + input[full] += load.power; + input_intervals[full].Add(load.power, 0); + if (!std::isfinite(input[full])) + return Result::Failure(Failure( + "Source-row heat-input accumulation is not finite.", load.location)); + } + for (std::size_t i = 0; i < input.Size(); ++i) { + input_error[i] = input_intervals[i].Error(input[i]); + if (!std::isfinite(input_error[i])) + return Result::Failure( + Failure("Heat-input summation enclosure is not finite.")); + } + std::vector kernels; + std::vector> scatters; + kernels.reserve(local.size()); + scatters.reserve(local.size()); + for (auto& contribution : local) { + kernels.push_back(std::move(*contribution.kernel)); + scatters.push_back(contribution.scatter); + } + return Result::Success(ThermalAssembly{ + domain, std::move(k.Value()), std::move(c.Value()), std::move(input), + std::move(ek.Value()), std::move(ec.Value()), std::move(input_error), + order, std::move(kernels), std::move(scatters)}); +} + +ThermalAssembly::ThermalAssembly( + const Domain& domain, SparseMatrix conductivity, SparseMatrix capacity, + Vector input, SparseMatrix conductivity_error, SparseMatrix capacity_error, + Vector input_error, std::vector element_order, + std::vector kernels, + std::vector> scatters) + : domain_{&domain}, + conductivity_{std::move(conductivity)}, + capacity_{std::move(capacity)}, + input_{std::move(input)}, + conductivity_error_{std::move(conductivity_error)}, + capacity_error_{std::move(capacity_error)}, + input_error_{std::move(input_error)}, + element_order_{std::move(element_order)}, + kernels_{std::move(kernels)}, + scatters_{std::move(scatters)} {} +const SparseMatrix& ThermalAssembly::Conductivity() const noexcept { + return conductivity_; +} +const SparseMatrix& ThermalAssembly::Capacity() const noexcept { + return capacity_; +} +const Vector& ThermalAssembly::HeatInput() const noexcept { return input_; } +const SparseMatrix& ThermalAssembly::ConductivityError() const noexcept { + return conductivity_error_; +} +const SparseMatrix& ThermalAssembly::CapacityError() const noexcept { + return capacity_error_; +} +const Vector& ThermalAssembly::HeatInputError() const noexcept { + return input_error_; +} +const std::vector& ThermalAssembly::ElementOrder() const noexcept { + return element_order_; +} +const std::vector& ThermalAssembly::Kernels() const noexcept { + return kernels_; +} +const std::vector>& ThermalAssembly::Scatters() + const noexcept { + return scatters_; +} +Status ThermalAssembly::ValidateFor(const AnalysisModel& model, + const DofManager& dofs) const { + if (&model.GetDomain() != domain_) + return Failure("Assembly belongs to a different Domain."); + return ValidateDofs(model, dofs); +} +Result> ThermalAssembly::RecoverFlux( + const Vector& temperature, const ParallelFor& parallel_for) const { + if (temperature.Size() != conductivity_.Rows()) + return Result>::Failure(Failure( + "Recovery temperature does not match the full thermal inventory.")); + for (std::size_t i = 0; i < temperature.Size(); ++i) + if (!std::isfinite(temperature[i])) + return Result>::Failure( + Failure("Recovery temperature is not finite.")); + std::vector flux(kernels_.size()); + std::vector status(kernels_.size(), Status::Ok()); + parallel_for.Execute(kernels_.size(), [&](std::size_t index) { + std::array local{}; + for (std::size_t p = 0; p < 20; ++p) + local[p] = temperature[scatters_[index][p]]; + const auto result = kernels_[index].Recover(local); + if (!result.HasValue()) { + status[index] = result.GetStatus(); + return; + } + flux[index] = result.Value().global_flux; + }); + for (const auto& result : status) + if (!result.IsOk()) + return Result>::Failure(result); + return Result>::Success(std::move(flux)); +} + +} // namespace fesa diff --git a/tests/unit/analysis/thermal_analysis_state_test.cpp b/tests/unit/analysis/thermal_analysis_state_test.cpp new file mode 100644 index 0000000..c1f58b8 --- /dev/null +++ b/tests/unit/analysis/thermal_analysis_state_test.cpp @@ -0,0 +1,301 @@ +#include "fesa/analysis/thermal_analysis_state.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { +fesa::ModelDefinition Definition(bool constrained = true, double period = .9, + double increment = .3, + std::uint32_t maximum = 3) { + fesa::ModelDefinition model{}; + model.nodes = {{{"unused", 99, "0099"}, {8, 8, 8}, {}}, + {{"plate", 7, "0007"}, {0, 0, 0}, {}}, + {{"plate", 2, "0002"}, {1, 0, 0}, {}}, + {{"plate", 9, "0009"}, {1, 1, 0}, {}}, + {{"plate", 1, "0001"}, {0, 1, 0}, {}}}; + model.thermal_materials = {{{"plate", 1, "1"}, "heat", 1, 1, 1, {}}}; + model.thermal_sections = {{{"plate", 1, "1"}, "section", 2, 0, {}}}; + model.thermal_elements = {{{"plate", 13, "0013"}, {{1, 2, 3, 4}}, 0, {}}}; + model.thermal_initial_conditions = {{1, -20, {}}, {4, 7.125, {}}}; + model.thermal_step = + fesa::ThermalStepDefinition{"named-step", increment, period, maximum}; + if (constrained) { + model.thermal_model_boundaries = {{1, 11, 100, {}}, {4, 15, 20, {}}}; + model.thermal_step->boundaries = {{1, 11, 200, {}}}; + } + return model; +} + +// Domain, model, DOFs, and assembly all die before any returned state is +// consumed. +fesa::Result State(bool constrained = true, + double period = .9, + double increment = .3, + std::uint32_t maximum = 3) { + auto domain = + fesa::Domain::Create(Definition(constrained, period, increment, maximum)); + if (!domain.HasValue()) + return fesa::Result::Failure( + domain.GetStatus()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + if (!model.HasValue()) + return fesa::Result::Failure(model.GetStatus()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + if (!dofs.HasValue()) + return fesa::Result::Failure(dofs.GetStatus()); + auto assembly = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::SerialParallelFor{}); + if (!assembly.HasValue()) + return fesa::Result::Failure( + assembly.GetStatus()); + return fesa::ThermalAnalysisState::Create(model.Value(), dofs.Value(), + assembly.Value()); +} + +fesa::Result Candidate( + const fesa::ThermalAnalysisState& state, std::uint32_t index) { + fesa::Vector temperature{state.NodeOrder().size() * 5}; + for (std::size_t i = 0; i < temperature.Size(); ++i) + temperature[i] = static_cast(i) - 4.5; + fesa::Vector reaction{state.ReactionDofs().size(), -3}; + std::vector flux(state.ElementOrder().size()); + for (std::size_t e = 0; e < flux.size(); ++e) + for (std::size_t p = 0; p < 24; ++p) + flux[e][p] = fesa::Vector3{static_cast(p), -2, 3}; + return state.MakeCompletedFrame(index, std::move(temperature), + std::move(reaction), std::move(flux)); +} + +std::vector Snapshot(const fesa::ThermalAnalysisState& state) { + std::vector result; + const auto scalar = [&](double x) { + std::uint64_t bits = 0; + std::memcpy(&bits, &x, sizeof(bits)); + result.push_back(bits); + }; + const auto vector = [&](const fesa::Vector& v) { + result.push_back(v.Size()); + for (std::size_t i = 0; i < v.Size(); ++i) scalar(v[i]); + }; + for (const auto& frame : state.AcceptedFrames()) { + result.insert(result.end(), + {frame.frame_index, frame.increment_index, frame.is_initial, + frame.completed.has_value()}); + scalar(frame.time); + scalar(frame.time_increment); + result.push_back(frame.node_order.size()); + for (const auto node : frame.node_order) result.push_back(node); + vector(frame.temperature); + if (frame.completed) { + result.push_back(frame.completed->reaction_dofs.size()); + for (const auto& dof : frame.completed->reaction_dofs) { + result.push_back(dof.node); + result.push_back(dof.source_dof); + } + vector(frame.completed->reaction_heat_rate); + result.push_back(frame.completed->element_order.size()); + for (const auto e : frame.completed->element_order) result.push_back(e); + result.push_back(frame.completed->global_heat_flux.size()); + for (const auto& element : frame.completed->global_heat_flux) + for (const auto& flux : element) + for (const auto value : flux.Components()) scalar(value); + } + } + return result; +} +} // namespace + +TEST(ThermalAnalysisState, + OwnsRawInitialInventoryAfterDomainAndAssemblyDestruction) { + auto state = State(); + ASSERT_TRUE(state.HasValue()); + EXPECT_EQ(state.Value().NodeOrder(), + (std::vector{1, 2, 3, 4})); + EXPECT_EQ(state.Value().NodeSourceIds()[0].source_label_text, "0007"); + EXPECT_EQ(state.Value().ElementSourceIds()[0].source_label_text, "0013"); + EXPECT_EQ(state.Value().ElementOrder(), (std::vector{0})); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 1U); + const auto& initial = state.Value().CommittedState(); + EXPECT_TRUE(initial.is_initial); + EXPECT_FALSE(initial.completed.has_value()); + EXPECT_EQ(initial.time, 0); + EXPECT_EQ(initial.time_increment, 0); + for (std::size_t a = 0; a < 5; ++a) { + EXPECT_EQ(initial.temperature[a], -20.); + EXPECT_EQ(initial.temperature[5 + a], 0.); + EXPECT_EQ(initial.temperature[10 + a], 0.); + EXPECT_EQ(initial.temperature[15 + a], 7.125); + } + ASSERT_EQ(state.Value().ReactionDofs().size(), 2U); + EXPECT_EQ(state.Value().ReactionDofs()[0].node, 1U); + EXPECT_EQ(state.Value().ReactionDofs()[0].source_dof, 11); + EXPECT_EQ(state.Value().ReactionDofs()[1].node, 4U); + EXPECT_EQ(state.Value().ReactionDofs()[1].source_dof, 15); + for (std::size_t row = 0; row < 4; ++row) + for (std::size_t a = 0; a < 5; ++a) + EXPECT_EQ(state.Value().ConstraintMask()[row][a], + ((row == 0 && a == 0) || (row == 3 && a == 4)) ? 1 : 0); + ASSERT_EQ(state.Value().FluxLocalFrames().size(), 1U); + for (const auto& frame : state.Value().FluxLocalFrames()[0]) { + EXPECT_EQ(frame[0], (fesa::Vector3{1, 0, 0})); + EXPECT_EQ(frame[1], (fesa::Vector3{0, 1, 0})); + EXPECT_EQ(frame[2], (fesa::Vector3{0, 0, 1})); + } +} + +TEST(ThermalAnalysisState, CompletedZeroReactionPayloadAndExactActualDelta) { + auto state = State(false); + ASSERT_TRUE(state.HasValue()); + EXPECT_TRUE(state.Value().ReactionDofs().empty()); + EXPECT_FALSE(state.Value().IsComplete()); + for (std::uint32_t i = 1; i <= 3; ++i) { + auto frame = Candidate(state.Value(), i); + ASSERT_TRUE(frame.HasValue()); + ASSERT_TRUE(frame.Value().completed.has_value()); + EXPECT_EQ(frame.Value().completed->reaction_heat_rate.Size(), 0U); + EXPECT_EQ(frame.Value().time_increment, .3); + ASSERT_TRUE(state.Value().CommitAcceptedFrame(frame.Value()).IsOk()); + } + EXPECT_TRUE(state.Value().IsComplete()); + EXPECT_EQ(state.Value().AcceptedFrames().size(), 4U); + EXPECT_EQ(state.Value().CommittedState().time, .9); + EXPECT_EQ(state.Value().CommittedState().time_increment, .3); + EXPECT_NE(state.Value().CommittedState().time - + state.Value().AcceptedFrames()[2].time, + .3); + EXPECT_EQ(state.Value().NominalTimeIncrement(), .3); + EXPECT_EQ(state.Value().TimePeriod(), .9); + EXPECT_FALSE( + state.Value().CommitAcceptedFrame(state.Value().CommittedState()).IsOk()); + auto short_state = State(false, 10, 3, 4); + ASSERT_TRUE(short_state.HasValue()); + auto last = Candidate(short_state.Value(), 4); + ASSERT_TRUE(last.HasValue()); + EXPECT_EQ(last.Value().time, 10.); + EXPECT_EQ(last.Value().time_increment, 1.); +} + +TEST(ThermalAnalysisState, + EveryRejectedCandidatePreservesEntireAcceptedHistory) { + auto created = State(); + ASSERT_TRUE(created.HasValue()); + auto& state = created.Value(); + auto first = Candidate(state, 1); + ASSERT_TRUE(first.HasValue()); + ASSERT_TRUE(state.CommitAcceptedFrame(first.Value()).IsOk()); + const auto snapshot = Snapshot(state); + const std::vector> + corruptions{ + [](auto& f) { + f.frame_index = 3; + f.increment_index = 3; + }, + [](auto& f) { f.increment_index = 1; }, + [](auto& f) { f.is_initial = true; }, + [](auto& f) { f.time += .01; }, + [](auto& f) { f.time_increment = .2; }, + [](auto& f) { f.time = std::numeric_limits::quiet_NaN(); }, + [](auto& f) { + f.time_increment = std::numeric_limits::infinity(); + }, + [](auto& f) { f.node_order.pop_back(); }, + [](auto& f) { f.node_order.push_back(0); }, + [](auto& f) { std::swap(f.node_order[0], f.node_order[1]); }, + [](auto& f) { f.node_order[1] = f.node_order[0]; }, + [](auto& f) { f.temperature = fesa::Vector{19}; }, + [](auto& f) { + f.temperature[0] = std::numeric_limits::quiet_NaN(); + }, + [](auto& f) { f.completed.reset(); }, + [](auto& f) { f.completed->reaction_dofs.pop_back(); }, + [](auto& f) { f.completed->reaction_dofs[0].source_dof = 12; }, + [](auto& f) { + std::swap(f.completed->reaction_dofs[0], + f.completed->reaction_dofs[1]); + }, + [](auto& f) { + f.completed->reaction_dofs[1] = f.completed->reaction_dofs[0]; + }, + [](auto& f) { f.completed->reaction_heat_rate = fesa::Vector{3}; }, + [](auto& f) { + f.completed->reaction_heat_rate[0] = + std::numeric_limits::infinity(); + }, + [](auto& f) { f.completed->element_order[0] = 1; }, + [](auto& f) { f.completed->element_order.push_back(0); }, + [](auto& f) { f.completed->global_heat_flux.clear(); }, + [](auto& f) { + f.completed->global_heat_flux[0][9] = { + 0, 0, std::numeric_limits::quiet_NaN()}; + }}; + for (const auto& corrupt : corruptions) { + auto candidate = Candidate(state, 2); + ASSERT_TRUE(candidate.HasValue()); + corrupt(candidate.Value()); + EXPECT_FALSE(state.CommitAcceptedFrame(candidate.Value()).IsOk()); + EXPECT_EQ(Snapshot(state), snapshot); + } + auto second = Candidate(state, 2); + ASSERT_TRUE(second.HasValue()); + ASSERT_TRUE(state.CommitAcceptedFrame(second.Value()).IsOk()); + second.Value().temperature[0] = 900; + EXPECT_NE(state.CommittedState().temperature[0], 900); +} + +TEST(ThermalAnalysisState, + CheckedCandidateRejectsMissingDataAndHugeGridDoesNotPreallocateHistory) { + auto state = State(false, 4294967295., 1., 4294967295U); + ASSERT_TRUE(state.HasValue()); + EXPECT_EQ(state.Value().TimeGrid().Count(), 4294967295U); + EXPECT_EQ(state.Value().AcceptedFrames().size(), 1U); + EXPECT_FALSE(state.Value().IsComplete()); + auto last = Candidate(state.Value(), 4294967295U); + ASSERT_TRUE(last.HasValue()); + EXPECT_EQ(last.Value().frame_index, 4294967295ULL); + EXPECT_FALSE(state.Value().CommitAcceptedFrame(last.Value()).IsOk()); + EXPECT_FALSE( + state.Value() + .MakeCompletedFrame(0, fesa::Vector{20}, fesa::Vector{0}, {{}}) + .HasValue()); + EXPECT_FALSE( + state.Value() + .MakeCompletedFrame(1, fesa::Vector{19}, fesa::Vector{0}, {{}}) + .HasValue()); + EXPECT_FALSE( + state.Value() + .MakeCompletedFrame(1, fesa::Vector{20}, fesa::Vector{1}, {{}}) + .HasValue()); + EXPECT_FALSE(state.Value() + .MakeCompletedFrame(1, fesa::Vector{20}, fesa::Vector{0}, {}) + .HasValue()); +} + +TEST(ThermalAnalysisState, RejectsAssemblyFromAnotherDomainAndInvalidSchedule) { + auto first = fesa::Domain::Create(Definition()); + ASSERT_TRUE(first.HasValue()); + auto second = fesa::Domain::Create(Definition()); + ASSERT_TRUE(second.HasValue()); + auto fm = fesa::AnalysisModel::Create(first.Value()); + auto sm = fesa::AnalysisModel::Create(second.Value()); + ASSERT_TRUE(fm.HasValue()); + ASSERT_TRUE(sm.HasValue()); + auto fd = fesa::DofManager::CreateThermal(fm.Value()); + auto sd = fesa::DofManager::CreateThermal(sm.Value()); + ASSERT_TRUE(fd.HasValue()); + ASSERT_TRUE(sd.HasValue()); + auto assembly = + fesa::AssembleThermal(fm.Value(), fd.Value(), fesa::SerialParallelFor{}); + ASSERT_TRUE(assembly.HasValue()); + EXPECT_FALSE(fesa::ThermalAnalysisState::Create(sm.Value(), sd.Value(), + assembly.Value()) + .HasValue()); + EXPECT_FALSE(State(false, 10, 3, 3).HasValue()); +} diff --git a/tests/unit/assembly/thermal_assembler_test.cpp b/tests/unit/assembly/thermal_assembler_test.cpp new file mode 100644 index 0000000..9445c77 --- /dev/null +++ b/tests/unit/assembly/thermal_assembler_test.cpp @@ -0,0 +1,323 @@ +#include "fesa/assembly/thermal_assembler.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { +fesa::ModelDefinition Model() { + fesa::ModelDefinition model{}; + const std::array, 7> points{{{90, 90, 90}, + {0, 0, 0}, + {1, 0, 0}, + {2, 0, 0}, + {0, 1, 0}, + {1, 1, 0}, + {2, 1, 0}}}; + const std::array labels{{99, 70, 3, 80, 1, 40, 7}}; + for (std::size_t i = 0; i < points.size(); ++i) + model.nodes.push_back({{"plate", labels[i], std::to_string(labels[i])}, + points[i], + {"mesh.inp", i + 1}}); + model.materials = {{"unused-elastic", 210e9, .3, {}}}; + model.thermal_materials = {{{"plate", 1, "1"}, "first", 2, 3, 4, {}}, + {{"plate", 2, "2"}, "second", 10, 1, 1, {}}}; + model.thermal_sections = {{{"plate", 1, "1"}, "first", 1, 0, {}}, + {{"plate", 2, "2"}, "second", 2, 1, {}}}; + model.thermal_elements = { + {{"plate", 9, "0009"}, {{1, 2, 5, 4}}, 0, {"mesh.inp", 20}}, + {{"plate", 2, "0002"}, {{2, 3, 6, 5}}, 1, {"mesh.inp", 21}}}; + model.thermal_step = fesa::ThermalStepDefinition{"heat", .3, .9, 3}; + model.thermal_step->boundaries = {{2, 11, 100, {}}}; + const double large = std::ldexp(1., 54); + model.thermal_step->loads = {{2, 11, large, {}}, + {2, 11, 1, {}}, + {2, 11, -large, {}}, + {2, 11, 5, {}}, + {1, 15, -7, {}}}; + return model; +} + +fesa::Matrix Dense(const fesa::SparseMatrix& sparse) { + fesa::Matrix dense{sparse.Rows(), sparse.Columns()}; + for (std::size_t i = 0; i < sparse.Rows(); ++i) + for (std::size_t j = sparse.RowOffsets()[i]; j < sparse.RowOffsets()[i + 1]; + ++j) + dense(i, sparse.ColumnIndices()[j]) = sparse.Values()[j]; + return dense; +} + +// Closed-form unit-square polynomial integrals, assembled with hand-fixed +// six-node scatter; no production kernel/layout is used by the expectation. +struct ExpectedMatrices { + fesa::Matrix k{30, 30}; + fesa::Matrix c{30, 30}; + fesa::Matrix pk{30, 30}; + fesa::Matrix pc{30, 30}; +}; +ExpectedMatrices Expected() { + ExpectedMatrices result; + auto& k = result.k; + auto& c = result.c; + constexpr int m[4][4] = { + {4, 2, 1, 2}, {2, 4, 2, 1}, {1, 2, 4, 2}, {2, 1, 2, 4}}; + constexpr int s[4][4] = { + {4, -1, -2, -1}, {-1, 4, -1, -2}, {-2, -1, 4, -1}, {-1, -2, -1, 4}}; + constexpr int z0[5][5] = {{4, 2, -1, 0, 0}, + {2, 16, 2, 0, 0}, + {-1, 2, 8, 2, -1}, + {0, 0, 2, 16, 2}, + {0, 0, -1, 2, 4}}; + constexpr int z1[5][5] = {{7, -8, 1, 0, 0}, + {-8, 16, -8, 0, 0}, + {1, -8, 14, -8, 1}, + {0, 0, -8, 16, -8}, + {0, 0, 1, -8, 7}}; + const int nodes[2][4] = {{0, 1, 4, 3}, {1, 2, 5, 4}}; + for (std::size_t e = 0; e < 2; ++e) { + const double h = e == 0 ? 1. : 2.; + const double conductivity = e == 0 ? 2. : 10.; + const double capacity = e == 0 ? 12. : 1.; + for (std::size_t i = 0; i < 4; ++i) + for (std::size_t j = 0; j < 4; ++j) + for (std::size_t a = 0; a < 5; ++a) + for (std::size_t b = 0; b < 5; ++b) { + const auto p = 5 * nodes[e][i] + a; + const auto q = 5 * nodes[e][j] + b; + c(p, q) += capacity * h * m[i][j] * z0[a][b] / 2160.; + k(p, q) += conductivity * (h * s[i][j] * z0[a][b] / 360. + + m[i][j] * z1[a][b] / (54 * h)); + const double sx[4] = {-1, 1, 1, -1}; + const double sy[4] = {-1, -1, 1, 1}; + const double ps = (sx[i] == sx[j] ? 1. / 3 : 1. / 6) + + (sy[i] == sy[j] ? 1. / 3 : 1. / 6); + result.pk(p, q) += + conductivity * (h * ps * std::abs(z0[a][b]) / 60 + + std::abs(m[i][j] * z1[a][b] / (54 * h))); + result.pc(p, q) += + std::abs(capacity * h * m[i][j] * z0[a][b] / 2160.); + } + } + return result; +} + +void SameBits(const std::vector& a, const std::vector& b) { + ASSERT_EQ(a.size(), b.size()); + EXPECT_EQ(std::memcmp(a.data(), b.data(), a.size() * sizeof(double)), 0); +} +} // namespace + +TEST(ThermalAssembler, + HandAssemblyPreservesSharedStationsMaterialOffsetAndThicknessJump) { + auto domain = fesa::Domain::Create(Model()); + ASSERT_TRUE(domain.HasValue()); + const auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + const auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + fesa::SerialParallelFor serial; + const auto assembled = + fesa::AssembleThermal(model.Value(), dofs.Value(), serial); + ASSERT_TRUE(assembled.HasValue()); + ASSERT_EQ(assembled.Value().Conductivity().Rows(), 30U); + const auto expected = Expected(); + const auto k = Dense(assembled.Value().Conductivity()); + const auto c = Dense(assembled.Value().Capacity()); + const auto ek = Dense(assembled.Value().ConductivityError()); + const auto ec = Dense(assembled.Value().CapacityError()); + const double gamma = + 4096 * std::ldexp(1., -53) / (1 - 4096 * std::ldexp(1., -53)); + for (std::size_t i = 0; i < 30; ++i) + for (std::size_t j = 0; j < 30; ++j) { + EXPECT_LE(std::abs(k(i, j) - expected.k(i, j)), + ek(i, j) + gamma * expected.pk(i, j)); + EXPECT_LE(std::abs(c(i, j) - expected.c(i, j)), + ec(i, j) + gamma * expected.pc(i, j)); + } + EXPECT_EQ(domain.Value().ThermalElements()[0].MaterialIndex(), 1U); + EXPECT_EQ(domain.Value().ThermalElements()[1].MaterialIndex(), 2U); + EXPECT_EQ(assembled.Value().ElementOrder(), + (std::vector{0, 1})); + EXPECT_EQ(assembled.Value().Kernels()[0].SourceId().source_label_text, + "0009"); + for (std::size_t a = 0; a < 5; ++a) + EXPECT_EQ(assembled.Value().Scatters()[0][5 + a], + assembled.Value().Scatters()[1][a]); +} + +TEST(ThermalAssembler, + SourceRowHeatInputRetainsPrescribedSignedDuplicateContributions) { + auto domain = fesa::Domain::Create(Model()); + ASSERT_TRUE(domain.HasValue()); + const auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + const auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + const auto result = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::SerialParallelFor{}); + ASSERT_TRUE(result.HasValue()); + const auto full = + dofs.Value().FullThermalDof(2, fesa::ThermalDofComponent::kBottom); + EXPECT_FALSE(dofs.Value().FreeEquation(full).has_value()); + EXPECT_EQ(result.Value().HeatInput()[full], 5.); + EXPECT_EQ(result.Value().HeatInput()[4], -7.); + EXPECT_GE(result.Value().HeatInputError()[full], 1.); + for (std::size_t i = 0; i < 30; ++i) + if (i != full && i != 4) { + EXPECT_EQ(result.Value().HeatInput()[i], 0.); + EXPECT_EQ(result.Value().HeatInputError()[i], 0.); + } + auto reordered = Model(); + std::swap(reordered.thermal_step->loads[1], reordered.thermal_step->loads[2]); + auto second_domain = fesa::Domain::Create(std::move(reordered)); + ASSERT_TRUE(second_domain.HasValue()); + auto second_model = fesa::AnalysisModel::Create(second_domain.Value()); + ASSERT_TRUE(second_model.HasValue()); + auto second_dofs = fesa::DofManager::CreateThermal(second_model.Value()); + ASSERT_TRUE(second_dofs.HasValue()); + auto second = fesa::AssembleThermal(second_model.Value(), second_dofs.Value(), + fesa::SerialParallelFor{}); + ASSERT_TRUE(second.HasValue()); + EXPECT_EQ(second.Value().HeatInput()[full], 6.); +} + +TEST(ThermalAssembler, + SerialAndActualTbbOperatorsBoundsAndFluxAreBitwiseEqual) { + auto domain = fesa::Domain::Create(Model()); + ASSERT_TRUE(domain.HasValue()); + const auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + const auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + const auto serial = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::SerialParallelFor{}); + const auto parallel = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::TbbParallelFor{}); + ASSERT_TRUE(serial.HasValue()); + ASSERT_TRUE(parallel.HasValue()); + EXPECT_EQ(serial.Value().Conductivity().RowOffsets(), + parallel.Value().Conductivity().RowOffsets()); + EXPECT_EQ(serial.Value().Conductivity().ColumnIndices(), + parallel.Value().Conductivity().ColumnIndices()); + SameBits(serial.Value().Conductivity().Values(), + parallel.Value().Conductivity().Values()); + SameBits(serial.Value().Capacity().Values(), + parallel.Value().Capacity().Values()); + SameBits(serial.Value().ConductivityError().Values(), + parallel.Value().ConductivityError().Values()); + SameBits(serial.Value().CapacityError().Values(), + parallel.Value().CapacityError().Values()); + EXPECT_EQ( + std::memcmp(serial.Value().HeatInput().Data(), + parallel.Value().HeatInput().Data(), 30 * sizeof(double)), + 0); + fesa::Vector temperature{30}; + for (std::size_t i = 0; i < 30; ++i) + temperature[i] = static_cast(i) - 10; + const auto sf = + serial.Value().RecoverFlux(temperature, fesa::SerialParallelFor{}); + const auto pf = + parallel.Value().RecoverFlux(temperature, fesa::TbbParallelFor{}); + ASSERT_TRUE(sf.HasValue()); + ASSERT_TRUE(pf.HasValue()); + for (std::size_t e = 0; e < 2; ++e) + for (std::size_t r = 0; r < 24; ++r) + EXPECT_EQ( + std::memcmp(sf.Value()[e][r].Components().data(), + pf.Value()[e][r].Components().data(), 3 * sizeof(double)), + 0); +} + +TEST(ThermalAssembler, DisconnectedComponentsHaveNoCrossBlockEntries) { + auto definition = Model(); + definition.thermal_step->loads.clear(); + definition.thermal_step->boundaries.clear(); + definition.thermal_elements[1].node_indices = {{7, 8, 9, 10}}; + for (const auto point : std::array, 4>{ + {{3, 0, 0}, {4, 0, 0}, {4, 1, 0}, {3, 1, 0}}}) { + const auto id = static_cast(definition.nodes.size() + 100); + definition.nodes.push_back({{"other", id, std::to_string(id)}, point, {}}); + } + definition.thermal_elements[1].source_id.instance_name = "other"; + auto domain = fesa::Domain::Create(std::move(definition)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + auto result = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::TbbParallelFor{}); + ASSERT_TRUE(result.HasValue()); + const auto& matrix = result.Value().Capacity(); + ASSERT_EQ(matrix.Rows(), 40U); + for (std::size_t row = 0; row < 40; ++row) + for (std::size_t p = matrix.RowOffsets()[row]; + p < matrix.RowOffsets()[row + 1]; ++p) + EXPECT_EQ(row < 20, matrix.ColumnIndices()[p] < 20); +} + +TEST(ThermalAssembler, RejectsWrongProfilesPartitionsAndNonfiniteAccumulation) { + auto domain = fesa::Domain::Create(Model()); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + EXPECT_FALSE(fesa::AssembleThermal(model.Value(), fesa::DofManager{}, + fesa::SerialParallelFor{}) + .HasValue()); + auto other_definition = Model(); + other_definition.thermal_step->boundaries.clear(); + auto other_domain = fesa::Domain::Create(other_definition); + ASSERT_TRUE(other_domain.HasValue()); + auto other_model = fesa::AnalysisModel::Create(other_domain.Value()); + ASSERT_TRUE(other_model.HasValue()); + auto other_dofs = fesa::DofManager::CreateThermal(other_model.Value()); + ASSERT_TRUE(other_dofs.HasValue()); + EXPECT_FALSE(fesa::AssembleThermal(model.Value(), other_dofs.Value(), + fesa::SerialParallelFor{}) + .HasValue()); + auto definition = Model(); + definition.thermal_step->loads = { + {2, 11, std::numeric_limits::max(), {}}, + {2, 11, std::numeric_limits::max(), {}}}; + auto overflow_domain = fesa::Domain::Create(definition); + ASSERT_TRUE(overflow_domain.HasValue()); + auto overflow_model = fesa::AnalysisModel::Create(overflow_domain.Value()); + ASSERT_TRUE(overflow_model.HasValue()); + auto overflow_dofs = fesa::DofManager::CreateThermal(overflow_model.Value()); + ASSERT_TRUE(overflow_dofs.HasValue()); + EXPECT_FALSE(fesa::AssembleThermal(overflow_model.Value(), + overflow_dofs.Value(), + fesa::TbbParallelFor{}) + .HasValue()); +} + +TEST(ThermalAssembler, ParallelKernelFailuresSelectStableSourceElement) { + auto definition = Model(); + for (auto& material : definition.thermal_materials) { + material.density = std::numeric_limits::max(); + material.specific_heat = 2; + } + auto domain = fesa::Domain::Create(definition); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + const auto serial = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::SerialParallelFor{}); + const auto parallel = fesa::AssembleThermal(model.Value(), dofs.Value(), + fesa::TbbParallelFor{}); + ASSERT_FALSE(serial.HasValue()); + ASSERT_FALSE(parallel.HasValue()); + EXPECT_EQ(serial.GetStatus().Diagnostics()[0].entity_identity, "plate:0009"); + EXPECT_EQ(parallel.GetStatus().Diagnostics()[0].entity_identity, + "plate:0009"); +}