feat(utils/serialize): System for serializing anything that can be serialized in SERiF (note LLM's made signifigant contributions to this module)
There are loads of things (solutions, models, solver descriptions, mesh's, etc..) which need to be exported from SERiF and then either reimported into SERiF or into some other program. This module is responsible for all of that. Generally with the goal of feeling like the fmt / new format module in C++23 and on (e.g. by specializing formatting structs) Note that LLM's made signifigant contributions to this module (ChatGPT 6 astra max and Claude Opus 5.5 Medium). All code has been reviewd manually
This commit is contained in:
424
src/lib/serif/utils/serialize/archive.cpp
Normal file
424
src/lib/serif/utils/serialize/archive.cpp
Normal file
@@ -0,0 +1,424 @@
|
||||
#include "serif/utils/serialize/archive.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <atomic>
|
||||
#include <charconv>
|
||||
#include <chrono>
|
||||
#include <cstdio>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <locale>
|
||||
#include <set>
|
||||
#include <sstream>
|
||||
#include <system_error>
|
||||
|
||||
#include "serif/config.hpp"
|
||||
#include "format.hpp"
|
||||
|
||||
namespace serif::utils::serialize {
|
||||
namespace {
|
||||
struct Header {
|
||||
Format format{Format::binary};
|
||||
std::string version;
|
||||
std::string timestamp;
|
||||
std::uint64_t index_offset{detail::header_bytes};
|
||||
std::uint64_t index_bytes{0};
|
||||
std::uint64_t data_offset{0};
|
||||
std::uint64_t file_bytes{0};
|
||||
std::uint64_t entries{0};
|
||||
std::string checksum;
|
||||
};
|
||||
|
||||
struct Payload {
|
||||
Entry entry;
|
||||
std::string bytes;
|
||||
};
|
||||
|
||||
/** @brief Keep the fixed header readable and independent of compiler structure padding. */
|
||||
std::string print_header(const Header& header) {
|
||||
std::ostringstream text;
|
||||
text.imbue(std::locale::classic());
|
||||
text << "SERIF_ARCHIVE 1\nformat " << (header.format == Format::binary ? "binary" : "text")
|
||||
<< "\nserif_version " << std::quoted(header.version) << "\nwritten_utc " << std::quoted(header.timestamp)
|
||||
<< "\nledger_offset " << header.index_offset << "\nledger_bytes " << header.index_bytes
|
||||
<< "\ndata_offset " << header.data_offset << "\nfile_bytes " << header.file_bytes
|
||||
<< "\nentries " << header.entries << "\nledger_checksum " << header.checksum << "\nEND_HEADER\n";
|
||||
std::string result = text.str();
|
||||
result.resize(detail::header_bytes, ' ');
|
||||
result.back() = '\n';
|
||||
return result;
|
||||
}
|
||||
|
||||
/** @brief Parse a fixed sequence so duplicated or unknown control fields cannot be ambiguous. */
|
||||
Result<Header> parse_header(std::string_view bytes) {
|
||||
std::istringstream input{std::string{bytes}};
|
||||
input.imbue(std::locale::classic());
|
||||
Header header;
|
||||
std::string name;
|
||||
std::string format;
|
||||
int version = 0;
|
||||
/** @brief Consume the next required control-field name. */
|
||||
const auto next = [&](std::string_view expected) { return static_cast<bool>(input >> name) && name == expected; };
|
||||
const bool valid = next("SERIF_ARCHIVE") && (input >> version) && version == 1 && next("format") && (input >> format)
|
||||
&& next("serif_version") && (input >> std::quoted(header.version)) && next("written_utc") && (input >> std::quoted(header.timestamp))
|
||||
&& next("ledger_offset") && (input >> header.index_offset) && next("ledger_bytes") && (input >> header.index_bytes)
|
||||
&& next("data_offset") && (input >> header.data_offset) && next("file_bytes") && (input >> header.file_bytes)
|
||||
&& next("entries") && (input >> header.entries) && next("ledger_checksum") && (input >> header.checksum) && next("END_HEADER");
|
||||
input >> std::ws;
|
||||
if (!valid || !input.eof() || (format != "binary" && format != "text")) {
|
||||
return std::unexpected(Error{"header", "Malformed header or unsupported archive format version."});
|
||||
}
|
||||
header.format = format == "binary" ? Format::binary : Format::text;
|
||||
return header;
|
||||
}
|
||||
|
||||
/** @brief Compare semantic major versions without interpreting arbitrary strings as zero. */
|
||||
std::optional<std::uint64_t> major_version(std::string_view version) {
|
||||
const std::size_t dot = version.find('.');
|
||||
if (dot == std::string_view::npos || dot == 0) {
|
||||
return {};
|
||||
}
|
||||
std::uint64_t major = 0;
|
||||
const auto parsed = std::from_chars(version.data(), version.data() + dot, major);
|
||||
if (parsed.ec != std::errc{} || parsed.ptr != version.data() + dot) {
|
||||
return {};
|
||||
}
|
||||
return major;
|
||||
}
|
||||
|
||||
/** @brief Flatten a component tree into independent payloads and reject ambiguous names. */
|
||||
Result<void> collect(const Node& node, std::string path, Format format, std::vector<Payload>& result, std::size_t depth = 0) {
|
||||
if (!detail::valid_path(path) || depth > 128 || node.version == 0 || node.type.empty() || node.type.size() > 1024
|
||||
|| node.type.find_first_of("\r\n\0", 0, 3) != std::string::npos) {
|
||||
return std::unexpected(Error{std::move(path), "Invalid component path, schema identifier, or nesting depth."});
|
||||
}
|
||||
if (!node.children.empty() && !std::holds_alternative<std::monostate>(node.value)) {
|
||||
return std::unexpected(Error{std::move(path), "A component must contain either fields or a payload."});
|
||||
}
|
||||
Payload payload{
|
||||
.entry = {.path = path, .type = node.type, .version = node.version, .kind = static_cast<std::uint32_t>(node.value.index()),
|
||||
.count = detail::value_count(node.value)},
|
||||
.bytes = detail::encode_value(node.value, format)};
|
||||
payload.entry.bytes = payload.bytes.size();
|
||||
payload.entry.checksum = fingerprint(payload.bytes);
|
||||
result.push_back(std::move(payload));
|
||||
for (const auto& [name, child] : node.children) {
|
||||
if (name.find('/') != std::string::npos) {
|
||||
return std::unexpected(Error{path, "A field name cannot contain a path separator."});
|
||||
}
|
||||
auto saved = collect(child, path + "/" + name, format, result, depth + 1);
|
||||
if (!saved) {
|
||||
return saved;
|
||||
}
|
||||
}
|
||||
return {};
|
||||
}
|
||||
|
||||
/** @brief Fixed-width offsets keep the ledger length unchanged when final locations are filled in. */
|
||||
std::string print_ledger(const std::vector<Payload>& payloads) {
|
||||
std::ostringstream output;
|
||||
output.imbue(std::locale::classic());
|
||||
output << "# path schema version kind count offset bytes checksum\n";
|
||||
for (const Payload& payload : payloads) {
|
||||
const Entry& entry = payload.entry;
|
||||
output << std::quoted(entry.path) << ' ' << std::quoted(entry.type) << ' ' << entry.version << ' ' << entry.kind << ' '
|
||||
<< entry.count << ' ' << std::setfill('0') << std::setw(20) << entry.offset << ' ' << std::setw(20) << entry.bytes
|
||||
<< ' ' << entry.checksum << '\n';
|
||||
}
|
||||
return output.str();
|
||||
}
|
||||
|
||||
/** @brief Give each text payload a visible section label outside its indexed bytes. */
|
||||
std::string section_label(const Entry& entry, Format format) {
|
||||
return format == Format::text ? "\n# " + entry.path + "\n" : "";
|
||||
}
|
||||
|
||||
/** @brief Remove a failed temporary write without ever removing the previous destination. */
|
||||
struct TemporaryFile {
|
||||
std::filesystem::path path;
|
||||
std::FILE* stream{nullptr};
|
||||
bool owned{false};
|
||||
/** @brief Clean up regardless of which validation or I/O operation failed. */
|
||||
~TemporaryFile() {
|
||||
if (stream != nullptr) {
|
||||
std::fclose(stream);
|
||||
}
|
||||
if (owned) {
|
||||
std::error_code error;
|
||||
std::filesystem::remove(path, error);
|
||||
}
|
||||
}
|
||||
/** @brief Check every byte of an output section. */
|
||||
bool append(std::string_view bytes) const { return std::fwrite(bytes.data(), 1, bytes.size(), stream) == bytes.size(); }
|
||||
};
|
||||
|
||||
/** @brief Read an exact bounded region, clearing EOF from any preceding failed read. */
|
||||
Result<std::string> read_bytes(std::ifstream& file, std::uint64_t offset, std::uint64_t bytes, std::string_view path) {
|
||||
if (bytes > static_cast<std::uint64_t>(std::numeric_limits<std::streamsize>::max())
|
||||
|| offset > static_cast<std::uint64_t>(std::numeric_limits<std::streamoff>::max())) {
|
||||
return std::unexpected(Error{std::string{path}, "File region exceeds this reader's addressable range."});
|
||||
}
|
||||
file.clear();
|
||||
file.seekg(static_cast<std::streamoff>(offset));
|
||||
std::string result(bytes, '\0');
|
||||
file.read(result.data(), static_cast<std::streamsize>(bytes));
|
||||
if (!file) {
|
||||
return std::unexpected(Error{std::string{path}, "Unable to read the complete indexed file region."});
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
/** @brief Validate all structural bounds before any payload is allocated or interpreted. */
|
||||
Result<std::vector<Entry>> parse_ledger(std::string_view bytes, const Header& header, const ReadOptions& options) {
|
||||
std::istringstream input{std::string{bytes}};
|
||||
input.imbue(std::locale::classic());
|
||||
std::string line;
|
||||
std::getline(input, line);
|
||||
if (line != "# path schema version kind count offset bytes checksum") {
|
||||
return std::unexpected(Error{"ledger", "Unrecognized ledger columns."});
|
||||
}
|
||||
std::vector<Entry> entries;
|
||||
std::map<std::string, std::uint32_t> kinds;
|
||||
std::uint64_t previous_end = header.data_offset;
|
||||
while (std::getline(input, line)) {
|
||||
Entry entry;
|
||||
std::istringstream row{line};
|
||||
row.imbue(std::locale::classic());
|
||||
const bool parsed = static_cast<bool>(row >> std::quoted(entry.path) >> std::quoted(entry.type) >> entry.version >> entry.kind
|
||||
>> entry.count >> entry.offset >> entry.bytes >> entry.checksum);
|
||||
row >> std::ws;
|
||||
if (!parsed || !row.eof() || !detail::valid_path(entry.path) || entry.type.empty() || entry.type.size() > 1024
|
||||
|| entry.kind > 8 || entry.version == 0 || entries.size() >= options.maximum_entries
|
||||
|| static_cast<std::uint64_t>(std::count(entry.path.begin(), entry.path.end(), '/')) > options.maximum_depth
|
||||
|| (!entries.empty() && entries.back().path >= entry.path)) {
|
||||
return std::unexpected(Error{entry.path, "Malformed, duplicated, unsorted, or oversized ledger entry."});
|
||||
}
|
||||
const std::size_t separator = entry.path.rfind('/');
|
||||
if (separator != std::string::npos) {
|
||||
const auto parent = kinds.find(entry.path.substr(0, separator));
|
||||
if (parent == kinds.end() || parent->second != 0) {
|
||||
return std::unexpected(Error{entry.path, "A field's parent must be an existing component directory."});
|
||||
}
|
||||
}
|
||||
if (entry.kind == 0) {
|
||||
if (entry.count != 0 || entry.offset != 0 || entry.bytes != 0 || entry.checksum != fingerprint("")) {
|
||||
return std::unexpected(Error{entry.path, "Directory entries cannot reference payload bytes."});
|
||||
}
|
||||
} else {
|
||||
if (entry.offset < previous_end || entry.offset > header.file_bytes || entry.bytes > header.file_bytes - entry.offset
|
||||
|| entry.bytes > options.maximum_payload_bytes || entry.count > options.maximum_payload_bytes) {
|
||||
return std::unexpected(Error{entry.path, "Payload overlaps another section, exceeds the file, or exceeds configured limits."});
|
||||
}
|
||||
previous_end = entry.offset + entry.bytes;
|
||||
}
|
||||
kinds.emplace(entry.path, entry.kind);
|
||||
entries.push_back(std::move(entry));
|
||||
}
|
||||
if (entries.size() != header.entries || !kinds.contains("provenance")) {
|
||||
return std::unexpected(Error{"ledger", "Entry count is inconsistent or required provenance is missing."});
|
||||
}
|
||||
return entries;
|
||||
}
|
||||
}
|
||||
|
||||
void add_component(Pack& contents, Field component) {
|
||||
if (contents.error) {
|
||||
return;
|
||||
}
|
||||
if (component.name == "provenance" || !detail::valid_path(component.name) || component.name.find('/') != std::string::npos) {
|
||||
contents.error = Error{component.name, "Use a non-reserved, single-component archive name."};
|
||||
return;
|
||||
}
|
||||
const auto found = contents.root.children.find(component.name);
|
||||
if (found != contents.root.children.end() && found->second != component.node) {
|
||||
contents.error = Error{component.name, "Two different components were packed under the same name."};
|
||||
return;
|
||||
}
|
||||
contents.root.children.emplace(std::move(component.name), std::move(component.node));
|
||||
}
|
||||
|
||||
Result<void> write(const std::filesystem::path& path, const Pack& contents, const WriteOptions& options) {
|
||||
if (contents.error) {
|
||||
return std::unexpected(*contents.error);
|
||||
}
|
||||
std::vector<Payload> payloads;
|
||||
const Node metadata = detail::provenance(options);
|
||||
auto saved = collect(metadata, "provenance", options.format, payloads);
|
||||
if (!saved) {
|
||||
return saved;
|
||||
}
|
||||
for (const auto& [name, node] : contents.root.children) {
|
||||
if (name == "provenance") {
|
||||
return std::unexpected(Error{name, "Provenance is written by the archive, not supplied as a model component."});
|
||||
}
|
||||
saved = collect(node, name, options.format, payloads);
|
||||
if (!saved) {
|
||||
return saved;
|
||||
}
|
||||
}
|
||||
std::ranges::sort(payloads, {}, [](const Payload& payload) { return payload.entry.path; });
|
||||
Header header{.format = options.format, .version = SERIF_VERSION,
|
||||
.timestamp = std::get<std::string>(metadata.children.at("written_utc").value),
|
||||
.index_bytes = print_ledger(payloads).size(), .entries = payloads.size()};
|
||||
header.data_offset = header.index_offset + header.index_bytes;
|
||||
std::uint64_t position = header.data_offset;
|
||||
for (Payload& payload : payloads) {
|
||||
if (payload.entry.kind != 0) {
|
||||
position += section_label(payload.entry, options.format).size();
|
||||
payload.entry.offset = position;
|
||||
position += payload.bytes.size();
|
||||
}
|
||||
}
|
||||
header.file_bytes = position;
|
||||
const std::string ledger = print_ledger(payloads);
|
||||
header.checksum = fingerprint(ledger);
|
||||
|
||||
static std::atomic_uint64_t sequence{0};
|
||||
TemporaryFile temporary;
|
||||
for (int attempt = 0; attempt < 10 && temporary.stream == nullptr; ++attempt) {
|
||||
temporary.path = path;
|
||||
temporary.path += ".tmp-" + std::to_string(std::chrono::steady_clock::now().time_since_epoch().count()) + "-" + std::to_string(sequence++);
|
||||
temporary.stream = std::fopen(temporary.path.string().c_str(), "wbx");
|
||||
}
|
||||
if (temporary.stream == nullptr) {
|
||||
return std::unexpected(Error{path.string(), "Cannot create an exclusive temporary file beside the destination."});
|
||||
}
|
||||
temporary.owned = true;
|
||||
bool written = temporary.append(print_header(header)) && temporary.append(ledger);
|
||||
for (const Payload& payload : payloads) {
|
||||
if (written && payload.entry.kind != 0) {
|
||||
written = temporary.append(section_label(payload.entry, options.format)) && temporary.append(payload.bytes);
|
||||
}
|
||||
}
|
||||
const int close_status = std::fclose(temporary.stream);
|
||||
temporary.stream = nullptr;
|
||||
if (!written || close_status != 0) {
|
||||
return std::unexpected(Error{path.string(), "The archive could not be written completely; the previous destination was retained."});
|
||||
}
|
||||
std::error_code error;
|
||||
std::filesystem::rename(temporary.path, path, error);
|
||||
if (error) {
|
||||
return std::unexpected(Error{path.string(), "Cannot publish the completed archive: " + error.message()});
|
||||
}
|
||||
temporary.owned = false;
|
||||
return {};
|
||||
}
|
||||
|
||||
Result<Reader> Reader::open(const std::filesystem::path& path, ReadOptions options) {
|
||||
Reader reader;
|
||||
reader.m_options = options;
|
||||
reader.m_file.open(path, std::ios::binary | std::ios::ate);
|
||||
if (!reader.m_file) {
|
||||
return std::unexpected(Error{path.string(), "Cannot open archive for reading."});
|
||||
}
|
||||
const std::streamoff file_size = reader.m_file.tellg();
|
||||
if (file_size < static_cast<std::streamoff>(detail::header_bytes)) {
|
||||
return std::unexpected(Error{"header", "Archive is truncated before the complete header."});
|
||||
}
|
||||
auto bytes = read_bytes(reader.m_file, 0, detail::header_bytes, "header");
|
||||
if (!bytes) {
|
||||
return std::unexpected(bytes.error());
|
||||
}
|
||||
auto header = parse_header(*bytes);
|
||||
if (!header) {
|
||||
return std::unexpected(header.error());
|
||||
}
|
||||
const bool same_major = major_version(header->version) && major_version(header->version) == major_version(SERIF_VERSION);
|
||||
if ((options.version_tolerance == VersionTolerance::exact && header->version != SERIF_VERSION)
|
||||
|| (options.version_tolerance == VersionTolerance::same_major && !same_major)) {
|
||||
return std::unexpected(Error{"header/serif_version", "SERiF " + header->version + " is incompatible with the selected reader tolerance (reader " SERIF_VERSION ")."});
|
||||
}
|
||||
|
||||
if (header->index_offset != detail::header_bytes || header->file_bytes != static_cast<std::uint64_t>(file_size)
|
||||
|| header->index_bytes > options.maximum_index_bytes || header->index_bytes > header->file_bytes - detail::header_bytes
|
||||
|| header->data_offset != header->index_offset + header->index_bytes || header->entries > options.maximum_entries) {
|
||||
|
||||
std::string cmp_string = std::format("header: [index_offset={}, index_bytes={}, data_offset={}, file_bytes={}, entries={}]", header->index_offset, header->index_bytes, header->data_offset, header->file_bytes, header->entries);
|
||||
cmp_string = std::format("{}, file: [size={}]", cmp_string, file_size);
|
||||
cmp_string = std::format("{}, limits: [max_index_bytes={}, max_entries={}]", cmp_string, options.maximum_index_bytes, options.maximum_entries);
|
||||
cmp_string = std::format("{}, version_tolerance={}", cmp_string, static_cast<int>(options.version_tolerance));
|
||||
cmp_string = std::format("{}, detail::header_bytes={}", cmp_string, detail::header_bytes);
|
||||
|
||||
return std::unexpected(Error{"header", std::format("Invalid file size, section offsets, or configured allocation limits. {}", cmp_string)});
|
||||
}
|
||||
auto ledger = read_bytes(reader.m_file, header->index_offset, header->index_bytes, "ledger");
|
||||
if (!ledger) {
|
||||
return std::unexpected(ledger.error());
|
||||
}
|
||||
if (fingerprint(*ledger) != header->checksum) {
|
||||
return std::unexpected(Error{"ledger", "Ledger checksum does not match the header."});
|
||||
}
|
||||
auto entries = parse_ledger(*ledger, *header, options);
|
||||
if (!entries) {
|
||||
return std::unexpected(entries.error());
|
||||
}
|
||||
reader.m_entries = std::move(*entries);
|
||||
reader.m_format = header->format;
|
||||
reader.m_version = std::move(header->version);
|
||||
return reader;
|
||||
}
|
||||
|
||||
bool Reader::contains(std::string_view path) const {
|
||||
const auto found = std::ranges::lower_bound(m_entries, path, {}, &Entry::path);
|
||||
return found != m_entries.end() && found->path == path;
|
||||
}
|
||||
|
||||
Result<Node> Reader::node(std::string_view path) const {
|
||||
const auto first = std::ranges::lower_bound(m_entries, path, {}, &Entry::path);
|
||||
if (first == m_entries.end() || first->path != path) {
|
||||
return std::unexpected(Error{std::string{path}, "Requested component is not present in the archive ledger."});
|
||||
}
|
||||
Node root;
|
||||
std::map<std::string, Node*> nodes;
|
||||
const std::string prefix = std::string{path} + "/";
|
||||
std::uint64_t total_bytes = 0;
|
||||
for (auto entry = first; entry != m_entries.end(); ++entry) {
|
||||
if (entry != first && !entry->path.starts_with(prefix)) {
|
||||
// Sibling names can sort before a parent's descendants (a.b precedes a/b).
|
||||
if (entry->path < prefix) {
|
||||
continue;
|
||||
}
|
||||
break;
|
||||
}
|
||||
if (entry->bytes > m_options.maximum_payload_bytes - total_bytes) {
|
||||
return std::unexpected(Error{entry->path, "Selected subtree exceeds the configured payload budget."});
|
||||
}
|
||||
total_bytes += entry->bytes;
|
||||
auto bytes = read_bytes(m_file, entry->offset, entry->bytes, entry->path);
|
||||
if (!bytes) {
|
||||
return std::unexpected(bytes.error());
|
||||
}
|
||||
if (fingerprint(*bytes) != entry->checksum) {
|
||||
return std::unexpected(Error{entry->path, "Payload checksum mismatch; the component is damaged."});
|
||||
}
|
||||
auto value = detail::decode_value(*bytes, *entry, m_format);
|
||||
if (!value) {
|
||||
return std::unexpected(value.error());
|
||||
}
|
||||
Node* destination = &root;
|
||||
if (entry != first) {
|
||||
const std::size_t separator = entry->path.rfind('/');
|
||||
const auto parent = nodes.find(entry->path.substr(0, separator));
|
||||
if (parent == nodes.end()) {
|
||||
return std::unexpected(Error{entry->path, "Requested subtree has a missing parent."});
|
||||
}
|
||||
destination = &parent->second->children[entry->path.substr(separator + 1)];
|
||||
}
|
||||
*destination = Node{.type = entry->type, .version = entry->version, .value = std::move(*value)};
|
||||
nodes.emplace(entry->path, destination);
|
||||
}
|
||||
return root;
|
||||
}
|
||||
|
||||
Result<void> Reader::validate() const {
|
||||
for (const Entry& entry : m_entries) {
|
||||
if (entry.kind != 0) {
|
||||
auto payload = node(entry.path);
|
||||
if (!payload) {
|
||||
return std::unexpected(payload.error());
|
||||
}
|
||||
}
|
||||
}
|
||||
return {};
|
||||
}
|
||||
}
|
||||
230
src/lib/serif/utils/serialize/format.cpp
Normal file
230
src/lib/serif/utils/serialize/format.cpp
Normal file
@@ -0,0 +1,230 @@
|
||||
#include "format.hpp"
|
||||
|
||||
#include <bit>
|
||||
#include <charconv>
|
||||
#include <cmath>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <locale>
|
||||
#include <sstream>
|
||||
#include <type_traits>
|
||||
|
||||
namespace serif::utils::serialize {
|
||||
std::string fingerprint(std::string_view bytes) {
|
||||
std::uint64_t hash = 14695981039346656037ULL;
|
||||
for (unsigned char byte : bytes) {
|
||||
hash ^= byte;
|
||||
hash *= 1099511628211ULL;
|
||||
}
|
||||
std::ostringstream output;
|
||||
output.imbue(std::locale::classic());
|
||||
output << "fnv1a64-v1:" << std::hex << std::setfill('0') << std::setw(16) << hash;
|
||||
return output.str();
|
||||
}
|
||||
}
|
||||
|
||||
namespace serif::utils::serialize::detail {
|
||||
namespace {
|
||||
/** @brief Emit one integer in little-endian order, independent of the writing host. */
|
||||
void append_word(std::string& output, std::uint64_t word) {
|
||||
for (int shift = 0; shift < 64; shift += 8) {
|
||||
output.push_back(static_cast<char>((word >> shift) & 0xff));
|
||||
}
|
||||
}
|
||||
|
||||
/** @brief Recover one unaligned little-endian word without aliasing native objects. */
|
||||
std::uint64_t read_word(std::string_view bytes, std::size_t index) {
|
||||
std::uint64_t word = 0;
|
||||
for (int byte = 0; byte < 8; ++byte) {
|
||||
word |= static_cast<std::uint64_t>(static_cast<unsigned char>(bytes[index * 8 + byte])) << (8 * byte);
|
||||
}
|
||||
return word;
|
||||
}
|
||||
|
||||
/** @brief Keep finite text readable and preserve even the payload bits of nonfinite doubles. */
|
||||
void write_real(std::ostream& output, double value) {
|
||||
if (std::isfinite(value)) {
|
||||
output << std::setprecision(std::numeric_limits<double>::max_digits10) << value;
|
||||
} else {
|
||||
output << "bits:" << std::hex << std::setfill('0') << std::setw(16) << std::bit_cast<std::uint64_t>(value) << std::dec;
|
||||
}
|
||||
}
|
||||
|
||||
/** @brief Parse text doubles with an explicit representation for nonfinite values. */
|
||||
bool read_real(std::istream& input, double& value) {
|
||||
std::string token;
|
||||
if (!(input >> token)) {
|
||||
return false;
|
||||
}
|
||||
if (token.starts_with("bits:")) {
|
||||
std::uint64_t bits = 0;
|
||||
const auto parsed = std::from_chars(token.data() + 5, token.data() + token.size(), bits, 16);
|
||||
if (token.size() != 21 || parsed.ec != std::errc{} || parsed.ptr != token.data() + token.size()) {
|
||||
return false;
|
||||
}
|
||||
value = std::bit_cast<double>(bits);
|
||||
return true;
|
||||
}
|
||||
const auto parsed = std::from_chars(token.data(), token.data() + token.size(), value);
|
||||
return parsed.ec == std::errc{} && parsed.ptr == token.data() + token.size();
|
||||
}
|
||||
|
||||
/** @brief Require an exact array extent, including no unconsumed trailing text. */
|
||||
template <typename T>
|
||||
Result<Value> read_array(std::string_view bytes, const Entry& entry, Format format) {
|
||||
if ((format == Format::binary && (bytes.size() % 8 != 0 || entry.count != bytes.size() / 8)) || entry.count > bytes.size()) {
|
||||
return std::unexpected(Error{entry.path, "Array extent does not match the payload length."});
|
||||
}
|
||||
std::vector<T> values(entry.count);
|
||||
if (format == Format::binary) {
|
||||
for (std::size_t index = 0; index < values.size(); ++index) {
|
||||
values[index] = std::bit_cast<T>(read_word(bytes, index));
|
||||
}
|
||||
} else {
|
||||
std::istringstream input{std::string{bytes}};
|
||||
input.imbue(std::locale::classic());
|
||||
for (T& value : values) {
|
||||
bool valid;
|
||||
if constexpr (std::same_as<T, double>) {
|
||||
valid = read_real(input, value);
|
||||
} else {
|
||||
std::string token;
|
||||
input >> token;
|
||||
const auto parsed = std::from_chars(token.data(), token.data() + token.size(), value);
|
||||
valid = !token.empty() && parsed.ec == std::errc{} && parsed.ptr == token.data() + token.size();
|
||||
}
|
||||
if (!valid) {
|
||||
return std::unexpected(Error{entry.path, "Invalid number in text payload."});
|
||||
}
|
||||
}
|
||||
input >> std::ws;
|
||||
if (!input.eof()) {
|
||||
return std::unexpected(Error{entry.path, "Unexpected trailing text in numeric payload."});
|
||||
}
|
||||
}
|
||||
return Value{std::move(values)};
|
||||
}
|
||||
|
||||
/** @brief Reuse checked array decoding for fixed-width scalars. */
|
||||
template <typename T>
|
||||
Result<Value> read_scalar(std::string_view bytes, const Entry& entry, Format format) {
|
||||
if (entry.count != 1) {
|
||||
return std::unexpected(Error{entry.path, "A scalar must have exactly one value."});
|
||||
}
|
||||
auto values = read_array<T>(bytes, entry, format);
|
||||
if (!values) {
|
||||
return std::unexpected(values.error());
|
||||
}
|
||||
return Value{std::get<std::vector<T>>(*values).front()};
|
||||
}
|
||||
}
|
||||
|
||||
std::uint64_t value_count(const Value& value) {
|
||||
return std::visit([](const auto& stored) -> std::uint64_t {
|
||||
using T = std::remove_cvref_t<decltype(stored)>;
|
||||
if constexpr (std::same_as<T, std::monostate>) {
|
||||
return 0;
|
||||
} else if constexpr (requires { stored.size(); }) {
|
||||
return stored.size();
|
||||
} else {
|
||||
return 1;
|
||||
}
|
||||
}, value);
|
||||
}
|
||||
|
||||
std::string encode_value(const Value& value, Format format) {
|
||||
std::string binary;
|
||||
std::ostringstream text;
|
||||
text.imbue(std::locale::classic());
|
||||
/** @brief Emit one numerical value through the selected representation. */
|
||||
const auto append_number = [&](auto number) {
|
||||
using T = decltype(number);
|
||||
if (format == Format::binary) {
|
||||
append_word(binary, std::bit_cast<std::uint64_t>(number));
|
||||
} else {
|
||||
if constexpr (std::same_as<T, double>) {
|
||||
write_real(text, number);
|
||||
} else {
|
||||
text << number;
|
||||
}
|
||||
text << '\n';
|
||||
}
|
||||
};
|
||||
std::visit([&](const auto& stored) {
|
||||
using T = std::remove_cvref_t<decltype(stored)>;
|
||||
if constexpr (std::same_as<T, bool>) {
|
||||
binary.push_back(stored ? '\1' : '\0');
|
||||
text << (stored ? '1' : '0');
|
||||
} else if constexpr (std::same_as<T, std::string>) {
|
||||
binary = stored;
|
||||
text << std::quoted(stored);
|
||||
} else if constexpr (std::is_arithmetic_v<T>) {
|
||||
append_number(stored);
|
||||
} else if constexpr (!std::same_as<T, std::monostate>) {
|
||||
for (auto number : stored) {
|
||||
append_number(number);
|
||||
}
|
||||
}
|
||||
}, value);
|
||||
return format == Format::binary ? binary : text.str();
|
||||
}
|
||||
|
||||
Result<Value> decode_value(std::string_view bytes, const Entry& entry, Format format) {
|
||||
switch (entry.kind) {
|
||||
case 0: {
|
||||
if (entry.count != 0 || !bytes.empty()) {
|
||||
return std::unexpected(Error{entry.path, "A component directory cannot contain a scalar payload."});
|
||||
}
|
||||
return Value{};
|
||||
}
|
||||
case 1: {
|
||||
const char zero = format == Format::binary ? '\0' : '0';
|
||||
const char one = format == Format::binary ? '\1' : '1';
|
||||
if (entry.count != 1 || bytes.size() != 1 || (bytes[0] != zero && bytes[0] != one)) {
|
||||
return std::unexpected(Error{entry.path, "Invalid boolean payload."});
|
||||
}
|
||||
return Value{bytes[0] == one};
|
||||
}
|
||||
case 2: { return read_scalar<std::int64_t>(bytes, entry, format); }
|
||||
case 3: { return read_scalar<std::uint64_t>(bytes, entry, format); }
|
||||
case 4: { return read_scalar<double>(bytes, entry, format); }
|
||||
case 5: {
|
||||
std::string value;
|
||||
if (format == Format::binary) {
|
||||
value = bytes;
|
||||
} else {
|
||||
std::istringstream input{std::string{bytes}};
|
||||
input.imbue(std::locale::classic());
|
||||
if (input.peek() != '"' || !(input >> std::quoted(value))) {
|
||||
return std::unexpected(Error{entry.path, "Invalid quoted text payload."});
|
||||
}
|
||||
input >> std::ws;
|
||||
if (!input.eof()) {
|
||||
return std::unexpected(Error{entry.path, "Unexpected trailing text after string."});
|
||||
}
|
||||
}
|
||||
if (value.size() != entry.count) {
|
||||
return std::unexpected(Error{entry.path, "String length does not match the ledger."});
|
||||
}
|
||||
return Value{std::move(value)};
|
||||
}
|
||||
case 6: { return read_array<double>(bytes, entry, format); }
|
||||
case 7: { return read_array<std::int64_t>(bytes, entry, format); }
|
||||
case 8: { return read_array<std::uint64_t>(bytes, entry, format); }
|
||||
default: { return std::unexpected(Error{entry.path, "Unknown payload representation."}); }
|
||||
}
|
||||
}
|
||||
|
||||
bool valid_path(std::string_view path) {
|
||||
if (path.empty() || path.size() > 4096 || path.front() == '/' || path.back() == '/' || path.find("//") != std::string_view::npos) {
|
||||
return false;
|
||||
}
|
||||
for (unsigned char character : path) {
|
||||
const bool alphabet = (character >= 'a' && character <= 'z') || (character >= 'A' && character <= 'Z');
|
||||
if (!alphabet && !(character >= '0' && character <= '9') && character != '_' && character != '-' && character != '.' && character != '/') {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
}
|
||||
21
src/lib/serif/utils/serialize/format.hpp
Normal file
21
src/lib/serif/utils/serialize/format.hpp
Normal file
@@ -0,0 +1,21 @@
|
||||
#pragma once
|
||||
|
||||
#include <string>
|
||||
#include <string_view>
|
||||
|
||||
#include "serif/utils/serialize/archive.hpp"
|
||||
|
||||
namespace serif::utils::serialize::detail {
|
||||
inline constexpr std::uint64_t header_bytes = 512;
|
||||
|
||||
/** @brief Scalar types have one value; numerical arrays and strings carry their explicit extent. */
|
||||
[[nodiscard]] std::uint64_t value_count(const Value& value);
|
||||
/** @brief Encode one indexed payload without relying on native object layout or byte order. */
|
||||
[[nodiscard]] std::string encode_value(const Value& value, Format format);
|
||||
/** @brief Validate the extent and representation of an independently read payload. */
|
||||
[[nodiscard]] Result<Value> decode_value(std::string_view bytes, const Entry& entry, Format format);
|
||||
/** @brief Capture provenance once at the beginning of a write. */
|
||||
[[nodiscard]] Node provenance(const WriteOptions& options);
|
||||
/** @brief Restrict field names so path separators and line breaks cannot corrupt the ledger. */
|
||||
[[nodiscard]] bool valid_path(std::string_view path);
|
||||
}
|
||||
86
src/lib/serif/utils/serialize/provenance.cpp
Normal file
86
src/lib/serif/utils/serialize/provenance.cpp
Normal file
@@ -0,0 +1,86 @@
|
||||
#include "format.hpp"
|
||||
|
||||
#include <bit>
|
||||
#include <chrono>
|
||||
#include <ctime>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <locale>
|
||||
#include <sstream>
|
||||
#include <thread>
|
||||
|
||||
#include <nameof.hpp>
|
||||
#include <magic_enum/magic_enum.hpp>
|
||||
#include "serif/config.hpp"
|
||||
|
||||
#if defined(__unix__) || defined(__APPLE__)
|
||||
#include <sys/utsname.h>
|
||||
#endif
|
||||
|
||||
namespace serif::utils::serialize::detail {
|
||||
namespace {
|
||||
/** @brief Record the writing instant in an unambiguous UTC representation. */
|
||||
std::string timestamp() {
|
||||
const std::time_t now = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now());
|
||||
std::tm utc{};
|
||||
#if defined(_WIN32)
|
||||
gmtime_s(&utc, &now);
|
||||
#else
|
||||
gmtime_r(&now, &utc);
|
||||
#endif
|
||||
std::ostringstream output;
|
||||
output.imbue(std::locale::classic());
|
||||
output << std::put_time(&utc, "%Y-%m-%dT%H:%M:%SZ");
|
||||
return output.str();
|
||||
}
|
||||
|
||||
/** @brief Keep semantic library versions readable without relying on implementation-specific formatters. */
|
||||
std::string version(int major, int minor, int patch) {
|
||||
return std::to_string(major) + "." + std::to_string(minor) + "." + std::to_string(patch);
|
||||
}
|
||||
|
||||
/** @brief Describe the host without storing a user's account, environment variables, or home directory. */
|
||||
Node host_information() {
|
||||
Node host = object("host",
|
||||
field("operating_system", std::string{SERIF_BUILD_SYSTEM}),
|
||||
field("architecture", std::string{SERIF_BUILD_ARCHITECTURE}),
|
||||
field("hardware_threads", std::thread::hardware_concurrency()),
|
||||
field("pointer_bits", static_cast<unsigned int>(8 * sizeof(void*))),
|
||||
field("byte_order", std::string{std::endian::native == std::endian::little ? "little" : "big"}));
|
||||
#if defined(__unix__) || defined(__APPLE__)
|
||||
struct utsname information{};
|
||||
if (uname(&information) == 0) {
|
||||
host.children.emplace("name", encode(std::string{information.nodename}));
|
||||
host.children.emplace("kernel", encode(std::string{information.sysname}));
|
||||
host.children.emplace("kernel_release", encode(std::string{information.release}));
|
||||
host.children.emplace("kernel_version", encode(std::string{information.version}));
|
||||
host.children.emplace("machine", encode(std::string{information.machine}));
|
||||
}
|
||||
#endif
|
||||
return host;
|
||||
}
|
||||
}
|
||||
|
||||
Node provenance(const WriteOptions& options) {
|
||||
static_assert(sizeof(double) == 8 && std::numeric_limits<double>::is_iec559,
|
||||
"The native archive currently represents IEEE-754 binary64 coefficients.");
|
||||
Node libraries = object("libraries",
|
||||
field("nameof", version(NAMEOF_VERSION_MAJOR, NAMEOF_VERSION_MINOR, NAMEOF_VERSION_PATCH)),
|
||||
field("magic_enum", version(MAGIC_ENUM_VERSION_MAJOR, MAGIC_ENUM_VERSION_MINOR, MAGIC_ENUM_VERSION_PATCH)),
|
||||
field("mfem_build_dependency", std::string{SERIF_MFEM_BUILD_VERSION}));
|
||||
#if defined(_LIBCPP_VERSION)
|
||||
libraries.children.emplace("libcxx", encode(std::to_string(_LIBCPP_VERSION)));
|
||||
#elif defined(__GLIBCXX__)
|
||||
libraries.children.emplace("libstdcxx", encode(std::to_string(__GLIBCXX__)));
|
||||
#endif
|
||||
Node annotations{.type = "metadata"};
|
||||
for (const auto& [key, value] : options.metadata) {
|
||||
annotations.children.emplace(key, encode(value));
|
||||
}
|
||||
return object("provenance",
|
||||
field("written_utc", timestamp()), field("serif_version", std::string{SERIF_VERSION}),
|
||||
field("host", host_information()), field("libraries", libraries),
|
||||
field("compiler", std::string{SERIF_BUILD_COMPILER}), field("build_type", std::string{SERIF_BUILD_TYPE}),
|
||||
field("author", options.author), field("comment", options.comment), field("metadata", annotations));
|
||||
}
|
||||
}
|
||||
425
src/lib/serif/utils/serialize/top.cpp
Normal file
425
src/lib/serif/utils/serialize/top.cpp
Normal file
@@ -0,0 +1,425 @@
|
||||
#include "serif/utils/serialize/top.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <atomic>
|
||||
#include <cerrno>
|
||||
#include <chrono>
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstring>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <locale>
|
||||
#include <numbers>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
#if defined(__linux__)
|
||||
#include <fcntl.h>
|
||||
#include <sys/syscall.h>
|
||||
#include <unistd.h>
|
||||
#endif
|
||||
|
||||
namespace serif::utils::serialize {
|
||||
namespace {
|
||||
constexpr double pi = std::numbers::pi;
|
||||
using Matrix = std::vector<std::vector<double>>;
|
||||
|
||||
struct Failure { Error error; };
|
||||
[[noreturn]] void fail(std::string path, std::string message) {
|
||||
throw Failure{Error{std::move(path), std::move(message)}};
|
||||
}
|
||||
double checked(Result<double> value, const std::string& path) {
|
||||
if (!value) {
|
||||
fail(path + (value.error().path.empty() ? "" : "/" + value.error().path), value.error().message);
|
||||
}
|
||||
if (!std::isfinite(*value)) { fail(path, "The source returned a nonfinite value."); }
|
||||
return *value;
|
||||
}
|
||||
bool positive(double value) { return std::isfinite(value) && value > 0.; }
|
||||
|
||||
/** @brief Ordinary Legendre polynomials; normalization is applied exactly once at projection. */
|
||||
std::vector<double> legendre(double mu, int degree) {
|
||||
std::vector<double> result(degree + 1, 1.);
|
||||
if (degree) { result[1] = mu; }
|
||||
for (int l = 2; l <= degree; ++l) {
|
||||
result[l] = ((2*l - 1)*mu*result[l-1] - (l-1)*result[l-2])/l;
|
||||
}
|
||||
return result;
|
||||
}
|
||||
std::pair<std::vector<double>, std::vector<double>> gauss(int count) {
|
||||
std::vector<double> nodes(count), weights(count);
|
||||
for (int i = 0; i < (count+1)/2; ++i) {
|
||||
double x = std::cos(pi*(i+.75)/(count+.5));
|
||||
double derivative = 0.;
|
||||
for (int step = 0; step < 100; ++step) {
|
||||
const auto p = legendre(x, count);
|
||||
derivative = count*(x*p[count]-p[count-1])/(x*x-1.);
|
||||
const double next = x-p[count]/derivative;
|
||||
if (std::abs(next-x) <= 4*std::numeric_limits<double>::epsilon()) {
|
||||
x = next;
|
||||
const auto final = legendre(x, count);
|
||||
derivative = count*(x*final[count]-final[count-1])/(x*x-1.);
|
||||
break;
|
||||
}
|
||||
x = next;
|
||||
if (step == 99) { fail("options/angular_points", "Gauss-Legendre nodes did not converge."); }
|
||||
}
|
||||
nodes[i] = -x; nodes[count-1-i] = x;
|
||||
weights[i] = weights[count-1-i] = 2./((1-x*x)*derivative*derivative);
|
||||
}
|
||||
return {std::move(nodes), std::move(weights)};
|
||||
}
|
||||
double reconstruct(const std::vector<double>& coefficients, double mu) {
|
||||
const auto p = legendre(mu, static_cast<int>(coefficients.size())-1);
|
||||
long double sum = 0.;
|
||||
for (std::size_t l = 0; l < coefficients.size(); l += 2) {
|
||||
sum += coefficients[l]*std::sqrt((2.*l+1)/(4*pi))*p[l];
|
||||
}
|
||||
return static_cast<double>(sum);
|
||||
}
|
||||
std::vector<double> project(const std::vector<double>& values, const std::vector<double>& mu,
|
||||
const std::vector<double>& weights, int degree) {
|
||||
std::vector<long double> sums(degree+1, 0.);
|
||||
for (std::size_t j = 0; j < mu.size(); ++j) {
|
||||
const auto p = legendre(mu[j], degree);
|
||||
for (int l = 0; l <= degree; l += 2) { sums[l] += values[j]*weights[j]*p[l]; }
|
||||
}
|
||||
std::vector<double> result(degree+1, 0.);
|
||||
for (int l = 0; l <= degree; l += 2) {
|
||||
result[l] = static_cast<double>(sums[l])*std::sqrt(pi*(2*l+1));
|
||||
}
|
||||
return result;
|
||||
}
|
||||
Matrix derivative_matrix(const std::vector<double>& zeta) {
|
||||
const int count = static_cast<int>(zeta.size());
|
||||
Matrix result(count, std::vector<double>(count));
|
||||
for (int i = 0; i < count; ++i) {
|
||||
const double wi = (i%2 ? -1. : 1.)*(i==0 || i==count-1 ? .5 : 1.);
|
||||
for (int j = 0; j < count; ++j) {
|
||||
if (i == j) { continue; }
|
||||
const double wj = (j%2 ? -1. : 1.)*(j==0 || j==count-1 ? .5 : 1.);
|
||||
result[i][j] = (wj/wi)/(zeta[i]-zeta[j]);
|
||||
}
|
||||
}
|
||||
return result;
|
||||
}
|
||||
Matrix differentiate(const Matrix& values, const Matrix& derivative) {
|
||||
Matrix result(values.size(), std::vector<double>(values.front().size()));
|
||||
for (std::size_t i = 0; i < values.size(); ++i) {
|
||||
for (std::size_t l = 0; l < values[i].size(); l += 2) {
|
||||
long double sum = 0.;
|
||||
for (std::size_t j = 0; j < values.size(); ++j) {
|
||||
sum += derivative[i][j]*(values[j][l]-values[i][l]);
|
||||
}
|
||||
result[i][l] = static_cast<double>(sum);
|
||||
if (!std::isfinite(result[i][l])) { fail("enthalpy_et_der", "A radial spectral derivative is nonfinite."); }
|
||||
}
|
||||
}
|
||||
return result;
|
||||
}
|
||||
struct TemporaryDirectory {
|
||||
std::filesystem::path path;
|
||||
~TemporaryDirectory() { if (!path.empty()) { std::error_code error; std::filesystem::remove_all(path, error); } }
|
||||
};
|
||||
std::ofstream open_output(const std::filesystem::path& path) {
|
||||
std::ofstream output(path);
|
||||
output.imbue(std::locale::classic());
|
||||
output << std::scientific << std::setprecision(std::numeric_limits<double>::max_digits10);
|
||||
if (!output) { fail(path.string(), "Could not create TOP output file."); }
|
||||
return output;
|
||||
}
|
||||
void close_output(std::ofstream& output, const std::filesystem::path& path) {
|
||||
output.close();
|
||||
if (!output) { fail(path.string(), "Could not finish writing TOP output file."); }
|
||||
}
|
||||
/** @brief Unlike rename(), these publication primitives cannot replace an existing empty directory. */
|
||||
void publish(const std::filesystem::path& temporary, const std::filesystem::path& path) {
|
||||
#if defined(__APPLE__)
|
||||
if (::renamex_np(temporary.c_str(), path.c_str(), RENAME_EXCL) == 0) { return; }
|
||||
fail(path.string(), "Could not publish TOP directory without replacement: " + std::string{std::strerror(errno)});
|
||||
#elif defined(__linux__) && defined(SYS_renameat2)
|
||||
if (::syscall(SYS_renameat2, AT_FDCWD, temporary.c_str(), AT_FDCWD, path.c_str(), 1U) == 0) { return; }
|
||||
fail(path.string(), "Could not publish TOP directory without replacement: " + std::string{std::strerror(errno)});
|
||||
#else
|
||||
fail(path.string(), "Atomic directory publication without replacement is unavailable on this platform.");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
Result<TopExportReport> export_top(const std::filesystem::path& requested_path,
|
||||
const TopPolytropeSource& source, const TopExportOptions& options) {
|
||||
try {
|
||||
if (requested_path.empty() || requested_path.filename().empty() || requested_path.filename()=="." || requested_path.filename()=="..") {
|
||||
fail(requested_path.string(), "Specify a new model directory, not an empty path or parent directory.");
|
||||
}
|
||||
const auto path = std::filesystem::absolute(requested_path).lexically_normal();
|
||||
if (std::filesystem::symlink_status(path).type() != std::filesystem::file_type::not_found) {
|
||||
fail(path.string(), "TOP export refuses to replace an existing destination.");
|
||||
}
|
||||
if (!std::filesystem::is_directory(path.parent_path())) { fail(path.string(), "The output parent directory does not exist."); }
|
||||
if (options.radial_points < 4 || options.radial_points > 4096) { fail("options/radial_points", "Use between 4 and 4096 radial points."); }
|
||||
if (options.maximum_degree < 0 || options.maximum_degree > 512 || options.maximum_degree%2) {
|
||||
fail("options/maximum_degree", "The maximum harmonic degree must be even and between 0 and 512.");
|
||||
}
|
||||
if (options.angular_points < std::max(2, options.maximum_degree+1) || options.angular_points > 4096) {
|
||||
fail("options/angular_points", "Angular quadrature requires at least maximum_degree+1 points and between 2 and 4096 points.");
|
||||
}
|
||||
if (options.azimuthal_points < 4 || options.azimuthal_points > 1024) { fail("options/azimuthal_points", "Use between 4 and 1024 azimuthal samples."); }
|
||||
for (const auto& [name, value] : {std::pair{"relative_symmetry_tolerance", options.relative_symmetry_tolerance},
|
||||
{"relative_surface_tolerance", options.relative_surface_tolerance},
|
||||
{"relative_surface_projection_tolerance", options.relative_surface_projection_tolerance}}) {
|
||||
if (!std::isfinite(value) || value < 0.) { fail("options/" + std::string{name}, "Tolerance must be finite and nonnegative."); }
|
||||
}
|
||||
for (const auto& [name, value] : {std::pair{"polytropic_index", source.polytropic_index},
|
||||
{"gravitational_constant", source.gravitational_constant}, {"mass", source.mass},
|
||||
{"central_density", source.central_density}, {"central_enthalpy", source.central_enthalpy}}) {
|
||||
if (!positive(value)) { fail("source/" + std::string{name}, "This physical scale must be finite and positive."); }
|
||||
}
|
||||
if (source.polytropic_index < 1.) { fail("source/polytropic_index", "TOP's vacuum polytrope equations require n>=1 in this exporter."); }
|
||||
if (!std::isfinite(source.angular_velocity)) { fail("source/angular_velocity", "Angular velocity must be finite."); }
|
||||
if (!source.surface_radius || !source.enthalpy) { fail("source", "Both surface-radius and enthalpy sampling callbacks are required."); }
|
||||
TopExportReport report{.directory=path, .radial_points=options.radial_points, .maximum_degree=options.maximum_degree,
|
||||
.angular_points=options.angular_points, .azimuthal_points=options.azimuthal_points};
|
||||
const auto [mu, weights] = gauss(options.angular_points);
|
||||
auto radii = [&](double cosine) {
|
||||
std::vector<double> values;
|
||||
values.reserve(2*options.azimuthal_points);
|
||||
for (double sign : {1., -1.}) {
|
||||
for (int k = 0; k < options.azimuthal_points; ++k) {
|
||||
const double phi = 2*pi*k/options.azimuthal_points;
|
||||
const double value = checked(source.surface_radius(sign*cosine, phi), "source/surface_radius");
|
||||
if (value <= 0.) { fail("source/surface_radius", "The surface must be star-shaped with positive radius in every sampled direction."); }
|
||||
values.push_back(value);
|
||||
}
|
||||
}
|
||||
return values;
|
||||
};
|
||||
auto mean = [](const std::vector<double>& values) {
|
||||
long double sum = 0.; for (double value : values) { sum += value; }
|
||||
return static_cast<double>(sum/values.size());
|
||||
};
|
||||
report.source_equatorial_radius = mean(radii(0.));
|
||||
report.source_polar_radius = mean(radii(1.));
|
||||
auto surface_samples = [&](double cosine) {
|
||||
auto values = radii(cosine);
|
||||
const double average = mean(values);
|
||||
for (double value : values) {
|
||||
report.maximum_relative_surface_symmetry_error = std::max(report.maximum_relative_surface_symmetry_error,
|
||||
std::abs(value-average)/report.source_equatorial_radius);
|
||||
}
|
||||
return values;
|
||||
};
|
||||
Matrix radius_values;
|
||||
std::vector<double> surface_means;
|
||||
for (double cosine : mu) {
|
||||
radius_values.push_back(surface_samples(cosine));
|
||||
surface_means.push_back(mean(radius_values.back()));
|
||||
}
|
||||
auto surface_coefficients = project(surface_means, mu, weights, options.maximum_degree);
|
||||
report.equatorial_radius = reconstruct(surface_coefficients, 0.);
|
||||
report.polar_radius = reconstruct(surface_coefficients, 1.);
|
||||
if (!positive(report.equatorial_radius) || !positive(report.polar_radius)) { fail("domains_boundaries", "Projected equatorial and polar radii must be positive."); }
|
||||
report.minimum_mapping_derivative = std::numeric_limits<double>::infinity();
|
||||
auto mapping_check = [&](double radius) {
|
||||
// A'(zeta) ranges from zero to 15/8, so this is the exact radial minimum.
|
||||
const double minimum = std::min(report.polar_radius,
|
||||
report.polar_radius + 1.875*(radius-report.polar_radius))/report.equatorial_radius;
|
||||
report.minimum_mapping_derivative = std::min(report.minimum_mapping_derivative, minimum);
|
||||
if (!positive(minimum)) { fail("domains_boundaries", "The source or projected surface gives a nonpositive TOP mapping derivative."); }
|
||||
};
|
||||
for (int j = 0; j <= 2*options.angular_points; ++j) {
|
||||
const double cosine = std::cos(pi*j/(2*options.angular_points));
|
||||
const auto values = surface_samples(cosine);
|
||||
const double projected = reconstruct(surface_coefficients, cosine);
|
||||
if (!positive(projected)) { fail("domains_boundaries", "Projected surface radius must stay positive."); }
|
||||
mapping_check(projected);
|
||||
for (double value : values) { mapping_check(value); }
|
||||
report.maximum_relative_surface_projection_error = std::max(report.maximum_relative_surface_projection_error,
|
||||
std::abs(projected-mean(values))/report.source_equatorial_radius);
|
||||
}
|
||||
if (report.maximum_relative_surface_symmetry_error > options.relative_symmetry_tolerance) {
|
||||
fail("source/surface_radius", "Surface axisymmetry or equatorial-symmetry defect exceeds relative_symmetry_tolerance.");
|
||||
}
|
||||
if (report.maximum_relative_surface_projection_error > options.relative_surface_projection_tolerance) {
|
||||
fail("domains_boundaries", "Surface harmonic truncation exceeds relative_surface_projection_tolerance; increase maximum_degree or angular_points.");
|
||||
}
|
||||
report.flattening = 1.-report.polar_radius/report.equatorial_radius;
|
||||
report.angular_frequency_unit = std::sqrt(4*pi*source.gravitational_constant*source.central_density);
|
||||
report.lambda = std::pow(report.angular_frequency_unit*report.equatorial_radius, 2)/source.central_enthalpy;
|
||||
report.alpha = 4*pi*source.central_density*std::pow(report.equatorial_radius, 3)/(3*source.mass);
|
||||
report.rotation = source.angular_velocity*report.equatorial_radius/std::sqrt(source.central_enthalpy);
|
||||
report.rotation_over_kepler = source.angular_velocity*std::sqrt(std::pow(report.equatorial_radius, 3)/(source.gravitational_constant*source.mass));
|
||||
if (!positive(report.angular_frequency_unit) || !positive(report.lambda) || !positive(report.alpha)
|
||||
|| !std::isfinite(report.rotation) || !std::isfinite(report.rotation_over_kepler)) {
|
||||
fail("lambda_eps", "Physical normalization overflowed or underflowed.");
|
||||
}
|
||||
std::vector<double> zeta(options.radial_points);
|
||||
for (int i = 0; i < options.radial_points; ++i) { zeta[i] = .5*(1-std::cos(pi*i/(options.radial_points-1))); }
|
||||
zeta.front() = 0.; zeta.back() = 1.;
|
||||
report.minimum_enthalpy = std::numeric_limits<double>::infinity();
|
||||
auto enthalpy_samples = [&](int i, double cosine, const std::vector<double>& surfaces) {
|
||||
const double z = zeta[i], a = (5*z*z*z-3*std::pow(z, 5))/2;
|
||||
std::vector<double> values;
|
||||
values.reserve(surfaces.size());
|
||||
std::size_t index = 0;
|
||||
for (double sign : {1., -1.}) {
|
||||
for (int k = 0; k < options.azimuthal_points; ++k) {
|
||||
const double radius = i==options.radial_points-1 ? surfaces[index] :
|
||||
report.polar_radius*z+a*(surfaces[index]-report.polar_radius);
|
||||
const double raw = checked(source.enthalpy(radius, sign*cosine, 2*pi*k/options.azimuthal_points),
|
||||
"source/enthalpy/radial_" + std::to_string(i))/source.central_enthalpy;
|
||||
if (!std::isfinite(raw)) { fail("source/enthalpy", "Normalized enthalpy is nonfinite."); }
|
||||
if (i==options.radial_points-1) {
|
||||
report.maximum_relative_surface_enthalpy = std::max(report.maximum_relative_surface_enthalpy, std::abs(raw));
|
||||
if (std::abs(raw) > options.relative_surface_tolerance) {
|
||||
fail("source/enthalpy/surface", "TOP's vacuum surface requires zero enthalpy within relative_surface_tolerance.");
|
||||
}
|
||||
} else if (raw <= 0.) {
|
||||
fail("source/enthalpy/radial_"+std::to_string(i), "TOP requires strictly positive interior enthalpy for its logarithmic derivatives.");
|
||||
}
|
||||
if (i==0 && std::abs(raw-1.) > 1e-8) { fail("source/central_enthalpy", "The supplied central enthalpy does not match the source center."); }
|
||||
values.push_back(raw);
|
||||
++index;
|
||||
}
|
||||
}
|
||||
const double average = mean(values);
|
||||
for (double value : values) {
|
||||
report.maximum_relative_enthalpy_symmetry_error = std::max(report.maximum_relative_enthalpy_symmetry_error, std::abs(value-average));
|
||||
}
|
||||
if (report.maximum_relative_enthalpy_symmetry_error > options.relative_symmetry_tolerance) {
|
||||
fail("source/enthalpy", "Enthalpy axisymmetry or equatorial-symmetry defect exceeds relative_symmetry_tolerance.");
|
||||
}
|
||||
// The vacuum boundary is imposed only after measuring and accepting its source defect.
|
||||
const double result = i==options.radial_points-1 ? 0. : average;
|
||||
report.minimum_enthalpy = std::min(report.minimum_enthalpy, result);
|
||||
return result;
|
||||
};
|
||||
Matrix coefficients;
|
||||
coefficients.reserve(options.radial_points);
|
||||
for (int i = 0; i < options.radial_points; ++i) {
|
||||
std::vector<double> values;
|
||||
for (int j = 0; j < options.angular_points; ++j) { values.push_back(enthalpy_samples(i, mu[j], radius_values[j])); }
|
||||
coefficients.push_back(project(values, mu, weights, options.maximum_degree));
|
||||
for (int j = 0; j < options.angular_points; ++j) {
|
||||
report.maximum_relative_enthalpy_projection_error = std::max(report.maximum_relative_enthalpy_projection_error,
|
||||
std::abs(reconstruct(coefficients.back(), mu[j])-values[j]));
|
||||
}
|
||||
}
|
||||
// Include poles and equator in physical boundary/symmetry checks, independently of Gauss nodes.
|
||||
for (double cosine : {0., 1.}) {
|
||||
const auto surfaces = surface_samples(cosine);
|
||||
for (int i = 0; i < options.radial_points; ++i) { enthalpy_samples(i, cosine, surfaces); }
|
||||
}
|
||||
const auto derivative = derivative_matrix(zeta);
|
||||
const auto first = differentiate(coefficients, derivative);
|
||||
const auto second = differentiate(first, derivative);
|
||||
// TOP uses h_zeta_theta/h_zeta at the surface, so a flat or outward-increasing
|
||||
// exported profile cannot be handed to that reader even when all samples are finite.
|
||||
for (int j = 0; j <= 2*options.angular_points; ++j) {
|
||||
const double slope = reconstruct(first.back(), std::cos(pi*j/(2*options.angular_points)));
|
||||
if (!std::isfinite(slope) || slope >= -128*std::numeric_limits<double>::epsilon()) {
|
||||
fail("enthalpy_et_der/surface_derivative", "TOP requires strictly negative, nonvanishing surface enthalpy derivative; refine the export or source model.");
|
||||
}
|
||||
}
|
||||
// TOP takes abs(h) silently and divides by h. Reject nonpositive projection values before either operation.
|
||||
for (int i = 0; i < options.radial_points-1; ++i) {
|
||||
for (int j = 0; j <= 2*options.angular_points; ++j) {
|
||||
const double value = reconstruct(coefficients[i], std::cos(pi*j/(2*options.angular_points)));
|
||||
if (!positive(value)) {
|
||||
fail("enthalpy_et_der", "Harmonic projection creates nonpositive interior enthalpy; increase angular resolution.");
|
||||
}
|
||||
}
|
||||
}
|
||||
if (source.roundoff_density_clamps) { report.roundoff_density_clamps = source.roundoff_density_clamps(); }
|
||||
for (double& coefficient : surface_coefficients) { coefficient /= report.equatorial_radius; }
|
||||
TemporaryDirectory temporary;
|
||||
static std::atomic<unsigned long long> sequence{0};
|
||||
for (int attempt = 0; attempt < 100; ++attempt) {
|
||||
const auto suffix = std::to_string(std::chrono::steady_clock::now().time_since_epoch().count()) + "-" + std::to_string(sequence++);
|
||||
auto candidate = path.parent_path()/("."+path.filename().string()+".top-tmp-"+suffix);
|
||||
if (std::filesystem::create_directory(candidate)) { temporary.path=std::move(candidate); break; }
|
||||
}
|
||||
if (temporary.path.empty()) { fail(path.string(), "Could not reserve a temporary TOP directory."); }
|
||||
{
|
||||
const auto file = temporary.path/"lambda_eps";
|
||||
auto output = open_output(file);
|
||||
for (double value : {report.lambda, report.flattening, report.alpha, 0., source.polytropic_index,
|
||||
source.polytropic_index, report.rotation, report.rotation_over_kepler}) { output << value << '\n'; }
|
||||
output << options.radial_points << '\n';
|
||||
close_output(output, file);
|
||||
}
|
||||
{
|
||||
const auto file = temporary.path/"domains_boundaries";
|
||||
auto output = open_output(file);
|
||||
output << options.maximum_degree << '\n';
|
||||
for (int l = 0; l <= options.maximum_degree; ++l) { output << 0. << ' '; } output << '\n';
|
||||
for (double value : surface_coefficients) { output << value << ' '; } output << '\n';
|
||||
close_output(output, file);
|
||||
}
|
||||
{
|
||||
const auto file = temporary.path/"enthalpy_et_der";
|
||||
auto output = open_output(file);
|
||||
output << options.radial_points << ' ' << options.maximum_degree+1 << '\n';
|
||||
output << "enthalpy in physical space: (one hemisphere)\n";
|
||||
const auto physical_mu = gauss(options.maximum_degree+1).first;
|
||||
for (int j = 0; j <= options.maximum_degree/2; ++j) {
|
||||
output << j+1 << '\n';
|
||||
for (const auto& row : coefficients) { output << reconstruct(row, physical_mu[j]) << '\n'; }
|
||||
}
|
||||
auto block = [&](const char* label, const Matrix& values) {
|
||||
output << label << '\n';
|
||||
for (int l = 0; l <= options.maximum_degree; l += 2) {
|
||||
output << l << '\n'; for (const auto& row : values) { output << row[l] << '\n'; }
|
||||
}
|
||||
};
|
||||
block("enthalpy in legendre space: l=0,2,4...lmax", coefficients);
|
||||
block("first radial derivative in legendre space", first);
|
||||
block("second radial derivative in legendre space", second);
|
||||
close_output(output, file);
|
||||
}
|
||||
{
|
||||
const auto file = temporary.path/"metadata.txt";
|
||||
auto output = open_output(file);
|
||||
output << "SERiF TOP poly_ester export version 1\n"
|
||||
<< "Coordinates: zeta=(1-cos(pi*i/(nr-1)))/2; mu=cos(theta) relative to the source rotation axis.\n"
|
||||
<< "Basis: Y_l0=sqrt((2*l+1)/(4*pi))*P_l(mu); even l only.\n"
|
||||
<< "Projection: source surface and enthalpy are averaged over azimuth and both hemispheres, then projected at common surface-fitted zeta.\n"
|
||||
<< "Physical sampling map: r=Rp*zeta+(5*zeta^3-3*zeta^5)/2*(source_surface_radius-Rp).\n"
|
||||
<< "TOP replaces source_surface_radius with its harmonic projection; this is a recorded geometry approximation, not an exact FE remap.\n"
|
||||
<< "Derivatives: Chebyshev interpolation derivatives at fixed angle; export refinement is required before interpreting mode accuracy.\n"
|
||||
<< "Surface: accepted source enthalpy defect is recorded, then boundary enthalpy is set to zero. Interior values are not clipped.\n"
|
||||
<< "Enthalpy projection diagnostic: maximum Gauss-node error of the averaged H/Hc after vacuum-boundary normalization; the source surface defect is reported separately.\n"
|
||||
<< "Diagnostic units: enthalpy and its defects are normalized by central_enthalpy; radial mapping derivative by equatorial_radius.\n"
|
||||
<< "Frequency: TOP angular frequency times angular_frequency_unit gives source-unit angular frequency; divide by 2*pi for cyclic frequency.\n"
|
||||
<< "Gamma1 is selected by TOP's equation set, not by these files. Diagnostic maxima are sampled, not rigorous global bounds.\n";
|
||||
#define SERIF_TOP_METADATA(member) output << #member << " " << report.member << '\n'
|
||||
SERIF_TOP_METADATA(radial_points); SERIF_TOP_METADATA(maximum_degree); SERIF_TOP_METADATA(angular_points); SERIF_TOP_METADATA(azimuthal_points);
|
||||
SERIF_TOP_METADATA(source_equatorial_radius); SERIF_TOP_METADATA(source_polar_radius);
|
||||
SERIF_TOP_METADATA(equatorial_radius); SERIF_TOP_METADATA(polar_radius); SERIF_TOP_METADATA(lambda); SERIF_TOP_METADATA(flattening);
|
||||
SERIF_TOP_METADATA(alpha); SERIF_TOP_METADATA(rotation); SERIF_TOP_METADATA(rotation_over_kepler); SERIF_TOP_METADATA(angular_frequency_unit);
|
||||
SERIF_TOP_METADATA(maximum_relative_surface_symmetry_error); SERIF_TOP_METADATA(maximum_relative_enthalpy_symmetry_error);
|
||||
SERIF_TOP_METADATA(maximum_relative_surface_projection_error); SERIF_TOP_METADATA(maximum_relative_enthalpy_projection_error);
|
||||
SERIF_TOP_METADATA(maximum_relative_surface_enthalpy);
|
||||
SERIF_TOP_METADATA(minimum_mapping_derivative); SERIF_TOP_METADATA(minimum_enthalpy); SERIF_TOP_METADATA(roundoff_density_clamps);
|
||||
#undef SERIF_TOP_METADATA
|
||||
output << "polytropic_index " << source.polytropic_index << '\n' << "gravitational_constant " << source.gravitational_constant << '\n'
|
||||
<< "mass " << source.mass << '\n' << "angular_velocity " << source.angular_velocity << '\n'
|
||||
<< "central_density " << source.central_density << '\n' << "central_enthalpy " << source.central_enthalpy << '\n'
|
||||
<< "relative_symmetry_tolerance " << options.relative_symmetry_tolerance << '\n'
|
||||
<< "relative_surface_tolerance " << options.relative_surface_tolerance << '\n'
|
||||
<< "relative_surface_projection_tolerance " << options.relative_surface_projection_tolerance << '\n';
|
||||
close_output(output, file);
|
||||
}
|
||||
publish(temporary.path, path);
|
||||
temporary.path.clear();
|
||||
return report;
|
||||
} catch (const Failure& failure) {
|
||||
return std::unexpected(failure.error);
|
||||
} catch (const std::exception& error) {
|
||||
return std::unexpected(Error{requested_path.string(), error.what()});
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user