baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
176d53cb67d8d2453cd0ce2ff9f9eead18fc6de1
FESADev/src/fesa/io/hdf5/eigen_hdf5_results_writer.cpp
T
KOKO\Mimi 7ebdce373c fix: use residual stopping for FEAST eigenanalysis
2026-09-10 01:52:49 +09:00

975 lines
45 KiB
C++
Raw Blame History

#include "fesa/io/hdf5/eigen_hdf5_results_writer.h"
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <exception>
#include <limits>
#include <set>
#include <string>
#include <utility>
#include <vector>
#include "fesa/analysis/analysis_model.h"
#include "fesa/build_info.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/model/domain.h"
#include "io/hdf5/eigen_hdf5_self_check.h"
#include "io/hdf5/hdf5_atomic_file.h"
#include "io/hdf5/hdf5_model_writer.h"
#include "io/hdf5/hdf5_primitives.h"
#include "io/hdf5/hdf5_raii.h"
namespace fesa {
namespace {
constexpr const char* kEigenRoot = "/steps/Step-1/eigen";
void WriteQuantityAttributes(hid_t dataset, const std::string& components,
const std::string& units,
const std::string& coordinates,
const std::string& location) {
hdf5_internal::WriteStringAttribute(dataset, "component_names", components);
hdf5_internal::WriteStringAttribute(dataset, "component_unit_dimensions",
units);
hdf5_internal::WriteStringAttribute(dataset, "coordinate_system",
coordinates);
hdf5_internal::WriteStringAttribute(dataset, "location", location);
}
void WriteModalAttributes(hid_t dataset, const std::string& source_step) {
hdf5_internal::WriteUint64Attribute(dataset, "eigen_schema_version", 1U);
hdf5_internal::WriteStringAttribute(dataset, "source_step_name", source_step);
hdf5_internal::WriteStringAttribute(dataset, "domain", "mode");
}
bool IsFinite(const Matrix& matrix) {
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
if (!std::isfinite(matrix(row, column))) {
return false;
}
}
}
return true;
}
Status Validate(const Domain& domain, const EigenAnalysisState& state,
const std::vector<Diagnostic>& diagnostics) {
const FrequencyStepDefinition* frequency = domain.Frequency();
if (frequency == nullptr || !hdf5_internal::IsShellDomain(domain)) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Eigen HDF5 output requires one shell frequency domain.");
}
std::set<EntityIndex> density_materials;
EntityIndex previous_density_material = 0U;
bool first_density = true;
for (const MaterialDensityDefinition& density : domain.Densities()) {
if (density.material_index >= domain.Materials().Size() ||
!std::isfinite(density.rho) || !(density.rho > 0.0) ||
(!first_density &&
density.material_index <= previous_density_material)) {
return hdf5_internal::OutputFailure("invalid-eigen-result-state",
"Eigen material density inventory "
"must be finite and source ordered.");
}
density_materials.insert(density.material_index);
previous_density_material = density.material_index;
first_density = false;
}
for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) {
if (density_materials.count(
domain.ShellElements()[index].MaterialIndex()) == 0U) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Every active shell material requires one positive density.");
}
}
const std::size_t node_count = domain.Nodes().size();
if (node_count > (std::numeric_limits<std::size_t>::max)() / 6U ||
state.full_dof_count != 6U * node_count ||
state.free_dof_count !=
state.mass_null_dimension + state.finite_dimension ||
state.modes.size() !=
static_cast<std::size_t>(frequency->requested_mode_count) ||
state.vectors.Rows() != state.modes.size() ||
state.vectors.Columns() != 6U * node_count ||
state.nodal_reduction.size() != node_count ||
state.constraint_mask.size() != 6U * node_count ||
state.prescribed_displacement.size() != 6U * node_count ||
state.eigenpair_checks.size() < state.modes.size() ||
state.eigenpair_checks.size() > state.finite_dimension ||
!IsFinite(state.vectors) || !std::isfinite(state.characteristic_length) ||
!(state.characteristic_length > 0.0) ||
!std::isfinite(state.backend.operator_scale) ||
!(state.backend.operator_scale > 0.0)) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Eigen state does not satisfy the durable result inventory.");
}
const auto metric_ok = [](double value, double limit) {
return std::isfinite(value) && value >= 0.0 && value <= limit;
};
const auto& backend = state.backend;
if (backend.base_fpm.size() != 128U || backend.backend_calls.empty() ||
backend.backend_calls.size() > 512U || backend.thread_count != 1U ||
backend.seed != 0x4D495443344D4645ULL ||
backend.runtime_version.empty() ||
(backend.integer_abi != "LP64" && backend.integer_abi != "ILP64") ||
backend.requested_modes != state.modes.size() ||
!metric_ok(backend.gram_error, 1e-10) ||
!std::isfinite(backend.enclosure_lower) ||
!std::isfinite(backend.enclosure_upper) ||
!std::isfinite(backend.solve_lower) ||
!std::isfinite(backend.solve_upper) ||
!(backend.enclosure_lower < backend.enclosure_upper) ||
backend.solve_lower != backend.enclosure_lower ||
!(backend.solve_lower < backend.solve_upper) ||
backend.solve_upper > backend.enclosure_upper) {
return hdf5_internal::OutputFailure("invalid-eigen-result-state",
"Backend inventory is invalid.");
}
for (const auto& control :
std::array<std::pair<std::size_t, std::int64_t>, 11>{{{0U, 0},
{1U, 16},
{2U, 12},
{3U, 20},
{4U, 1},
{5U, 0},
{13U, 0},
{26U, 1},
{27U, 1},
{49U, 0},
{63U, 0}}}) {
if (backend.base_fpm[control.first] != control.second)
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Backend controls differ from the numerical contract.");
}
std::size_t solve_calls = 0U;
bool enclosed = false;
for (std::size_t i = 0U; i < backend.backend_calls.size(); ++i) {
const auto& call = backend.backend_calls[i];
const bool count = call.kind == "count", solve = call.kind == "solve";
if (solve) ++solve_calls;
if (call.ordinal != i || (!count && !solve) || !std::isfinite(call.lower) ||
!std::isfinite(call.upper) || !(call.lower < call.upper) ||
call.allocated_m0 == 0U || call.allocated_m0 > state.finite_dimension ||
(count && call.allocated_m0 != state.finite_dimension) ||
(call.count_available && call.count > call.allocated_m0) ||
(!call.count_available && call.count != 0U) ||
(call.loop_available && call.loop < 0) ||
(!call.loop_available && call.loop != 0) ||
(call.epsout_available && !std::isfinite(call.epsout)) ||
(!call.epsout_available && call.epsout != 0.0) ||
(call.outcome != "accepted" && call.outcome != "retry" &&
call.outcome != "optional-abandoned") ||
(call.outcome == "accepted" &&
(call.info != 0 || !call.count_available)) ||
call.decision_reason.empty() ||
(enclosed && call.lower != backend.enclosure_lower)) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Backend call history is inconsistent.");
}
if (call.purpose == "total-enclosure" && call.outcome == "accepted" &&
call.count == state.finite_dimension)
enclosed = true;
}
const auto& last_call = backend.backend_calls.back();
if (!enclosed || solve_calls == 0U || solve_calls > 3U ||
last_call.kind != "solve" || last_call.outcome != "accepted" ||
last_call.count != state.eigenpair_checks.size() ||
!last_call.loop_available || !last_call.epsout_available ||
last_call.lower != backend.solve_lower ||
last_call.upper != backend.solve_upper) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state", "Final all-c solve evidence is invalid.");
}
if (!metric_ok(state.eta_gram, 1e-10) ||
!metric_ok(state.max_eta_orth, 1e-10) ||
!metric_ok(state.eta_mz, 1e-12) || !metric_ok(state.eta_kz, 1e-12) ||
!metric_ok(state.eta_pz, 1e-12) ||
state.null_checks_applicable != (state.mass_null_dimension != 0U) ||
(!state.null_checks_applicable &&
(state.eta_mz != 0.0 || state.eta_kz != 0.0 || state.eta_pz != 0.0))) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state", "Global modal verification is invalid.");
}
const auto model = AnalysisModel::Create(domain);
if (!model.HasValue()) return model.GetStatus();
const auto dofs = DofManager::Create(model.Value());
if (!dofs.HasValue()) return dofs.GetStatus();
for (std::size_t i = 0U; i < state.full_dof_count; ++i) {
const auto constrained =
static_cast<std::uint8_t>(!dofs.Value().FreeEquation(i).has_value());
if (state.constraint_mask[i] != constrained)
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state", "Constraint mask differs from Domain.");
if (constrained != 0U) {
for (std::size_t mode = 0U; mode < state.vectors.Rows(); ++mode)
if (state.vectors(mode, i) != 0.0)
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Constrained modal component is nonzero.");
}
}
std::size_t cluster = 1U;
for (std::size_t index = 0U; index < state.modes.size(); ++index) {
const EigenMode& mode = state.modes[index];
const double mu = mode.eigenvalue / state.backend.operator_scale;
if (index > 0U &&
mu - state.modes[index - 1U].eigenvalue / state.backend.operator_scale >
1e-8)
++cluster;
if (mode.mode_id != index + 1U || !std::isfinite(mode.eigenvalue) ||
!std::isfinite(mu) || mu < -1e-9 || mode.cluster_id != cluster ||
mode.classification != static_cast<std::uint8_t>(std::abs(mu) > 1e-9) ||
(index > 0U && mode.eigenvalue < state.modes[index - 1U].eigenvalue) ||
(mode.eigenvalue < 0.0 && mode.frequency.has_value()) ||
(mode.eigenvalue >= 0.0 && !mode.frequency.has_value()) ||
(mode.frequency &&
(!std::isfinite(*mode.frequency) || *mode.frequency < 0.0 ||
(mode.eigenvalue > 0.0 && *mode.frequency == 0.0) ||
*mode.frequency !=
(mode.eigenvalue == 0.0
? 0.0
: std::sqrt(mode.eigenvalue) / (2.0 * std::acos(-1.0)))))) {
return hdf5_internal::OutputFailure("invalid-eigen-result-state",
"Eigen mode inventory is invalid.");
}
}
std::size_t free_count = 0U;
std::size_t null_count = 0U;
std::size_t positive_count = 0U;
for (const EigenNodalReductionRecord& record : state.nodal_reduction) {
if (record.free_count != record.null_count + record.positive_count ||
record.free_count > 6U || record.null_count > 1U ||
record.rcond_available != (record.positive_count != 0U) ||
(record.rcond_available &&
(!std::isfinite(record.rcond_inf) ||
record.rcond_inf < std::ldexp(1.0, -32) || record.rcond_inf > 1.0)) ||
(!record.rcond_available && record.rcond_inf != 0.0)) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Nodal mass-reduction evidence is invalid.");
}
free_count += record.free_count;
null_count += record.null_count;
positive_count += record.positive_count;
}
if (free_count != state.free_dof_count ||
null_count != state.mass_null_dimension ||
positive_count != state.finite_dimension) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Nodal mass-reduction counts do not match state dimensions.");
}
std::set<std::size_t> columns;
for (std::size_t index = 0U; index < state.eigenpair_checks.size(); ++index) {
const EigenpairCheck& check = state.eigenpair_checks[index];
if (check.sorted_rank != index + 1U || !std::isfinite(check.eigenvalue) ||
check.solver_column >= state.eigenpair_checks.size() ||
!columns.insert(check.solver_column).second ||
(index < state.modes.size() &&
check.eigenvalue != state.modes[index].eigenvalue) ||
(index > 0U &&
check.eigenvalue < state.eigenpair_checks[index - 1U].eigenvalue) ||
!std::isfinite(check.modal_mass) ||
std::abs(check.modal_mass - 1.0) > 1e-10 ||
!std::isfinite(check.rayleigh_quotient) ||
!std::isfinite(check.eta_r) || !std::isfinite(check.eta_w) ||
!std::isfinite(check.rho_b) || !std::isfinite(check.eta_norm) ||
!metric_ok(check.eta_r, 1e-10) || !metric_ok(check.eta_w, 1e-10) ||
!metric_ok(check.rho_b, 1e-10) || !metric_ok(check.eta_norm, 1e-10) ||
!metric_ok(check.eta_energy, 1e-10) ||
(check.backend_res_available && !std::isfinite(check.backend_res)) ||
(!check.backend_res_available && check.backend_res != 0.0)) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Eigenpair verification evidence is invalid.");
}
}
for (std::uint8_t value : state.constraint_mask) {
if (value > 1U) {
return hdf5_internal::OutputFailure("invalid-eigen-result-state",
"Constraint mask must be binary.");
}
}
for (double value : state.prescribed_displacement) {
if (value != 0.0) {
return hdf5_internal::OutputFailure(
"invalid-eigen-result-state",
"Eigen prescribed displacements must be exact zero.");
}
}
for (const Diagnostic& diagnostic : diagnostics) {
if (diagnostic.severity != Severity::kWarning) {
return hdf5_internal::OutputFailure(
"invalid-eigen-diagnostics",
"Successful eigen output stores warnings only.");
}
}
return Status::Ok();
}
void ReplaceStringAttribute(hid_t object, const char* name,
const std::string& value) {
if (H5Aexists(object, name) > 0) {
hdf5_internal::RequireHdf5(H5Adelete(object, name),
"Unable to replace HDF5 metadata attribute.");
}
hdf5_internal::WriteStringAttribute(object, name, value);
}
void WriteEigenMetadata(hid_t file) {
hdf5_internal::Hdf5Handle metadata{
hdf5_internal::RequireHdf5Id(H5Gopen2(file, "/metadata", H5P_DEFAULT),
"Unable to open metadata."),
H5Gclose};
ReplaceStringAttribute(metadata.Get(), "feature_id",
"mitc4-shell-feast-eigen2");
ReplaceStringAttribute(metadata.Get(), "analysis_procedure",
"linear-frequency");
ReplaceStringAttribute(metadata.Get(), "coordinate_convention",
"global-cartesian; shell-local=initial-reference; "
"positive-thickness=+zeta");
hdf5_internal::WriteUint64Attribute(metadata.Get(), "eigen_schema_version",
1U);
hdf5_internal::WriteStringAttribute(metadata.Get(), "mass_integration_rule",
"3x3x3-gauss");
hdf5_internal::WriteStringAttribute(metadata.Get(), "stabilization_rule",
"initial-physical-rotary-trace-mean");
hdf5_internal::WriteStringAttribute(metadata.Get(), "mass_variant",
"derived-lumped");
hdf5_internal::WriteStringAttribute(metadata.Get(), "mass_origin",
"consistent-moment-row-sum");
hdf5_internal::WriteStringAttribute(metadata.Get(), "normalization", "mass");
hdf5_internal::WriteStringAttribute(metadata.Get(), "eigensolver_backend",
"MKL FEAST");
hdf5_internal::WriteStringAttribute(metadata.Get(), "eigensolver_routine",
"dfeast_scsrgv");
hdf5_internal::WriteStringAttribute(
metadata.Get(), "backend_problem",
"scaled-whitened-lumped-with-identity-mass");
hdf5_internal::WriteStringAttribute(metadata.Get(), "numerical_contract",
"mitc4-shell-feast-eigen2:r1");
}
void WriteStep(hid_t file, const Domain& domain,
const EigenAnalysisState& state) {
const FrequencyStepDefinition& source = *domain.Frequency();
auto step = hdf5_internal::CreateGroup(file, "/steps/Step-1");
hdf5_internal::WriteStringAttribute(step.Get(), "source_step_name",
source.source_step_name_raw);
hdf5_internal::WriteStringAttribute(step.Get(), "source_step_name_canonical",
source.source_step_name_canonical);
hdf5_internal::WriteUint64Attribute(step.Get(), "step_source_line",
source.step_location.line);
hdf5_internal::WriteUint64Attribute(step.Get(), "frequency_source_line",
source.frequency_location.line);
hdf5_internal::WriteUint64Attribute(step.Get(), "frequency_data_line",
source.data_location.line);
hdf5_internal::WriteUint8Attribute(step.Get(), "source_nlgeom", 0U);
hdf5_internal::WriteUint8Attribute(
step.Get(), "source_perturbation_present",
static_cast<std::uint8_t>(source.perturbation_explicit));
const std::array<std::uint8_t, 4> presence{
static_cast<std::uint8_t>(source.eigensolver_explicit),
static_cast<std::uint8_t>(source.sim_explicit),
static_cast<std::uint8_t>(source.acoustic_coupling_explicit),
static_cast<std::uint8_t>(source.normalization_explicit)};
auto presence_space = hdf5_internal::CreateDatasetSpace({4U});
hdf5_internal::Hdf5Handle presence_attr{
hdf5_internal::RequireHdf5Id(
H5Acreate2(step.Get(), "source_option_presence", H5T_STD_U8LE,
presence_space.Get(), H5P_DEFAULT, H5P_DEFAULT),
"Unable to create source option presence."),
H5Aclose};
hdf5_internal::RequireHdf5(
H5Awrite(presence_attr.Get(), H5T_NATIVE_UINT8, presence.data()),
"Unable to write source option presence.");
hdf5_internal::WriteStringAttribute(step.Get(), "source_eigensolver",
"LANCZOS");
hdf5_internal::WriteStringAttribute(step.Get(), "source_sim", "YES");
hdf5_internal::WriteStringAttribute(step.Get(), "source_acoustic_coupling",
"ON");
hdf5_internal::WriteStringAttribute(step.Get(), "source_normalization",
"MASS");
hdf5_internal::WriteUint64Attribute(step.Get(), "requested_mode_count",
state.modes.size());
hdf5_internal::WriteUint64Attribute(step.Get(), "returned_mode_count",
state.modes.size());
hdf5_internal::WriteUint64Attribute(step.Get(), "full_dof_count",
state.full_dof_count);
hdf5_internal::WriteUint64Attribute(step.Get(), "free_dof_count",
state.free_dof_count);
hdf5_internal::WriteUint64Attribute(step.Get(), "mass_null_dimension",
state.mass_null_dimension);
hdf5_internal::WriteUint64Attribute(step.Get(), "finite_dimension",
state.finite_dimension);
hdf5_internal::WriteUint64Attribute(step.Get(), "validated_eigenpair_count",
state.eigenpair_checks.size());
hdf5_internal::WriteStringAttribute(step.Get(), "frame_indexing",
"mode-rank-one-based");
}
struct ModeRow {
std::uint64_t mode_id;
double eigenvalue;
std::uint8_t classification;
std::uint64_t cluster_id;
};
struct FrequencyRow {
std::uint64_t mode_id;
double frequency;
};
struct BackendCallRow {
std::uint64_t ordinal;
const char* kind;
const char* purpose;
double lower;
double upper;
std::uint64_t allocated_m0;
std::int64_t info;
std::uint8_t count_available;
std::uint64_t count;
std::uint8_t loop_available;
std::uint64_t loop;
std::uint8_t epsout_available;
double epsout;
const char* outcome;
const char* decision_reason;
};
struct DensityRow {
std::uint64_t material_internal_id;
double rho;
const char* source_file;
std::uint64_t source_line;
};
struct NodalReductionRow {
std::uint8_t free_count;
std::uint8_t null_count;
std::uint8_t positive_count;
std::uint8_t rcond_available;
double rcond_inf;
};
struct EigenpairCheckRow {
std::uint64_t sorted_rank;
std::uint64_t solver_column;
double eigenvalue;
double modal_mass;
double rayleigh_quotient;
double eta_r;
double eta_w;
double rho_b;
double eta_norm;
double eta_energy;
std::uint8_t backend_res_available;
double backend_res;
};
struct DiagnosticRow {
const char* severity;
const char* code;
const char* file;
std::uint64_t line;
const char* keyword;
const char* entity_identity;
const char* message;
};
template <typename Row>
hdf5_internal::Hdf5Handle MakeCompoundType() {
return hdf5_internal::Hdf5Handle{
hdf5_internal::RequireHdf5Id(H5Tcreate(H5T_COMPOUND, sizeof(Row)),
"Unable to create compound type."),
H5Tclose};
}
void WriteModes(hid_t file, const EigenAnalysisState& state,
const std::string& source_step) {
std::vector<ModeRow> modes;
std::vector<FrequencyRow> frequencies;
modes.reserve(state.modes.size());
for (const EigenMode& mode : state.modes) {
modes.push_back(
{mode.mode_id, mode.eigenvalue, mode.classification, mode.cluster_id});
if (mode.frequency) frequencies.push_back({mode.mode_id, *mode.frequency});
}
auto mode_file = MakeCompoundType<ModeRow>();
auto mode_memory = MakeCompoundType<ModeRow>();
for (hid_t type : {mode_file.Get(), mode_memory.Get()}) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "mode_id", HOFFSET(ModeRow, mode_id), H5T_STD_U64LE),
"Unable to define mode ID.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "eigenvalue", HOFFSET(ModeRow, eigenvalue),
H5T_IEEE_F64LE),
"Unable to define eigenvalue.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "classification", HOFFSET(ModeRow, classification),
H5T_STD_U8LE),
"Unable to define classification.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "cluster_id", HOFFSET(ModeRow, cluster_id),
H5T_STD_U64LE),
"Unable to define cluster ID.");
}
auto dataset = hdf5_internal::WriteCompoundDataset(
file, "/steps/Step-1/eigen/modes", modes.size(), mode_file.Get(),
mode_memory.Get(), modes.data());
WriteQuantityAttributes(dataset.Get(), "Value", "time^-2", "scalar",
"modal-frame");
WriteModalAttributes(dataset.Get(), source_step);
hdf5_internal::WriteStringAttribute(dataset.Get(), "quantity", "eigenvalue");
auto frequency_file = MakeCompoundType<FrequencyRow>();
auto frequency_memory = MakeCompoundType<FrequencyRow>();
for (hid_t type : {frequency_file.Get(), frequency_memory.Get()}) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "mode_id", HOFFSET(FrequencyRow, mode_id),
H5T_STD_U64LE),
"Unable to define frequency mode ID.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "frequency", HOFFSET(FrequencyRow, frequency),
H5T_IEEE_F64LE),
"Unable to define frequency.");
}
auto frequency_dataset = hdf5_internal::WriteCompoundDataset(
file, "/steps/Step-1/eigen/frequencies", frequencies.size(),
frequency_file.Get(), frequency_memory.Get(), frequencies.data());
WriteQuantityAttributes(frequency_dataset.Get(), "Freq", "time^-1", "scalar",
"modal-frame");
WriteModalAttributes(frequency_dataset.Get(), source_step);
hdf5_internal::WriteStringAttribute(frequency_dataset.Get(), "quantity",
"cycle-frequency");
hdf5_internal::WriteStringAttribute(
frequency_dataset.Get(), "unavailable_rule", "negative-raw-eigenvalue");
}
void WriteBaseFpm(hid_t file, const std::vector<std::int64_t>& values) {
const std::vector<hsize_t> dimensions{static_cast<hsize_t>(values.size())};
auto space = hdf5_internal::CreateDatasetSpace(dimensions);
hdf5_internal::Hdf5Handle dataset{
hdf5_internal::RequireHdf5Id(
H5Dcreate2(file, "/steps/Step-1/eigen/backend/base_fpm",
H5T_STD_I64LE, space.Get(), H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT),
"Unable to create eigen backend control dataset."),
H5Dclose};
if (!values.empty()) {
hdf5_internal::RequireHdf5(
H5Dwrite(dataset.Get(), H5T_NATIVE_INT64, H5S_ALL, H5S_ALL, H5P_DEFAULT,
values.data()),
"Unable to write eigen backend control dataset.");
}
WriteQuantityAttributes(dataset.Get(), "fpm", "1", "none",
"backend-parameter");
}
void WriteBackendCalls(hid_t file, const std::vector<FeastBackendCall>& calls) {
std::vector<BackendCallRow> rows;
rows.reserve(calls.size());
for (const FeastBackendCall& call : calls) {
rows.push_back({call.ordinal, call.kind.c_str(), call.purpose.c_str(),
call.lower, call.upper, call.allocated_m0, call.info,
static_cast<std::uint8_t>(call.count_available), call.count,
static_cast<std::uint8_t>(call.loop_available),
static_cast<std::uint64_t>(call.loop),
static_cast<std::uint8_t>(call.epsout_available),
call.epsout, call.outcome.c_str(),
call.decision_reason.c_str()});
}
auto string_type = hdf5_internal::MakeUtf8StringType();
auto file_type = MakeCompoundType<BackendCallRow>();
auto memory_type = MakeCompoundType<BackendCallRow>();
const auto add_fields = [&](hid_t type, hid_t u64, hid_t i64, hid_t f64) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "ordinal", HOFFSET(BackendCallRow, ordinal), u64),
"Unable to define backend ordinal.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "kind", HOFFSET(BackendCallRow, kind),
string_type.Get()),
"Unable to define backend kind.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "purpose", HOFFSET(BackendCallRow, purpose),
string_type.Get()),
"Unable to define backend purpose.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "lower", HOFFSET(BackendCallRow, lower), f64),
"Unable to define backend lower.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "upper", HOFFSET(BackendCallRow, upper), f64),
"Unable to define backend upper.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "allocated_m0", HOFFSET(BackendCallRow, allocated_m0),
u64),
"Unable to define backend capacity.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "info", HOFFSET(BackendCallRow, info), i64),
"Unable to define backend info.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "count_available",
HOFFSET(BackendCallRow, count_available), H5T_STD_U8LE),
"Unable to define count flag.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "count", HOFFSET(BackendCallRow, count), u64),
"Unable to define count.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "loop_available",
HOFFSET(BackendCallRow, loop_available), H5T_STD_U8LE),
"Unable to define loop flag.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "loop", HOFFSET(BackendCallRow, loop), u64),
"Unable to define loop.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "epsout_available",
HOFFSET(BackendCallRow, epsout_available), H5T_STD_U8LE),
"Unable to define epsout flag.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "epsout", HOFFSET(BackendCallRow, epsout), f64),
"Unable to define epsout.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "outcome", HOFFSET(BackendCallRow, outcome),
string_type.Get()),
"Unable to define outcome.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "decision_reason",
HOFFSET(BackendCallRow, decision_reason), string_type.Get()),
"Unable to define decision reason.");
};
add_fields(file_type.Get(), H5T_STD_U64LE, H5T_STD_I64LE, H5T_IEEE_F64LE);
add_fields(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_INT64,
H5T_NATIVE_DOUBLE);
auto dataset = hdf5_internal::WriteCompoundDataset(
file, "/steps/Step-1/eigen/backend/calls", rows.size(), file_type.Get(),
memory_type.Get(), rows.data());
WriteQuantityAttributes(dataset.Get(), "lower,upper,epsout", "1,1,1",
"scalar", "backend-call");
}
void WriteBackendAttributes(hid_t file, const EigenAnalysisState& state) {
hdf5_internal::Hdf5Handle backend{
hdf5_internal::RequireHdf5Id(
H5Gopen2(file, "/steps/Step-1/eigen/backend", H5P_DEFAULT),
"Unable to open eigen backend group."),
H5Gclose};
const FeastEigenResult& result = state.backend;
hdf5_internal::WriteStringAttribute(backend.Get(), "runtime_version",
result.runtime_version);
hdf5_internal::WriteStringAttribute(backend.Get(), "integer_abi",
result.integer_abi);
hdf5_internal::WriteUint64Attribute(backend.Get(), "thread_count",
result.thread_count);
hdf5_internal::WriteUint64Attribute(backend.Get(), "seed", result.seed);
hdf5_internal::WriteDoubleAttribute(backend.Get(), "enclosure_lower",
result.enclosure_lower);
hdf5_internal::WriteDoubleAttribute(backend.Get(), "enclosure_upper",
result.enclosure_upper);
hdf5_internal::WriteDoubleAttribute(backend.Get(), "solve_lower",
result.solve_lower);
hdf5_internal::WriteDoubleAttribute(backend.Get(), "solve_upper",
result.solve_upper);
std::size_t count_calls = 0U;
std::size_t solve_calls = 0U;
std::size_t final_solve_ordinal = 0U;
for (const FeastBackendCall& call : result.backend_calls) {
if (call.kind == "count") ++count_calls;
if (call.kind == "solve") {
++solve_calls;
final_solve_ordinal = call.ordinal;
}
}
hdf5_internal::WriteUint64Attribute(backend.Get(), "count_calls",
count_calls);
hdf5_internal::WriteUint64Attribute(backend.Get(), "solve_calls",
solve_calls);
hdf5_internal::WriteUint64Attribute(backend.Get(), "final_solve_ordinal",
final_solve_ordinal);
hdf5_internal::WriteUint64Attribute(backend.Get(), "final_count",
state.eigenpair_checks.size());
hdf5_internal::WriteStringAttribute(backend.Get(), "upper_search_termination",
result.upper_search_termination);
}
void WriteDensities(hid_t file, const Domain& domain) {
std::vector<DensityRow> rows;
std::vector<std::string> source_paths;
rows.reserve(domain.Densities().size());
source_paths.reserve(domain.Densities().size());
for (const MaterialDensityDefinition& density : domain.Densities()) {
source_paths.push_back(
hdf5_internal::NormalizedPathString(density.location.file));
rows.push_back({density.material_index, density.rho,
source_paths.back().c_str(), density.location.line});
}
auto string_type = hdf5_internal::MakeUtf8StringType();
auto file_type = MakeCompoundType<DensityRow>();
auto memory_type = MakeCompoundType<DensityRow>();
const auto add_fields = [&](hid_t type, hid_t u64, hid_t f64) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "material_internal_id",
HOFFSET(DensityRow, material_internal_id), u64),
"Unable to define density material ID.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "rho", HOFFSET(DensityRow, rho), f64),
"Unable to define density value.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "source_file", HOFFSET(DensityRow, source_file),
string_type.Get()),
"Unable to define density source file.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "source_line", HOFFSET(DensityRow, source_line), u64),
"Unable to define density source line.");
};
add_fields(file_type.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE);
add_fields(memory_type.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE);
auto dataset = hdf5_internal::WriteCompoundDataset(
file, "/model/shell/densities", rows.size(), file_type.Get(),
memory_type.Get(), rows.data());
WriteQuantityAttributes(dataset.Get(), "rho", "mass/length^3", "scalar",
"model-material");
}
void WriteVerification(hid_t file, const EigenAnalysisState& state) {
auto verification =
hdf5_internal::CreateGroup(file, "/steps/Step-1/eigen/verification");
hdf5_internal::WriteDoubleAttribute(
verification.Get(), "characteristic_length", state.characteristic_length);
hdf5_internal::WriteDoubleAttribute(verification.Get(), "operator_scale",
state.backend.operator_scale);
hdf5_internal::WriteStringAttribute(verification.Get(), "validation_scope",
"all-computed-c-columns");
hdf5_internal::WriteDoubleAttribute(verification.Get(), "eta_gram",
state.eta_gram);
hdf5_internal::WriteDoubleAttribute(verification.Get(), "max_eta_orth",
state.max_eta_orth);
hdf5_internal::WriteUint8Attribute(
verification.Get(), "null_checks_applicable",
static_cast<std::uint8_t>(state.null_checks_applicable));
hdf5_internal::WriteDoubleAttribute(verification.Get(), "eta_mz",
state.eta_mz);
hdf5_internal::WriteDoubleAttribute(verification.Get(), "eta_kz",
state.eta_kz);
hdf5_internal::WriteDoubleAttribute(verification.Get(), "eta_pz",
state.eta_pz);
std::vector<NodalReductionRow> reductions;
reductions.reserve(state.nodal_reduction.size());
for (const EigenNodalReductionRecord& record : state.nodal_reduction) {
reductions.push_back(
{record.free_count, record.null_count, record.positive_count,
static_cast<std::uint8_t>(record.rcond_available), record.rcond_inf});
}
auto reduction_file = MakeCompoundType<NodalReductionRow>();
auto reduction_memory = MakeCompoundType<NodalReductionRow>();
const auto add_reduction = [](hid_t type, hid_t f64) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "free_count", HOFFSET(NodalReductionRow, free_count),
H5T_STD_U8LE),
"Unable to define free count.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "null_count", HOFFSET(NodalReductionRow, null_count),
H5T_STD_U8LE),
"Unable to define null count.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "positive_count",
HOFFSET(NodalReductionRow, positive_count), H5T_STD_U8LE),
"Unable to define positive count.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "rcond_available",
HOFFSET(NodalReductionRow, rcond_available), H5T_STD_U8LE),
"Unable to define rcond availability.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "rcond_inf", HOFFSET(NodalReductionRow, rcond_inf),
f64),
"Unable to define rcond.");
};
add_reduction(reduction_file.Get(), H5T_IEEE_F64LE);
add_reduction(reduction_memory.Get(), H5T_NATIVE_DOUBLE);
auto reduction_dataset = hdf5_internal::WriteCompoundDataset(
file, "/steps/Step-1/eigen/verification/nodal_reduction",
reductions.size(), reduction_file.Get(), reduction_memory.Get(),
reductions.data());
WriteQuantityAttributes(reduction_dataset.Get(), "rcond_inf", "1", "scalar",
"nodal-mass-reduction");
std::vector<EigenpairCheckRow> checks;
checks.reserve(state.eigenpair_checks.size());
for (const EigenpairCheck& check : state.eigenpair_checks) {
checks.push_back({check.sorted_rank, check.solver_column, check.eigenvalue,
check.modal_mass, check.rayleigh_quotient, check.eta_r,
check.eta_w, check.rho_b, check.eta_norm,
check.eta_energy,
static_cast<std::uint8_t>(check.backend_res_available),
check.backend_res});
}
auto check_file = MakeCompoundType<EigenpairCheckRow>();
auto check_memory = MakeCompoundType<EigenpairCheckRow>();
const auto add_check = [](hid_t type, hid_t u64, hid_t f64) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, "sorted_rank", HOFFSET(EigenpairCheckRow, sorted_rank),
u64),
"Unable to define sorted rank.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "solver_column",
HOFFSET(EigenpairCheckRow, solver_column), u64),
"Unable to define solver column.");
for (const auto& field : std::array<std::pair<const char*, std::size_t>, 8>{
{{"eigenvalue", HOFFSET(EigenpairCheckRow, eigenvalue)},
{"modal_mass", HOFFSET(EigenpairCheckRow, modal_mass)},
{"rayleigh_quotient",
HOFFSET(EigenpairCheckRow, rayleigh_quotient)},
{"eta_r", HOFFSET(EigenpairCheckRow, eta_r)},
{"eta_w", HOFFSET(EigenpairCheckRow, eta_w)},
{"rho_b", HOFFSET(EigenpairCheckRow, rho_b)},
{"eta_norm", HOFFSET(EigenpairCheckRow, eta_norm)},
{"eta_energy", HOFFSET(EigenpairCheckRow, eta_energy)}}}) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, field.first, field.second, f64),
"Unable to define eigenpair metric.");
}
hdf5_internal::RequireHdf5(
H5Tinsert(type, "backend_res_available",
HOFFSET(EigenpairCheckRow, backend_res_available),
H5T_STD_U8LE),
"Unable to define backend residual flag.");
hdf5_internal::RequireHdf5(
H5Tinsert(type, "backend_res", HOFFSET(EigenpairCheckRow, backend_res),
f64),
"Unable to define backend residual.");
};
add_check(check_file.Get(), H5T_STD_U64LE, H5T_IEEE_F64LE);
add_check(check_memory.Get(), H5T_NATIVE_UINT64, H5T_NATIVE_DOUBLE);
auto check_dataset = hdf5_internal::WriteCompoundDataset(
file, "/steps/Step-1/eigen/verification/eigenpair_checks", checks.size(),
check_file.Get(), check_memory.Get(), checks.data());
WriteQuantityAttributes(check_dataset.Get(),
"eigenvalue,modal_mass,rayleigh_quotient,eta_r,"
"eta_w,rho_b,eta_norm,eta_energy,backend_res",
"time^-2,1,time^-2,1,1,1,1,1,1", "scalar",
"eigenpair-verification");
}
void WriteDiagnostics(hid_t file, const std::vector<Diagnostic>& diagnostics) {
std::vector<std::string> files;
files.reserve(diagnostics.size());
std::vector<DiagnosticRow> rows;
rows.reserve(diagnostics.size());
for (const Diagnostic& diagnostic : diagnostics) {
files.push_back(
hdf5_internal::NormalizedPathString(diagnostic.location.file));
rows.push_back({"warning", diagnostic.code.c_str(), files.back().c_str(),
diagnostic.location.line, diagnostic.keyword.c_str(),
diagnostic.entity_identity.c_str(),
diagnostic.message.c_str()});
}
auto string_type = hdf5_internal::MakeUtf8StringType();
auto file_type = MakeCompoundType<DiagnosticRow>();
auto memory_type = MakeCompoundType<DiagnosticRow>();
const auto add_fields = [&](hid_t type, hid_t u64) {
for (const auto& field : std::array<std::pair<const char*, std::size_t>, 6>{
{{"severity", HOFFSET(DiagnosticRow, severity)},
{"code", HOFFSET(DiagnosticRow, code)},
{"file", HOFFSET(DiagnosticRow, file)},
{"keyword", HOFFSET(DiagnosticRow, keyword)},
{"entity_identity", HOFFSET(DiagnosticRow, entity_identity)},
{"message", HOFFSET(DiagnosticRow, message)}}}) {
hdf5_internal::RequireHdf5(
H5Tinsert(type, field.first, field.second, string_type.Get()),
"Unable to define diagnostic text field.");
}
hdf5_internal::RequireHdf5(
H5Tinsert(type, "line", HOFFSET(DiagnosticRow, line), u64),
"Unable to define diagnostic line.");
};
add_fields(file_type.Get(), H5T_STD_U64LE);
add_fields(memory_type.Get(), H5T_NATIVE_UINT64);
(void)hdf5_internal::WriteCompoundDataset(file, "/diagnostics", rows.size(),
file_type.Get(), memory_type.Get(),
rows.data());
}
} // namespace
Status EigenHdf5ResultsWriter::Write(
const std::filesystem::path& output_path, const Domain& domain,
const EigenAnalysisState& state,
const std::vector<Diagnostic>& diagnostics) {
const Status valid = Validate(domain, state, diagnostics);
if (!valid.IsOk()) return valid;
try {
hdf5_internal::Hdf5ErrorSilencer silence;
const auto candidate_path = hdf5_internal::MakeCandidatePath(output_path);
hdf5_internal::CandidateFileGuard cleanup{candidate_path};
hdf5_internal::Hdf5Handle file{
hdf5_internal::RequireHdf5Id(
H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL,
H5P_DEFAULT, H5P_DEFAULT),
"Unable to create eigen candidate."),
H5Fclose};
hdf5_internal::WriterModelData model{
{}, state.constraint_mask, state.prescribed_displacement};
(void)hdf5_internal::CreateGroup(file.Get(), "/model");
(void)hdf5_internal::CreateGroup(file.Get(), "/model/shell");
hdf5_internal::WriteModel(file.Get(), domain, model);
WriteDensities(file.Get(), domain);
WriteDiagnostics(file.Get(), diagnostics);
WriteEigenMetadata(file.Get());
WriteStep(file.Get(), domain, state);
(void)hdf5_internal::CreateGroup(file.Get(), kEigenRoot);
WriteModes(file.Get(), state, domain.Frequency()->source_step_name_raw);
std::vector<double> vector_values;
vector_values.reserve(state.vectors.Rows() * state.vectors.Columns());
for (std::size_t row = 0U; row < state.vectors.Rows(); ++row)
for (std::size_t column = 0U; column < state.vectors.Columns(); ++column)
vector_values.push_back(state.vectors(row, column));
hdf5_internal::WriteModelDoubleDataset(
file.Get(), "/steps/Step-1/eigen/vectors",
{static_cast<hsize_t>(state.modes.size()),
static_cast<hsize_t>(domain.Nodes().size()), 6U},
vector_values.data(), vector_values.size(), "UX,UY,UZ,URX,URY,URZ",
"mass^-1/2,mass^-1/2,mass^-1/2,mass^-1/2*length^-1,mass^-1/"
"2*length^-1,mass^-1/2*length^-1",
"global-cartesian", "nodal-modal-vector");
{
hdf5_internal::Hdf5Handle vectors{
hdf5_internal::RequireHdf5Id(
H5Dopen2(file.Get(), "/steps/Step-1/eigen/vectors", H5P_DEFAULT),
"Unable to open modal vectors."),
H5Dclose};
WriteModalAttributes(vectors.Get(),
domain.Frequency()->source_step_name_raw);
hdf5_internal::WriteStringAttribute(vectors.Get(), "normalization",
"mass");
hdf5_internal::WriteStringAttribute(vectors.Get(), "mass_variant",
"derived-lumped");
hdf5_internal::WriteStringAttribute(
vectors.Get(), "interpretation",
"normalized-mode-not-response-displacement");
}
WriteVerification(file.Get(), state);
(void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1/eigen/backend");
WriteBackendAttributes(file.Get(), state);
WriteBaseFpm(file.Get(), state.backend.base_fpm);
WriteBackendCalls(file.Get(), state.backend.backend_calls);
hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL),
"Unable to flush eigen candidate.");
hdf5_internal::RequireHdf5(file.CloseChecked(),
"Unable to close eigen candidate.");
const Status checked = hdf5_internal::SelfCheckEigenFile(
candidate_path, domain, state, diagnostics);
const Status result = hdf5_internal::AtomicFinalizeValidatedCandidate(
candidate_path, output_path, checked);
if (result.IsOk()) cleanup.Release();
return result;
} catch (const hdf5_internal::Hdf5Failure& failure) {
return hdf5_internal::OutputFailure("hdf5-write-failure", failure.what());
} catch (const std::exception& failure) {
return hdf5_internal::OutputFailure("hdf5-write-failure", failure.what());
}
}
} // namespace fesa
Reference in New Issue View Git Blame Copy Permalink