tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(pythonInterface/eos): fast forward

Browse Source
This commit is contained in:
Emily Boudreaux
2025-06-12 14:04:11 -04:00
parent b5980ea57a 353e2072e7
commit 54e408a4af
273 changed files with 42783 additions and 12196 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
src/composition/private/composition.cpp
View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "composition.h"
#include "quill/LogMacros.h"
@@ -10,7 +30,7 @@
#include "atomicSpecies.h"
using namespace composition;
namespace serif::composition {
CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")), m_initialized(false) {}
@@ -527,10 +547,24 @@ Composition Composition::operator+(const Composition& other) const {
return mix(other, 0.5);
}
std::ostream& composition::operator<<(std::ostream& os, const Composition& composition) {
std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp) {
os << "Global Composition: \n";
os << "\tSpecific Number Density: " << comp.specificNumberDensity << "\n";
os << "\tMean Particle Mass: " << comp.meanParticleMass << "\n";
return os;
}
std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry) {
os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">";
return os;
}
std::ostream& operator<<(std::ostream& os, const Composition& composition) {
os << "Composition: \n";
for (const auto& [symbol, entry] : composition.m_compositions) {
os << entry << "\n";
}
return os;
}
}
} // namespace serif::composition

39
src/composition/public/composition.h
View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef COMPOSITION_H
#define COMPOSITION_H
@@ -13,7 +33,7 @@
#include "atomicSpecies.h"
namespace composition{
namespace serif::composition {
/**
* @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use.
*/
@@ -22,12 +42,7 @@ namespace composition{
double meanParticleMass; ///< The mean particle mass of the composition (\sum_{i} \frac{n_i}{m_i}. where n_i is the number fraction of the ith species and m_i is the mass of the ith species).
// Overload the output stream operator for GlobalComposition
friend std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp) {
os << "Global Composition: \n";
os << "\tSpecific Number Density: " << comp.specificNumberDensity << "\n";
os << "\tMean Particle Mass: " << comp.meanParticleMass << "\n";
return os;
}
friend std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp);
};
/**
@@ -154,10 +169,7 @@ namespace composition{
* @param entry The CompositionEntry to output.
* @return The output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry) {
os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">";
return os;
}
friend std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry);
};
/**
@@ -190,7 +202,8 @@ namespace composition{
*/
class Composition {
private:
Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
bool m_finalized = false; ///< True if the composition is finalized.
@@ -450,6 +463,6 @@ namespace composition{
Composition operator+(const Composition& other) const;
};
};
}; // namespace serif::composition
#endif // COMPOSITION_H

134
src/config/private/config.cpp
View File

@@ -2,7 +2,7 @@
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 21, 2025
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -29,77 +29,81 @@
#include "config.h"
Config::Config() {}
namespace serif::config {
Config::~Config() {}
Config::Config() {}
Config& Config::getInstance() {
static Config instance;
return instance;
}
Config::~Config() {}
bool Config::loadConfig(const std::string& configFile) {
configFilePath = configFile;
try {
yamlRoot = YAML::LoadFile(configFile);
} catch (YAML::BadFile& e) {
std::cerr << "Error: " << e.what() << std::endl;
return false;
Config& Config::getInstance() {
static Config instance;
return instance;
}
m_loaded = true;
return true;
}
bool Config::isKeyInCache(const std::string &key) {
return configMap.find(key) != configMap.end();
}
void Config::addToCache(const std::string &key, const YAML::Node &node) {
configMap[key] = node;
}
void Config::registerUnknownKey(const std::string &key) {
unknownKeys.push_back(key);
}
bool Config::has(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
}
if (isKeyInCache(key)) { return true; }
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
registerUnknownKey(key);
bool Config::loadConfig(const std::string& configFile) {
configFilePath = configFile;
try {
yamlRoot = YAML::LoadFile(configFile);
} catch (YAML::BadFile& e) {
std::cerr << "Error: " << e.what() << std::endl;
return false;
}
node = node[subKey]; // go deeper
}
// Key exists and is of the requested type
addToCache(key, node);
return true;
}
void recurse_keys(const YAML::Node& node, std::vector<std::string>& keyList, const std::string& path = "") {
if (node.IsMap()) {
for (const auto& it : node) {
std::string key = it.first.as<std::string>();
std::string new_path = path.empty() ? key : path + ":" + key;
recurse_keys(it.second, keyList, new_path);
}
} else {
keyList.push_back(path);
m_loaded = true;
return true;
}
bool Config::isKeyInCache(const std::string &key) {
return configMap.find(key) != configMap.end();
}
void Config::addToCache(const std::string &key, const YAML::Node &node) {
configMap[key] = node;
}
void Config::registerUnknownKey(const std::string &key) {
unknownKeys.push_back(key);
}
bool Config::has(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
}
if (isKeyInCache(key)) { return true; }
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
registerUnknownKey(key);
return false;
}
node = node[subKey]; // go deeper
}
// Key exists and is of the requested type
addToCache(key, node);
return true;
}
void recurse_keys(const YAML::Node& node, std::vector<std::string>& keyList, const std::string& path = "") {
if (node.IsMap()) {
for (const auto& it : node) {
auto key = it.first.as<std::string>();
auto new_path = path.empty() ? key : path + ":" + key;
recurse_keys(it.second, keyList, new_path);
}
} else {
keyList.push_back(path);
}
}
std::vector<std::string> Config::keys() const {
std::vector<std::string> keyList;
YAML::Node node = YAML::Clone(yamlRoot);
recurse_keys(node, keyList);
return keyList;
}
}
std::vector<std::string> Config::keys() const {
std::vector<std::string> keyList;
YAML::Node node = YAML::Clone(yamlRoot);
recurse_keys(node, keyList);
return keyList;
}

223
src/config/public/config.h
View File

@@ -2,7 +2,7 @@
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 20, 2025
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -32,201 +32,214 @@
#include "yaml-cpp/yaml.h"
// -- Forward Def of Resource manager to let it act as a friend of Config --
class ResourceManager;
namespace serif::resource { class ResourceManager; }
/**
class configTestPrivateAccessor; // Forward declaration for test utility
namespace serif::config {
/**
* @class Config
* @brief Singleton class to manage configuration settings loaded from a YAML file.
*/
class Config {
private:
/**
class Config {
private:
/**
* @brief Private constructor to prevent instantiation.
*/
Config();
Config();
/**
/**
* @brief Destructor.
*/
~Config();
~Config();
YAML::Node yamlRoot; ///< Root node of the YAML configuration.
std::string configFilePath; ///< Path to the configuration file.
bool debug = false; ///< Flag to enable debug output.
YAML::Node yamlRoot; ///< Root node of the YAML configuration.
std::string configFilePath; ///< Path to the configuration file.
bool debug = false; ///< Flag to enable debug output.
bool loaded = false; ///< Flag to indicate if the configuration has been loaded.
std::map<std::string, YAML::Node> configMap; ///< Cache for the location of configuration settings.
std::vector<std::string> unknownKeys; ///< Cache for the existence of configuration settings.
std::map<std::string, YAML::Node> configMap; ///< Cache for the location of configuration settings.
std::vector<std::string> unknownKeys; ///< Cache for the existence of configuration settings.
/**
/**
* @brief Get a value from the configuration cache.
* @tparam T Type of the value to retrieve.
* @param key Key of the configuration value.
* @param defaultValue Default value to return if the key does not exist.
* @return Configuration value of type T.
*/
template <typename T>
T getFromCache(const std::string &key, T defaultValue) {
if (configMap.find(key) != configMap.end()) {
try {
template <typename T>
T getFromCache(const std::string &key, T defaultValue) {
if (configMap.find(key) != configMap.end()) {
try {
return configMap[key].as<T>();
} catch (const YAML::Exception& e) {
} catch (const YAML::Exception& e) {
return defaultValue;
}
}
return defaultValue;
}
return defaultValue;
}
/**
/**
* @brief Check if a key exists in the configuration cache.
* @param key Key to check.
* @return True if the key exists in the cache, false otherwise.
*/
bool isKeyInCache(const std::string &key);
bool isKeyInCache(const std::string &key);
/**
/**
* @brief Add a key-value pair to the configuration cache.
* @param key Key of the configuration value.
* @param node YAML node containing the configuration value.
*/
void addToCache(const std::string &key, const YAML::Node &node);
void addToCache(const std::string &key, const YAML::Node &node);
/**
/**
* @brief Register a key as not found in the configuration.
* @param key Key that was not found.
*/
void registerUnknownKey(const std::string &key);
void registerUnknownKey(const std::string &key);
bool m_loaded = false;
bool m_loaded = false;
// Only friends can access get without a default value
template <typename T>
T get(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
// Only friends can access get without a default value
template <typename T>
T get(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
}
if (has(key)) {
return getFromCache<T>(key, T());
} else {
throw std::runtime_error("Error! Key not found in config file");
}
}
if (has(key)) {
return getFromCache<T>(key, T());
} else {
throw std::runtime_error("Error! Key not found in config file");
}
}
public:
/**
public:
/**
* @brief Get the singleton instance of the Config class.
* @return Reference to the Config instance.
*/
static Config& getInstance();
static Config& getInstance();
Config (const Config&) = delete;
Config& operator= (const Config&) = delete;
Config (Config&&) = delete;
Config& operator= (Config&&) = delete;
Config (const Config&) = delete;
Config& operator= (const Config&) = delete;
Config (Config&&) = delete;
Config& operator= (Config&&) = delete;
void setDebug(bool debug) { this->debug = debug; }
void setDebug(bool debug) { this->debug = debug; }
/**
/**
* @brief Load configuration from a YAML file.
* @param configFilePath Path to the YAML configuration file.
* @return True if the configuration was loaded successfully, false otherwise.
*/
bool loadConfig(const std::string& configFilePath);
bool loadConfig(const std::string& configFilePath);
/**
/**
* @brief Get the input table from the configuration.
* @return Input table as a string.
*/
std::string getInputTable() const;
/**
* @brief Get a configuration value by key.
* @tparam T Type of the value to retrieve.
* @param key Key of the configuration value.
* @param defaultValue Default value to return if the key does not exist.
* @return Configuration value of type T.
*
*
* @example
* @code
* Config& config = Config::getInstance();
* config.loadConfig("example.yaml");
* int maxIter = config.get<int>("opac:lowTemp:numeric:maxIter", 10);
*/
template <typename T>
T get(const std::string &key, T defaultValue) {
if (!m_loaded) {
// ONLY THROW ERROR IF HARSH OR WARN CONFIGURATION
template <typename T>
T get(const std::string &key, T defaultValue) {
if (!m_loaded) {
// ONLY THROW ERROR IF HARSH OR WARN CONFIGURATION
#if defined(CONFIG_HARSH)
throw std::runtime_error("Error! Config file not loaded. To disable this error, recompile with CONFIG_HARSH=0");
#elif defined(CONFIG_WARN)
std::cerr << "Warning! Config file not loaded. This instance of 4DSSE was compiled with CONFIG_WARN so the code will continue using only default values" << std::endl;
#endif
}
// --- Check if the key has already been checked for existence
if (std::find(unknownKeys.begin(), unknownKeys.end(), key) != unknownKeys.end()) {
return defaultValue; // If the key has already been added to the unknown cache do not traverse the YAML tree or hit the cache
}
}
// --- Check if the key has already been checked for existence
if (std::find(unknownKeys.begin(), unknownKeys.end(), key) != unknownKeys.end()) {
return defaultValue; // If the key has already been added to the unknown cache do not traverse the YAML tree or hit the cache
}
// --- Check if the key is already in the cache (avoid traversing YAML nodes)
if (isKeyInCache(key)) {
return getFromCache<T>(key, defaultValue);
}
// --- If the key is not in the cache, check the YAML file
else {
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
// --- Check if the key is already in the cache (avoid traversing YAML nodes)
if (isKeyInCache(key)) {
return getFromCache<T>(key, defaultValue);
}
// --- If the key is not in the cache, check the YAML file
else {
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
// Key does not exist
registerUnknownKey(key);
return defaultValue;
}
node = node[subKey]; // go deeper
}
node = node[subKey]; // go deeper
}
try {
// Key exists and is of the requested type
addToCache(key, node);
return node.as<T>();
} catch (const YAML::Exception& e) {
// Key is not of the requested type
registerUnknownKey(key);
return defaultValue; // return default value if the key does not exist
try {
// Key exists and is of the requested type
addToCache(key, node);
return node.as<T>();
} catch (const YAML::Exception& e) {
// Key is not of the requested type
registerUnknownKey(key);
return defaultValue; // return default value if the key does not exist
}
}
}
}
/**
/**
* @brief Check if the key exists in the given config file
* @param key Key to check;
* @return boolean true or false
*/
bool has(const std::string &key);
bool has(const std::string &key);
/**
/**
* @brief Get all keys defined in the configuration file.
* @return Vector of all keys in the configuration file.
*/
std::vector<std::string> keys() const;
std::vector<std::string> keys() const;
/**
/**
* @brief Print the configuration file path and the YAML root node.
* @param os Output stream.
* @param config Config object to print.
* @return Output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const Config& config) {
if (!config.m_loaded) {
os << "Config file not loaded" << std::endl;
friend std::ostream& operator<<(std::ostream& os, const Config& config) {
if (!config.m_loaded) {
os << "Config file not loaded" << std::endl;
return os;
}
if (!config.debug) {
os << "Config file: " << config.configFilePath << std::endl;
} else{
// Print entire YAML file from root
os << "Config file: " << config.configFilePath << std::endl;
os << config.yamlRoot << std::endl;
}
return os;
}
if (!config.debug) {
os << "Config file: " << config.configFilePath << std::endl;
} else{
// Print entire YAML file from root
os << "Config file: " << config.configFilePath << std::endl;
os << config.yamlRoot << std::endl;
}
return os;
}
// Setup gTest class as a friend
friend class configTestPrivateAccessor;
// --- Resource Manager is a friend of config so it can create a separate instance
friend class ResourceManager;
};
// Setup gTest class as a friend
friend class ::configTestPrivateAccessor; // Friend declaration for global test accessor
// -- Resource Manager is a friend of config so it can create a seperate instance
friend class serif::resource::ResourceManager; // Adjusted friend declaration
};
}

124
src/const/private/const.cpp
View File

@@ -1,124 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 12, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <iostream>
#include <map>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>
#include <set>
#include "const.h"
#include "embedded_constants.h" // Generated at build time by meson
Constants::Constants() {
loaded_ = initialize();
}
bool Constants::initialize() {
return load();
}
Constant Constants::get(const std::string& name) const {
auto it = constants_.find(name);
if (it != constants_.end()) {
return it->second;
} else {
throw std::out_of_range("Constant '" + name + "' not found.");
}
}
Constant Constants::operator[](const std::string& name) const {
return this->get(name);
}
bool Constants::has(const std::string& name) const {
return constants_.find(name) != constants_.end();
}
std::set<std::string> Constants::keys() const {
std::set<std::string> keys;
for (const auto& pair : constants_) {
keys.insert(pair.first);
}
return keys;
}
std::string Constants::trim(const std::string& str) {
size_t first = str.find_first_not_of(" \t");
if (first == std::string::npos) return "";
size_t last = str.find_last_not_of(" \t");
return str.substr(first, last - first + 1);
}
bool Constants::load() {
std::istringstream fileStream(embeddedConstants);
std::string line;
bool data_section = false;
int line_count = 0;
while (std::getline(fileStream, line)) {
line_count++;
// Detect start of data section (double divider line)
if (!data_section) {
if (line.find("Symbol") != std::string::npos) { // Find header row
std::getline(fileStream, line); // Skip dashed divider
std::getline(fileStream, line); // Skip second dashed divider
data_section = true;
}
continue;
}
// Ensure the line is long enough to contain all fields
if (line.length() < 125) continue;
// Define exact column widths from Python script
int start = 0;
const std::string symbol = trim(line.substr(start, col_widths_[0])); start += col_widths_[0];
const std::string name = trim(line.substr(start, col_widths_[1])); start += col_widths_[1];
const std::string valueStr = line.substr(start, col_widths_[2]); start += col_widths_[2];
const std::string unit = trim(line.substr(start, col_widths_[3])); start += col_widths_[3]; // Only trim the unit
const std::string uncertaintyStr = line.substr(start, col_widths_[4]); start += col_widths_[4];
const std::string reference = trim(line.substr(start, col_widths_[5])); // Only trim reference
// Convert numerical fields safely
double value = 0.0, uncertainty = 0.0;
try {
value = std::stod(valueStr);
} catch (...) {
std::cerr << "Warning: Invalid value in line " << line_count << ": " << valueStr << std::endl;
}
try {
uncertainty = std::stod(uncertaintyStr);
} catch (...) {
std::cerr << "Warning: Invalid uncertainty in line " << line_count << ": " << uncertaintyStr << std::endl;
}
// Store in map
constants_.emplace(symbol, Constant{name, value, uncertainty, unit, reference});
}
loaded_ = true;
return true;
}

0
src/const/meson.build → src/constants/meson.build
View File

126
src/constants/private/const.cpp Normal file
View File

@@ -0,0 +1,126 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 17, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <iostream>
#include <map>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>
#include <set>
#include "const.h"
#include "embedded_constants.h" // Generated at build time by meson
namespace serif::constant {
Constants::Constants() {
loaded_ = initialize();
}
bool Constants::initialize() {
return load();
}
Constant Constants::get(const std::string& name) const {
auto it = constants_.find(name);
if (it != constants_.end()) {
return it->second;
} else {
throw std::out_of_range("Constant '" + name + "' not found.");
}
}
Constant Constants::operator[](const std::string& name) const {
return this->get(name);
}
bool Constants::has(const std::string& name) const {
return constants_.find(name) != constants_.end();
}
std::set<std::string> Constants::keys() const {
std::set<std::string> keys;
for (const auto& pair : constants_) {
keys.insert(pair.first);
}
return keys;
}
std::string Constants::trim(const std::string& str) {
size_t first = str.find_first_not_of(" \t");
if (first == std::string::npos) return "";
size_t last = str.find_last_not_of(" \t");
return str.substr(first, last - first + 1);
}
bool Constants::load() {
std::istringstream fileStream(embeddedConstants);
std::string line;
bool data_section = false;
int line_count = 0;
while (std::getline(fileStream, line)) {
line_count++;
// Detect start of data section (double divider line)
if (!data_section) {
if (line.find("Symbol") != std::string::npos) { // Find header row
std::getline(fileStream, line); // Skip dashed divider
std::getline(fileStream, line); // Skip second dashed divider
data_section = true;
}
continue;
}
// Ensure the line is long enough to contain all fields
if (line.length() < 125) continue;
// Define exact column widths from Python script
int start = 0;
const std::string symbol = trim(line.substr(start, col_widths_[0])); start += col_widths_[0];
const std::string name = trim(line.substr(start, col_widths_[1])); start += col_widths_[1];
const std::string valueStr = line.substr(start, col_widths_[2]); start += col_widths_[2];
const std::string unit = trim(line.substr(start, col_widths_[3])); start += col_widths_[3]; // Only trim the unit
const std::string uncertaintyStr = line.substr(start, col_widths_[4]); start += col_widths_[4];
const std::string reference = trim(line.substr(start, col_widths_[5])); // Only trim reference
// Convert numerical fields safely
double value = 0.0, uncertainty = 0.0;
try {
value = std::stod(valueStr);
} catch (...) {
std::cerr << "Warning: Invalid value in line " << line_count << ": " << valueStr << std::endl;
}
try {
uncertainty = std::stod(uncertaintyStr);
} catch (...) {
std::cerr << "Warning: Invalid uncertainty in line " << line_count << ": " << uncertaintyStr << std::endl;
}
// Store in map
constants_.emplace(symbol, Constant{name, value, uncertainty, unit, reference});
}
loaded_ = true;
return true;
}
}

7
src/const/public/const.h → src/constants/public/const.h
View File

@@ -2,7 +2,7 @@
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 05, 2025
// Last Modified: March 17, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -19,12 +19,12 @@
//
// *********************************************************************** */
#pragma once
#include <string>
#include <iostream>
#include <set>
#include <map>
namespace serif::constant {
/**
* @brief Structure to hold a constant's details.
*/
@@ -136,3 +136,6 @@ public:
std::set<std::string> keys() const;
};
} // namespace serif::const

22
src/dobj/meson.build
View File

@@ -1,22 +0,0 @@
# Define the library
dobj_sources = files(
'private/Metadata.cpp',
'private/DObject.cpp',
'private/LockableDObject.cpp'
)
dobj_headers = files(
'public/Metadata.h',
'public/DObject.h',
'public/LockableDObject.h'
)
# Define the libdobj library so it can be linked against by other parts of the build system
libdobj = library('dobj',
dobj_sources,
include_directories: include_directories('public'),
cpp_args: ['-fvisibility=default'],
install : true)
# Make headers accessible
install_headers(dobj_headers, subdir : '4DSSE/dobj')

239
src/dobj/private/DObject.cpp
View File

@@ -1,239 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 04, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "DObject.h"
#include <iostream>
#include <stdexcept>
/**
* @file DObject.cpp
* @brief Implementation of the DObject class.
*
* Provides the implementation for a universal data container with plugin support.
*/
// Default constructor
DObject::DObject() : data_(std::monostate()), debugEnabled_(false) {
// The data is initialized to an empty state using std::monostate (a variant placeholder).
}
/**
* @brief Constructor to initialize a DObject with data.
*/
DObject::DObject(const DataType& data)
: data_(data), debugEnabled_(false) {}
/**
* @brief Retrieves the data stored in the DObject.
*/
const DObject::DataType& DObject::getData() const noexcept {
return data_;
}
/**
* @brief Sets the data for the DObject.
*/
void DObject::setData(const DataType& data) {
data_ = data;
}
/**
* @brief Enables or disables debugging and tracing for the DObject.
*/
void DObject::setDebugging(bool enableDebug) {
debugEnabled_ = enableDebug;
}
/**
* @brief Checks if debugging is enabled for the DObject.
*/
bool DObject::isDebuggingEnabled() const noexcept {
return debugEnabled_;
}
/**
* @breif Sets an error code and returns the old one
*/
int DObject::setErrorCode(int code) noexcept {
int oldCode = errorCode_;
errorCode_ = code;
return oldCode;
}
int DObject::getErrorCode() const noexcept {
return errorCode_;
}
/**
* @brief Registers a plugin with the DObject.
*/
void DObject::registerPlugin(const std::string& id, Plugin plugin) {
if (plugins_.find(id) != plugins_.end()) {
throw std::runtime_error("Plugin with ID '" + id + "' already exists.");
}
plugins_[id] = std::move(plugin);
}
/**
* @brief Unregisters a plugin by its identifier.
*/
void DObject::unregisterPlugin(const std::string& id) {
if (plugins_.erase(id) == 0) {
throw std::runtime_error("Plugin with ID '" + id + "' not found.");
}
}
/**
* @brief Executes a plugin by its identifier.
*/
void DObject::runPlugin(const std::string& id) {
auto it = plugins_.find(id);
if (it == plugins_.end()) {
throw std::runtime_error("Plugin with ID '" + id + "' not found.");
}
it->second(*this); // Invoke the plugin, passing this DObject as a reference.
}
/**
* @brief Executes all registered plugins in the registry.
*/
void DObject::runAllPlugins() {
for (auto& [id, plugin] : plugins_) {
plugin(*this); // Invoke each plugin, passing this DObject as a reference.
}
}
/**
* @brief Provides a human-readable summary of the DObject.
*/
std::ostream& operator<<(std::ostream& os, const DObject& obj) {
os << "DObject Summary:\n";
os << " Debugging Enabled: " << (obj.debugEnabled_ ? "Yes" : "No") << "\n";
os << " Data: ";
std::visit([&os](const auto& value) {
using T = std::decay_t<decltype(value)>;
if constexpr (std::is_same_v<T, std::monostate>) {
os << "EMPTY";
} else if constexpr (std::is_same_v<T, int>) {
os << value;
} else if constexpr (std::is_same_v<T, float>) {
os << value;
} else if constexpr (std::is_same_v<T, double>) {
os << value;
} else if constexpr (std::is_same_v<T, std::string>) {
os << value;
} else if constexpr (std::is_same_v<T, std::vector<int>>) {
os << "[";
for (const auto& elem : value) {
os << elem << ", ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<float>>) {
os << "[";
for (const auto& elem : value) {
os << elem << ", ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<double>>) {
os << "[";
for (const auto& elem : value) {
os << elem << ", ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<int>>>) {
os << "\n[";
for (const auto& row : value) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<float>>>) {
os << "\n[";
for (const auto& row : value) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<double>>>) {
os << "\n[";
for (const auto& row : value) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<std::vector<int>>>>) {
os << "\n[";
for (const auto& plane : value) {
os << "[";
for (const auto& row : plane) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "], ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<std::vector<float>>>>) {
os << "\n[";
for (const auto& plane : value) {
os << "[";
for (const auto& row : plane) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "], ";
}
os << "]";
} else if constexpr (std::is_same_v<T, std::vector<std::vector<std::vector<double>>>>) {
os << "\n[";
for (const auto& plane : value) {
os << "[";
for (const auto& row : plane) {
os << "[";
for (const auto& elem : row) {
os << elem << ", ";
}
os << "], ";
}
os << "], ";
}
os << "]";
} else {
os << "UNKNOWN";
}
}, obj.data_);
os << "\n Data Type: " << obj.dataTypeMap.at(obj.data_.index());
os << "\n Plugins Registered: " << obj.plugins_.size() << "\n";
return os;
}

42
src/dobj/private/LockableDObject.cpp
View File

@@ -1,42 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: January 19, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "LockableDObject.h"
/**
* @brief Access the underlying DObject.
*/
DObject& LockableDObject::get() {
return object_;
}
/**
* @brief Locks the mutex to ensure thread-safe access.
*/
void LockableDObject::lock() {
mutex_.lock();
}
/**
* @brief Unlocks the mutex after thread-safe access.
*/
void LockableDObject::unlock() {
mutex_.unlock();
}

158
src/dobj/private/Metadata.cpp
View File

@@ -1,158 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: January 19, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "Metadata.h"
/**
* @file Metadata.cpp
* @brief Implementation of the Metadata class used in the dobj module.
*
* Provides methods to manage metadata for data objects, including size, type,
* dimensions, and debugging flags.
*/
/**
* @brief Constructor to initialize Metadata with specific attributes.
*
* @param byteSize The total size of the data in bytes.
* @param dataType The type of the data (e.g., "float", "double").
* @param dimensions A vector representing the size of each dimension (e.g., {3, 4} for a 3x4 matrix).
* @param debugFlag Whether debugging information is enabled for this Metadata instance.
*/
Metadata::Metadata(std::size_t byteSize, std::string dataType, std::vector<std::size_t> dimensions, bool debugFlag)
: byteSize_(byteSize), dataType_(std::move(dataType)), dimensions_(std::move(dimensions)), debugFlag_(debugFlag) {}
/**
* @brief Gets the total size of the data in bytes.
*
* The size is often required for memory allocation and validation in numerical routines.
*
* @return The total byte size of the data.
*/
std::size_t Metadata::getByteSize() const noexcept {
return byteSize_;
}
/**
* @brief Sets the total size of the data in bytes.
*
* It's important to ensure this matches the actual data size in memory to prevent overflows
* or incorrect computations downstream.
*
* @param byteSize The total byte size to set.
*/
void Metadata::setByteSize(std::size_t byteSize) noexcept {
byteSize_ = byteSize;
}
/**
* @brief Gets the type of the data.
*
* The type (e.g., "float", "double") is critical for casting raw data or interfacing with libraries
* that require specific types.
*
* @return A string representing the data type.
*/
const std::string& Metadata::getDataType() const noexcept {
return dataType_;
}
/**
* @brief Sets the type of the data.
*
* When setting the data type, ensure it aligns with the underlying data representation.
* Mismatched types can lead to undefined behavior in numerical calculations.
*
* @param dataType A string representing the data type.
*/
void Metadata::setDataType(const std::string& dataType) {
dataType_ = dataType;
}
/**
* @brief Gets the dimensions of the data.
*
* Dimensions define the shape of the data (e.g., 2D arrays, 3D matrices). This is essential
* for ensuring that operations (e.g., matrix multiplication) are performed correctly.
*
* @return A vector representing the size of each dimension.
*/
const std::vector<std::size_t>& Metadata::getDimensions() const noexcept {
return dimensions_;
}
/**
* @brief Sets the dimensions of the data.
*
* When modifying dimensions, verify that they are consistent with the actual data representation.
*
* @param dimensions A vector representing the size of each dimension.
*/
void Metadata::setDimensions(const std::vector<std::size_t>& dimensions) {
dimensions_ = dimensions;
}
/**
* @brief Checks if debugging information is enabled.
*
* Debugging flags can be useful for tracking performance metrics or error provenance.
*
* @return True if debugging is enabled, false otherwise.
*/
bool Metadata::isDebugEnabled() const noexcept {
return debugFlag_;
}
/**
* @brief Sets the debugging flag.
*
* Enabling debugging can introduce performance overhead but provides valuable insights
* during development or testing.
*
* @param debugFlag Whether debugging is enabled.
*/
void Metadata::setDebugEnabled(bool debugFlag) noexcept {
debugFlag_ = debugFlag;
}
/**
* @brief Prints the metadata information for debugging purposes.
*
* This function provides a human-readable summary of the metadata. Useful for quick checks
* or logging during debugging sessions.
*
* @param os The output stream to print to.
* @param metadata The Metadata object to print.
* @return A reference to the output stream.
*/
std::ostream& operator<<(std::ostream& os, const Metadata& metadata) {
os << "Metadata Information:\n";
os << " Byte Size: " << metadata.byteSize_ << " bytes\n";
os << " Data Type: " << metadata.dataType_ << "\n";
os << " Dimensions: [";
for (size_t i = 0; i < metadata.dimensions_.size(); ++i) {
os << metadata.dimensions_[i];
if (i < metadata.dimensions_.size() - 1) os << ", ";
}
os << "]\n";
os << " Debug Enabled: " << (metadata.debugFlag_ ? "Yes" : "No") << "\n";
return os;
}

257
src/dobj/public/DObject.h
View File

@@ -1,257 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 04, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef DOBJECT_H
#define DOBJECT_H
#include <variant>
#include <memory>
#include <vector>
#include <string>
#include <mutex>
#include <map>
#include <functional>
/**
* @file DObject.h
* @brief Defines the DObject class, a universal data container for the project.
*
* The DObject class encapsulates arbitrary data and its associated metadata,
* providing a consistent interface for public-facing functions. It includes
* support for dynamically registered plugins.
*/
/**
* @class DObject
* @brief A universal data container class.
*
* The DObject class is designed to store arbitrary data alongside descriptive metadata.
* It supports plugin registration to allow extensible functionality.
*
* The general purpose of this is to simplify memory managment, function interfaces,
* and interoperability with other languages by wrapping all of that up inside of
* one location (DObject). There are still a limited number of types that DOBject can
* represent these include
* - bool
* - short int
* - int
* - long int
* - float
* - double
* - long double
* - std::string
* - std::monostate
* - std::vector<int>
* - std::vector<float>
* - std::vector<double>
* - std::vector<std::string>
* - std::vector<std::vector<int>>
* - std::vector<std::vector<float>>
* - std::vector<std::vector<double>>
* - std::vector<std::vector<std::vector<int>>>
* - std::vector<std::vector<std::vector<float>>>
* - std::vector<std::vector<std::vector<double>>
*
*
* If more types are needed, they can be added to the DataType variant.
*
* In general the usage of this might look like
* @code
* DObject obj;
* std::vector<int> data = {1, 2, 3, 4, 5};
* obj.setData(data);
* std::cout << "Data is " << obj << std::end;
* someFunction(&obj);
* @endcode
*
* All memory managment should be taken care of
* as all datatypes accepted by the varient are either primatives or
* include their own built in memory managment (i.e. this means that
* wherever possible vectors should be used in place of raw arrays).
*/
class DObject {
public:
/**
* @brief Supported data types for the DObject.
*
* This type alias uses `std::variant` to store different types of data,
* ensuring type safety and flexibility.
*/
using DataType = std::variant<
bool, short int, int, long int, float, double,
long double, std::string, std::monostate, std::vector<int>,
std::vector<float>, std::vector<double>, std::vector<std::string>,
std::vector<std::vector<int>>, std::vector<std::vector<float>>,
std::vector<std::vector<double>>,
std::vector<std::vector<std::vector<int>>>,
std::vector<std::vector<std::vector<float>>>,
std::vector<std::vector<std::vector<double>>>
>;
// Mapping of data types to their string representations
std::map<int, std::string> dataTypeMap = {
{0, "bool"}, {1, "short int"}, {2, "int"}, {3, "long int"}, {4, "float"},
{5, "double"}, {6, "long double"}, {7, "string"}, {8, "std::monostate"},
{9, "vector<int>"}, {10, "vector<float>"}, {11, "vector<double>"},
{12, "vector<string>"}, {13, "vector<vector<int>"},
{14, "vector<vector<float>"}, {15, "vector<vector<double>"},
{16, "vector<vector<vector<int>>"}, {17, "vector<vector<vector<float>>"},
{18, "vector<vector<vector<double>>"}
};
/**
* @brief Placeholder type for plugins.
*
* In the future, this will be replaced with a concrete interface.
*/
using Plugin = std::function<void(DObject&)>;
/**
* @brief Default constructor.
*
* Creates an empty DObject.
*/
DObject();
/**
* @brief Constructor to initialize a DObject with data.
*
* @param data The data to be stored in the DObject.
*/
DObject(const DataType& data);
/**
* @brief Retrieves the data stored in the DObject.
*
* Use the appropriate type (matching the stored data) with `std::get<T>()`.
*
* @return A constant reference to the stored data.
*/
const DataType& getData() const noexcept;
/**
* @brief Retrieves the data stored in the DObject as a typed object so that std::get<T>() is not needed.
*
* @return Data as type T
*/
template <typename T>
T getDataAs() const {
if (std::holds_alternative<T>(data_)) {
return std::get<T>(data_);
} else {
throw std::bad_variant_access();
}
}
/**
* @brief Sets the data for the DObject.
*
* Updates the stored data and optionally updates metadata.
*
* @param data The new data to store in the DObject.
*/
void setData(const DataType& data);
/**
* @brief Enables or disables debugging and tracing for the DObject.
*
* When debugging is enabled, the DObject tracks creation and modification history.
*
* @param enableDebug True to enable debugging, false to disable it.
*/
void setDebugging(bool enableDebug);
/**
* @brief Checks if debugging is enabled for the DObject.
*
* @return True if debugging is enabled, false otherwise.
*/
bool isDebuggingEnabled() const noexcept;
/**
* @breif Set error code tracked by the DOBject
*
* @return The previous error code
*/
int setErrorCode(int code) noexcept;
/**
* @brief Get the error code tracked by the DObject
*
* @return The error code
*/
int getErrorCode() const noexcept;
/**
* @brief Registers a plugin with the DObject.
*
* Plugins are stored in a registry and can add custom functionality to the DObject.
*
* @param id A unique identifier for the plugin.
* @param plugin The plugin function to register.
*/
void registerPlugin(const std::string& id, Plugin plugin);
/**
* @brief Unregisters a plugin by its identifier.
*
* Removes the plugin from the registry if it exists.
*
* @param id The unique identifier of the plugin to unregister.
*/
void unregisterPlugin(const std::string& id);
/**
* @brief Executes a plugin by its identifier.
*
* Invokes the registered plugin function. If the plugin is not found, no action is taken.
*
* @param id The unique identifier of the plugin to execute.
*/
void runPlugin(const std::string& id);
/**
* @brief Executes all registered plugins in the registry.
*
* Iterates through all plugins and invokes them on the current DObject.
*/
void runAllPlugins();
/**
* @brief Provides a human-readable summary of the DObject.
*
* Useful for quick inspection or logging during debugging sessions.
*
* @param os The output stream to write to.
* @param obj The DObject to summarize.
* @return A reference to the output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const DObject& obj);
private:
DataType data_; ///< The main data stored in the DObject.
bool debugEnabled_ = false; ///< Indicates whether debugging is enabled.
int errorCode_ = 0; ///< Error code tracked by the DObject.
std::map<std::string, Plugin> plugins_; ///< Registry for dynamically registered plugins.
};
#endif // DOBJECT_H

67
src/dobj/public/LockableDObject.h
View File

@@ -1,67 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: January 19, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef LOCKABLE_DOBJECT_H
#define LOCKABLE_DOBJECT_H
#include "DObject.h"
#include <mutex>
/**
* @file LockableDObject.h
* @brief A lightweight wrapper for DObject that adds locking capabilities.
*
* This class allows safe concurrent access to a DObject by providing
* locking and unlocking methods.
*/
/**
* @class LockableDObject
* @brief Wrapper for DObject with thread-safe access.
*/
class LockableDObject {
public:
/**
* @brief Default constructor.
*/
LockableDObject() = default;
/**
* @brief Access the underlying DObject.
* @return A reference to the wrapped DObject.
*/
DObject& get();
/**
* @brief Locks the mutex to ensure thread-safe access.
*/
void lock();
/**
* @brief Unlocks the mutex after thread-safe access.
*/
void unlock();
private:
DObject object_; ///< The underlying DObject instance.
std::mutex mutex_; ///< Mutex for thread-safe access.
};
#endif // LOCKABLE_DOBJECT_H

132
src/dobj/public/Metadata.h
View File

@@ -1,132 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: January 22, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef METADATA_H
#define METADATA_H
#if defined(__APPLE__) || defined(__linux__)
#define EXPORT_SYMBOL __attribute__((visibility("default")))
#else
#define EXPORT_SYMBOL
#endif
#include <string>
#include <vector>
#include <cstddef>
#include <iostream>
/**
* @file Metadata.h
* @brief Public interface for the Metadata class used in the dobj module.
*
* The Metadata class provides descriptive information about the data encapsulated
* within a dobj, including size, type, dimensions, and debugging flags.
*/
/**
* @class Metadata
* @brief Represents metadata information for data objects in the dobj module.
*
* The Metadata class encapsulates details such as data size, type, dimensions,
* and optional debugging flags. It is designed to provide descriptive attributes
* in a lightweight and efficient manner.
*/
class EXPORT_SYMBOL Metadata {
public:
/**
* @brief Default constructor for Metadata.
*/
Metadata() = default;
/**
* @brief Constructor to initialize Metadata with specific attributes.
*
* @param byteSize The total size of the data in bytes.
* @param dataType The type of the data (e.g., "float", "double").
* @param dimensions A vector representing the size of each dimension (e.g., {3, 4} for a 3x4 matrix).
* @param debugFlag Whether debugging information is enabled for this Metadata instance.
*/
Metadata(std::size_t byteSize, std::string dataType, std::vector<std::size_t> dimensions, bool debugFlag = false);
/**
* @brief Gets the total size of the data in bytes.
* @return The total byte size of the data.
*/
[[nodiscard]] std::size_t getByteSize() const noexcept;
/**
* @brief Sets the total size of the data in bytes.
* @param byteSize The total byte size to set.
*/
void setByteSize(std::size_t byteSize) noexcept;
/**
* @brief Gets the type of the data.
* @return A string representing the data type.
*/
[[nodiscard]] const std::string& getDataType() const noexcept;
/**
* @brief Sets the type of the data.
* @param dataType A string representing the data type.
*/
void setDataType(const std::string& dataType);
/**
* @brief Gets the dimensions of the data.
* @return A vector representing the size of each dimension.
*/
[[nodiscard]] const std::vector<std::size_t>& getDimensions() const noexcept;
/**
* @brief Sets the dimensions of the data.
* @param dimensions A vector representing the size of each dimension.
*/
void setDimensions(const std::vector<std::size_t>& dimensions);
/**
* @brief Checks if debugging information is enabled.
* @return True if debugging is enabled, false otherwise.
*/
[[nodiscard]] bool isDebugEnabled() const noexcept;
/**
* @brief Sets the debugging flag.
* @param debugFlag Whether debugging is enabled.
*/
void setDebugEnabled(bool debugFlag) noexcept;
/**
* @brief Prints the metadata information for debugging purposes.
*
* @param os The output stream to print to.
* @param metadata The Metadata object to print.
* @return A reference to the output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const Metadata& metadata);
private:
int byteSize_ = 0; ///< Total size of the data in bytes.
std::string dataType_; ///< Type of the data (e.g., "float", "double").
std::vector<std::size_t> dimensions_; ///< Dimensions of the data (e.g., {3, 4} for a 3x4 matrix).
bool debugFlag_ = false; ///< Indicates whether debugging is enabled.
};
#endif // METADATA_H

75
src/eos/private/EOSio.cpp
View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <string>
#include <utility>
@@ -5,33 +25,36 @@
#include "helm.h"
#include "debug.h"
EOSio::EOSio(std::string filename) : m_filename(std::move(filename)) {
load();
}
#include <string>
std::string EOSio::getFormat() const {
return m_format;
}
EOSTable& EOSio::getTable() {
return m_table;
}
void EOSio::load() {
// Load the EOS table from the file
// For now, just set the format to HELM
m_format = "helm";
if (m_format == "helm") {
loadHelm();
namespace serif::eos {
EOSio::EOSio(const std::string filename) : m_filename(filename) {
load();
}
}
void EOSio::loadHelm() {
// Load the HELM table from the file
auto helmTabPtr = helmholtz::read_helm_table(m_filename);
m_table = std::move(helmTabPtr);
m_loaded = true;
}
std::string EOSio::getFormat() const {
return m_format;
}
EOSTable& EOSio::getTable() {
return m_table;
}
void EOSio::load() {
// Load the EOS table from the file
// For now, just set the format to HELM
m_format = "helm";
if (m_format == "helm") {
loadHelm();
}
}
void EOSio::loadHelm() {
// Load the HELM table from the file
auto helmTabptr = serif::eos::helmholtz::read_helm_table(m_filename);
m_table = std::move(helmTabptr);
m_loaded = true;
}
}

108
src/eos/private/helm.cpp
View File

@@ -1,8 +1,8 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 17, 2025
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -41,56 +41,60 @@
#include "config.h"
#include "quill/LogMacros.h"
namespace serif::constant {
class Constants;
}
using namespace std;
double** heap_allocate_contiguous_2D_memory(const int rows, const int cols) {
// interpolating polynomila function definitions
namespace serif::eos::helmholtz {
double** heap_allocate_contiguous_2D_memory(const int rows, const int cols) {
if (rows <= 0 || cols <= 0) {
throw std::invalid_argument("Rows and columns must be positive integers.");
throw std::invalid_argument("Rows and columns must be positive integers.");
}
auto array = new double*[rows];
const auto array = new double*[rows];
array[0] = new double[rows * cols];
for (int i = 1; i < rows; i++) {
array[i] = array[0] + i * cols;
array[i] = array[0] + i * cols;
}
return array;
}
}
void heap_deallocate_contiguous_2D_memory(double **array) {
void heap_deallocate_contiguous_2D_memory(double **array) {
delete[] array[0];
delete[] array;
}
// interpolating polynomila function definitions
namespace helmholtz {
double psi0(double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }
}
double psi0(const double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }
double dpsi0(double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }
double dpsi0(const double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }
double ddpsi0(double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }
double ddpsi0(const double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }
double psi1(double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }
double psi1(const double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }
double dpsi1(double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }
double dpsi1(const double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }
double ddpsi1(double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }
double ddpsi1(const double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }
double psi2(double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }
double psi2(const double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }
double dpsi2(double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }
double dpsi2(const double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }
double ddpsi2(double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }
double ddpsi2(const double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }
double xpsi0(double z) { return z * z * (2.0 * z - 3.0) + 1.0; }
double xpsi0(const double z) { return z * z * (2.0 * z - 3.0) + 1.0; }
double xdpsi0(double z) { return z * (6.0*z - 6.0); }
double xdpsi0(const double z) { return z * (6.0*z - 6.0); }
double xpsi1(double z) { return z * ( z * (z - 2.0) + 1.0); }
double xpsi1(const double z) { return z * ( z * (z - 2.0) + 1.0); }
double xdpsi1(double z) { return z * (3.0 * z - 4.0) + 1.0; }
double xdpsi1(const double z) { return z * (3.0 * z - 4.0) + 1.0; }
double h3(std::array<double, 36> fi, double w0t, double w1t, double w0mt, double w1mt,
double w0d, double w1d, double w0md, double w1md) {
double h3(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w0mt, const double w1mt,
const double w0d, const double w1d, const double w0md, const double w1md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
@@ -101,9 +105,9 @@ namespace helmholtz {
+ fi[14] * w1d * w1mt + fi[15] * w1md * w1mt;
}
double h5(std::array<double, 36> fi, double w0t, double w1t, double w2t,double w0mt,
double w1mt, double w2mt, double w0d, double w1d, double w2d,
double w0md, double w1md, double w2md) {
double h5(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w2t, const double w0mt,
const double w1mt, const double w2mt, const double w0d, const double w1d, const double w2d,
const double w0md, const double w1md, const double w2md) {
return fi[0] * w0d * w0t + fi[1] * w0md * w0t
+ fi[2] * w0d * w0mt + fi[3] * w0md * w0mt
+ fi[4] * w0d * w1t + fi[5] * w0md * w1t
@@ -125,15 +129,15 @@ namespace helmholtz {
}
// this function reads in the HELM table and stores in the above arrays
std::unique_ptr<HELMTable> read_helm_table(const std::string filename) {
Config& config = Config::getInstance();
std::string logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
Probe::LogManager& logManager = Probe::LogManager::getInstance();
std::unique_ptr<HELMTable> read_helm_table(const std::string &filename) {
serif::config::Config& config = serif::config::Config::getInstance();
auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);
// Make a unique pointer to the HELMTable
std::unique_ptr<HELMTable> table = std::make_unique<HELMTable>();
auto table = std::make_unique<serif::eos::helmholtz::HELMTable>();
string data;
int i, j;
@@ -236,13 +240,13 @@ namespace helmholtz {
ion, radiation, electron-positron and Coulomb interaction
and returns the calculated quantities in the input
***/
EOS get_helm_EOS(EOSInput &q, const HELMTable &table) {
Config& config = Config::getInstance();
std::string logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
Probe::LogManager& logManager = Probe::LogManager::getInstance();
serif::eos::helmholtz::EOS get_helm_EOS(serif::eos::helmholtz::EOSInput &q, const serif::eos::helmholtz::HELMTable &table) {
serif::config::Config& config = serif::config::Config::getInstance();
auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
Constants& constants = Constants::getInstance();
serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double pi = std::numbers::pi;
const double amu = constants.get("u").value;
const double h = constants.get("h").value;
@@ -252,13 +256,13 @@ namespace helmholtz {
const double kergavo = kerg*avo;
const double clight2 = clight * clight;
const double ssol = 5.6704e-5;
constexpr double ssol = 5.6704e-5;
const double asol = 4 * ssol / clight;
const double asoli3 = asol / 3;
const double sioncon = (2.0 * pi * amu * kerg)/(h*h);
const double qe = 4.8032042712e-10;
const double esqu = qe * qe;
const double third = 1./3.;
constexpr double qe = 4.8032042712e-10;
constexpr double esqu = qe * qe;
constexpr double third = 1./3.;
std::array<double, 36> fi;
double T, rho, s, free, abar, zbar, ytot1, ye, ptot;
double prad, srad, erad, etot, stot;
@@ -630,14 +634,14 @@ namespace helmholtz {
// coulomb
const double a1 = -0.898004;
const double b1 = 0.96786;
const double c1 = 0.220703;
const double d1 = -0.86097;
const double e1 = 2.5269;
const double a2 = 0.29561;
const double b2 = 1.9885;
const double c2 = 0.288675;
constexpr double a1 = -0.898004;
constexpr double b1 = 0.96786;
constexpr double c1 = 0.220703;
constexpr double d1 = -0.86097;
constexpr double e1 = 2.5269;
constexpr double a2 = 0.29561;
constexpr double b2 = 1.9885;
constexpr double c2 = 0.288675;
z = 4 * pi * third;
s = z * xni;
@@ -825,7 +829,7 @@ namespace helmholtz {
double csound = clight * sqrt(gamma1/z);
// package in q:
EOS eos;
serif::eos::helmholtz::EOS eos;
eos.ptot = ptot; eos.etot = etot; eos.stot = stot;
eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas;
eos.prad = prad; eos.erad = erad; eos.srad = srad;

135
src/eos/public/EOSio.h
View File

@@ -1,69 +1,92 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <string>
#include <variant>
#include <memory>
// EOS table format includes
#include "helm.h"
using EOSTable = std::variant<
std::unique_ptr<helmholtz::HELMTable>
>;
namespace serif::eos {
/**
* @class EOSio
* @brief Handles the input/output operations for EOS tables.
*
* The EOSio class is responsible for loading and managing EOS tables from files.
* It supports different formats, currently only HELM format.
*
* Example usage:
* @code
* EOSio eosIO("path/to/file");
* std::string format = eosIO.getFormat();
* EOSTable& table = eosIO.getTable();
* @endcode
*/
class EOSio {
private:
std::string m_filename; ///< The filename of the EOS table.
bool m_loaded = false; ///< Flag indicating if the table is loaded.
std::string m_format; ///< The format of the EOS table.
EOSTable m_table; ///< The EOS table data.
/**
* @brief Loads the EOS table from the file.
*/
void load();
// EOS table format includes
using EOSTable = std::variant<
std::unique_ptr<serif::eos::helmholtz::HELMTable>
>;
/**
* @brief Loads the HELM format EOS table.
* @class EosIO
* @brief Handles the input/output operations for EOS tables.
*
* The EosIO class is responsible for loading and managing EOS tables from files.
* It supports different formats, currently only HELM format.
*
* Example usage:
* @code
* EosIO eosIO("path/to/file");
* std::string format = eosIO.getFormat();
* EOSTable& table = eosIO.getTable();
* @endcode
*/
void loadHelm();
public:
/**
* @brief Constructs an EosIO object with the given filename.
* @param filename The filename of the EOS table.
*/
explicit EOSio(std::string filename);
/**
* @brief Default destructor.
*/
~EOSio() = default;
class EOSio {
private:
std::string m_filename; ///< The filename of the EOS table.
bool m_loaded = false; ///< Flag indicating if the table is loaded.
std::string m_format; ///< The format of the EOS table.
EOSTable m_table; ///< The EOS table data.
/**
* @brief Gets the format of the EOS table.
* @return The format of the EOS table as a string.
*/
std::string getFormat() const;
/**
* @brief Loads the EOS table from the file.
*/
void load();
/**
* @brief Gets the EOS table.
* @return A reference to the EOS table.
*/
EOSTable& getTable();
/**
* @brief Loads the HELM format EOS table.
*/
void loadHelm();
public:
/**
* @brief Constructs an EosIO object with the given filename.
* @param filename The filename of the EOS table.
*/
explicit EOSio(const std::string filename);
/**
* @brief Default destructor.
*/
~EOSio() = default;
/**
* @brief Gets the format of the EOS table.
* @return The format of the EOS table as a string.
*/
std::string getFormat() const;
/**
* @brief Gets the EOS table.
* @return A reference to the EOS table.
*/
EOSTable& getTable();
std::string getFilename() const { return m_filename; }
};
}
std::string getFilename() const { return m_filename; }
};

30
src/eos/public/helm.h
View File

@@ -2,7 +2,7 @@
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 06, 2025
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -21,10 +21,9 @@
// this file contains a C++ conversion of Frank Timmes' fortran code
// helmholtz.f90, which implements the Helmholtz Free Energy EOS described
// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025
#pragma once
#define IMAX 541
#define JMAX 201
#include <array>
#include <iostream>
@@ -35,21 +34,22 @@
#include "debug.h"
/**
* @brief 2D array template alias.
* @tparam T Type of the array elements.
* @tparam J Number of columns.
* @tparam I Number of rows.
*/
template <typename T, std::size_t J, std::size_t I>
using array2D = std::array<std::array<T, J>, I>;
namespace serif::eos::helmholtz {
constexpr int IMAX = 541;
constexpr int JMAX = 201;
/**
* @brief 2D array template alias.
* @tparam T Type of the array elements.
* @tparam J Number of columns.
* @tparam I Number of rows.
*/
template <typename T, std::size_t J, std::size_t I>
using array2D = std::array<std::array<T, J>, I>;
double** heap_allocate_contiguous_2D_memory(int rows, int cols);
double** heap_allocate_contiguous_2D_memory(int rows, int cols);
void heap_deallocate_contiguous_2D_memory(double **array);
void heap_deallocate_contiguous_2D_memory(double **array);
namespace helmholtz
{
const double tlo = 3.0, thi = 13.0;
const double dlo = -12.0, dhi = 15.0;
const double tstp = (thi - tlo) / (JMAX - 1), tstpi = 1 / tstp;

87
src/meshIO/private/meshIO.cpp
View File

@@ -26,60 +26,57 @@
#include "meshIO.h"
namespace serif::mesh {
MeshIO::MeshIO(const std::string &mesh_file, double scale_factor)
{
mesh_file_ = mesh_file;
std::ifstream mesh_stream(mesh_file);
if (!mesh_stream)
{
throw std::runtime_error("Mesh file not found: " + mesh_file);
}
m_mesh = std::make_unique<mfem::Mesh>(mesh_stream, 1, 10);
MeshIO::MeshIO(const std::string &mesh_file, double scale_factor)
{
mesh_file_ = mesh_file;
std::ifstream mesh_stream(mesh_file);
if (!mesh_stream)
{
throw std::runtime_error("Mesh file not found: " + mesh_file);
loaded_ = false;
}
else
{
mesh_ = mfem::Mesh(mesh_stream, 1, 10);
bool discont = false;
int space_dim = m_mesh->SpaceDimension();
int ordering = mfem::Ordering::byVDIM;
// m_mesh->SetCurvature(3, discont, space_dim, ordering);
m_mesh->Finalize(false, true);
loaded_ = true;
if (scale_factor != 1.0) { LinearRescale(scale_factor); }
}
}
MeshIO::~MeshIO()
{
}
MeshIO::~MeshIO() = default;
void MeshIO::LinearRescale(double scale_factor) {
if (!loaded_) {
throw std::runtime_error("Mesh not loaded before rescaling.");
void MeshIO::LinearRescale(const double scale_factor) const {
if (!loaded_) {
throw std::runtime_error("Mesh not loaded before rescaling.");
}
if (scale_factor <= 0.0) {
throw std::invalid_argument("Scale factor must be positive.");
}
// Ensure there are nodes even for linear order meshes
m_mesh->EnsureNodes();
const mfem::GridFunction *nodes = m_mesh->GetNodes();
double *node_data = nodes->GetData();
const int data_size = nodes->Size();
for (int i = 0; i < data_size; ++i) {
node_data[i] *= scale_factor;
}
// nodes->Update(); /updates the fespace
m_mesh->NodesUpdated();
}
if (scale_factor <= 0.0) {
throw std::invalid_argument("Scale factor must be positive.");
bool MeshIO::IsLoaded() const
{
return loaded_;
}
// Ensure there are nodes even for linear order meshes
mesh_.EnsureNodes();
const mfem::GridFunction *nodes = mesh_.GetNodes();
double *node_data = nodes->GetData(); // THIS IS KEY
int data_size = nodes->Size();
for (int i = 0; i < data_size; ++i) {
node_data[i] *= scale_factor;
mfem::Mesh &MeshIO::GetMesh() const noexcept {
return *m_mesh;
}
// nodes->Update(); /updates the fespace
mesh_.NodesUpdated();
}
bool MeshIO::IsLoaded() const
{
return loaded_;
}
mfem::Mesh& MeshIO::GetMesh()
{
return mesh_;
}

74
src/meshIO/public/meshIO.h
View File

@@ -18,51 +18,51 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef MESHIO_H
#define MESHIO_H
#pragma once
#include "mfem.hpp"
#include <string>
/**
* @brief Class for handling mesh input/output operations.
*/
class MeshIO
{
private:
bool loaded_; ///< Flag to indicate if the mesh is loaded
std::string mesh_file_; ///< Filename of the mesh file
mfem::Mesh mesh_; ///< The mesh object
public:
namespace serif::mesh {
/**
* @brief Constructor that initializes the MeshIO object with a mesh file.
* @param mesh_file The name of the mesh file.
* @brief Class for handling mesh input/output operations.
*/
MeshIO(const std::string &mesh_file, double scale_factor = 1.0);
class MeshIO
{
private:
bool loaded_; ///< Flag to indicate if the mesh is loaded
std::string mesh_file_; ///< Filename of the mesh file
std::unique_ptr<mfem::Mesh> m_mesh; ///< The mesh object
/**
* @brief Destructor for the MeshIO class.
*/
~MeshIO();
public:
/**
* @brief Constructor that initializes the MeshIO object with a mesh file.
* @param mesh_file The name of the mesh file.
* @param scale_factor
*/
explicit MeshIO(const std::string &mesh_file, double scale_factor = 1.0);
/**
* @brief Rescale the mesh by a linear factor.
* @param scale_factor The factor by which to scale the mesh.
*/
void LinearRescale(double scale_factor);
/**
* @brief Destructor for the MeshIO class.
*/
~MeshIO();
/**
* @brief Get the mesh object.
* @return Reference to the mesh object.
*/
mfem::Mesh& GetMesh();
/**
* @brief Rescale the mesh by a linear factor.
* @param scale_factor The factor by which to scale the mesh.
*/
void LinearRescale(double scale_factor) const;
/**
* @brief Check if the mesh is loaded.
* @return True if the mesh is loaded, false otherwise.
*/
bool IsLoaded() const;
};
/**
* @brief Get the mesh object.
* @return Pointer to the mesh object.
*/
[[nodiscard]] mfem::Mesh & GetMesh() const noexcept;
#endif // MESHIO_H
/**
* @brief Check if the mesh is loaded.
* @return True if the mesh is loaded, false otherwise.
*/
[[nodiscard]] bool IsLoaded() const;
};
}

5
src/meson.build
View File

@@ -7,10 +7,9 @@ subdir('types')
subdir('misc')
subdir('config')
subdir('probe')
subdir('dobj')
# Physically Informed Libraries
subdir('const')
subdir('constants')
subdir('composition')
# Asset Libraries
@@ -22,7 +21,7 @@ subdir('resource')
# Physics Libraries
subdir('network')
subdir('poly')
subdir('polytrope')
# Python Bindings
subdir('python')

21
src/misc/macros/debug.h
View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
/**
* @file debug.h
* @brief Defines a macro for triggering a breakpoint in different compilers and platforms.
@@ -10,6 +30,7 @@
* BREAKPOINT(); // Triggers a breakpoint in the debugger
* @endcode
*/
#pragma once
#ifdef __GNUC__ // GCC and Clang
/**

26
src/misc/macros/warning_control.h
View File

@@ -1,6 +1,24 @@
#ifndef WARNING_CONTROL_H
#define WARNING_CONTROL_H
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#if defined(__GNUC__) || defined(__clang__)
#define DEPRECATION_WARNING_OFF _Pragma("GCC diagnostic push") \
_Pragma("GCC diagnostic ignored \"-Wdeprecated-declarations\"")
@@ -12,5 +30,3 @@
#define DEPRECATION_WARNING_OFF
#define DEPRECATION_WARNING_ON
#endif
#endif // WARNING_CONTROL_H

32
src/network/private/approx8.cpp
View File

@@ -1,3 +1,23 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <cmath>
#include <stdexcept>
#include <array>
@@ -52,7 +72,7 @@ The coefficients to the fit are from reaclib.jinaweb.org .
*/
namespace nnApprox8{
namespace serif::network::approx8{
// using namespace std;
using namespace boost::numeric::odeint;
@@ -225,7 +245,7 @@ namespace nnApprox8{
// a Jacobian matrix for implicit solvers
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) {
Constants& constants = Constants::getInstance();
serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
// EOS
@@ -330,7 +350,7 @@ namespace nnApprox8{
}
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) {
Constants& constants = Constants::getInstance();
serif::constant::Constants& constants = serif::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
@@ -424,7 +444,7 @@ namespace nnApprox8{
dydt[Net::iener] = -enuc*avo*clight*clight;
}
nuclearNetwork::NetOut Approx8Network::evaluate(const nuclearNetwork::NetIn &netIn) {
NetOut Approx8Network::evaluate(const NetIn &netIn) {
m_y = convert_netIn(netIn);
m_tmax = netIn.tmax;
m_dt0 = netIn.dt0;
@@ -468,7 +488,7 @@ namespace nnApprox8{
m_y[i] /= ysum;
}
nuclearNetwork::NetOut netOut;
NetOut netOut;
std::vector<double> outComposition;
outComposition.reserve(Net::nvar);
@@ -486,7 +506,7 @@ namespace nnApprox8{
m_stiff = stiff;
}
vector_type Approx8Network::convert_netIn(const nuclearNetwork::NetIn &netIn) {
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
if (netIn.composition.size() != Net::niso) {
LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn");
throw std::runtime_error("Error: composition size mismatch in convert_netIn");

28
src/network/private/network.cpp
View File

@@ -1,14 +1,34 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "network.h"
#include "probe.h"
#include "quill/LogMacros.h"
namespace nuclearNetwork {
namespace serif::network {
Network::Network() :
m_config(Config::getInstance()),
m_logManager(Probe::LogManager::getInstance()),
m_config(serif::config::Config::getInstance()),
m_logManager(serif::probe::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")) {
}
nuclearNetwork::NetOut nuclearNetwork::Network::evaluate(const NetIn &netIn) {
NetOut Network::evaluate(const NetIn &netIn) {
// You can throw an exception here or log a warning if it should never be used.
LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented");
throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented");

33
src/network/public/approx8.h
View File

@@ -1,5 +1,24 @@
#ifndef APPROX8_H
#define APPROX8_H
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <array>
@@ -34,7 +53,7 @@ typedef boost::numeric::ublas::matrix< double > matrix_type;
*/
typedef std::array<double,7> vec7;
namespace nnApprox8{
namespace serif::network::approx8{
using namespace boost::numeric::odeint;
@@ -275,14 +294,14 @@ namespace nnApprox8{
* @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network.
*/
class Approx8Network : public nuclearNetwork::Network {
class Approx8Network : public Network {
public:
/**
* @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network.
* @return Output results from the network.
*/
virtual nuclearNetwork::NetOut evaluate(const nuclearNetwork::NetIn &netIn);
virtual NetOut evaluate(const NetIn &netIn);
/**
* @brief Sets whether the solver should use a stiff method.
@@ -306,9 +325,7 @@ namespace nnApprox8{
* @param netIn Input parameters for the network.
* @return Internal state vector.
*/
vector_type convert_netIn(const nuclearNetwork::NetIn &netIn);
vector_type convert_netIn(const NetIn &netIn);
};
} // namespace nnApprox8
#endif

31
src/network/public/network.h
View File

@@ -1,5 +1,24 @@
#ifndef NETWORK_H
#define NETWORK_H
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Emily Boudreaux, Aaron Dotter
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <vector>
@@ -7,7 +26,7 @@
#include "config.h"
#include "quill/Logger.h"
namespace nuclearNetwork {
namespace serif::network {
/**
* @struct NetIn
@@ -83,11 +102,9 @@ namespace nuclearNetwork {
virtual NetOut evaluate(const NetIn &netIn);
protected:
Config& m_config; ///< Configuration instance
Probe::LogManager& m_logManager; ///< Log manager instance
serif::config::Config& m_config; ///< Configuration instance
serif::probe::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance
};
} // namespace nuclearNetwork
#endif // NETWORK_H

43
src/opac/public/opac.h
View File

@@ -1,43 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: February 04, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef OPAC_H
#define OPAC_H
#include <string>
#include <map>
#include "DObject.h"
/**
* @breif initlaize the opacity module
* @param args: a hash map of all arguments needed to intialize opac
* @return error code in a DObject
**/
DObject initlaize_opac(const std::map<std::string, DObject>& args);
/**
* @breif opacity given a temperature, density, and composition
* @param args: a hash map of all arguments needed to calculate opac
* @return error code in a DObject
**/
DObject
#endif

0
src/poly/coeff/meson.build
View File

0
src/poly/solver/meson.build
View File

135
src/poly/utils/private/integrators.cpp
View File

@@ -1,135 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "mfem.hpp"
#include <cmath>
#include "integrators.h"
#include <string>
// static std::ofstream debugOut("gradient.csv", std::ios::trunc);
namespace polyMFEMUtils {
NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
m_polytropicIndex(n) {}
void NonlinearPowerIntegrator::AssembleElementVector(
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::Vector &elvect) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
int dof = el.GetDof();
elvect.SetSize(dof);
elvect = 0.0;
mfem::Vector shape(dof);
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
mfem::IntegrationPoint ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
const double u_safe = std::max(u_val, 0.0);
const double u_nl = std::pow(u_safe, m_polytropicIndex);
const double x2_u_nl = u_nl;
for (int i = 0; i < dof; i++){
elvect(i) += shape(i) * x2_u_nl * weight;
}
}
}
void NonlinearPowerIntegrator::AssembleElementGrad (
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::DenseMatrix &elmat) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elmat.SetSize(dof);
elmat = 0.0;
mfem::Vector shape(dof);
mfem::DenseMatrix dshape(dof, 3);
mfem::DenseMatrix invJ(3, 3);
mfem::Vector gradPhys(3);
mfem::Vector physCoord(3);
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
// Calculate the Jacobian
// TODO: Check for when theta ~ 0?
const double u_safe = std::max(u_val, 0.0);
const double d_u_nl = m_polytropicIndex * std::pow(u_safe, m_polytropicIndex - 1);
const double x2_d_u_nl = d_u_nl;
for (int i = 0; i < dof; i++) {
for (int j = 0; j < dof; j++) {
elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight;
}
}
// // --- Debug Code to write out gradient ---
// Trans.Transform(ip,physCoord);
// el.CalcDShape(ip, dshape);
//
// mfem::CalcInverse(Trans.Jacobian(), invJ);
//
// mfem::DenseMatrix dshapePhys;
// dshapePhys.SetSize(dof, physCoord.Size());
// mfem::Mult(dshape, invJ, dshapePhys);
//
// gradPhys = 0.0;
// for (int j = 0; j < dof; j++) {
// for (int d = 0; d < gradPhys.Size(); d++) {
// gradPhys(d) += elfun(j)*dshapePhys(j, d);
// }
// }
//
// debugOut
// << physCoord(0) << ", " << physCoord(1) << ", " << physCoord(2)
// << ", " << gradPhys(0) << ", " << gradPhys(1) << ", " << gradPhys(2) << '\n';
}
// debugOut.flush();
}
} // namespace polyMFEMUtils

336
src/poly/utils/private/operator.cpp
View File

@@ -1,336 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "operator.h"
#include "4DSTARTypes.h"
#include "mfem.hpp"
#include <memory>
static int s_PolyOperatorMultCalledCount = 0;
void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfem::DenseMatrix *mat) {
if (!mat) {
throw std::runtime_error("The operator is not a SparseMatrix.");
}
std::ofstream outfile(filename);
if (!outfile.is_open()) {
throw std::runtime_error("Failed to open file: " + filename);
}
int height = mat->Height();
int width = mat->Width();
// Set precision for floating-point output
outfile << std::fixed << std::setprecision(precision);
for (int i = 0; i < width; i++) {
outfile << i;
if (i < width - 1) {
outfile << ",";
}
else {
outfile << "\n";
}
}
// Iterate through rows
for (int i = 0; i < height; ++i) {
for (int j = 0; j < width; ++j) {
outfile << mat->Elem(i, j);
if (j < width - 1) {
outfile << ",";
}
}
outfile << std::endl;
}
outfile.close();
}
/**
* @brief Writes the dense representation of an MFEM Operator (if it's a SparseMatrix) to a CSV file.
*
* @param op The MFEM Operator to write.
* @param filename The name of the output CSV file.
* @param precision Number of decimal places for floating-point values.
*/
void writeOperatorToCSV(const mfem::Operator &op,
const std::string &filename,
const int precision = 6) // Add precision argument
{
// Attempt to cast the Operator to a SparseMatrix
const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&op);
if (!sparse_mat) {
throw std::runtime_error("The operator is not a SparseMatrix.");
}
const mfem::DenseMatrix *mat = sparse_mat->ToDenseMatrix();
writeDenseMatrixToCSV(filename, precision, mat);
}
void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::SparseMatrix>& invMat, const std::string& name="matrix") {
// PERF: This likely can be made much more efficient and will probably be called in tight loops, a good
// PERF: place for some easy optimization might be here.
// Confirm that mat is a square matrix
MFEM_ASSERT(mat.Height() == mat.Width(), "Matrix " + name + " is not square, cannot invert.");
mfem::Vector diag;
mat.GetDiag(diag);
// Invert the diagonal
for (int i = 0; i < diag.Size(); i++) {
MFEM_ASSERT(diag(i) != 0, "Diagonal element (" + std::to_string(i) +") in " + name + " is zero, cannot invert.");
diag(i) = 1.0 / diag(i);
}
// If the matrix is already inverted, just set the diagonal to avoid reallocation
if (invMat != nullptr) {
MFEM_ASSERT(invMat->Height() == invMat->Width(), "invMat (result matrix) is not square, cannot invert " + name + " into it.");
MFEM_ASSERT(invMat->Height() == mat.Height(), "Incompatible matrix sizes for inversion of " + name + ", expected " + std::to_string(mat.Height()) + " but got " + std::to_string(invMat->Height()));
for (int i = 0; i < diag.Size(); i++) {
MFEM_ASSERT(diag(i) != 0, "Diagonal element (" + std::to_string(i) +") in " + name + " is zero, resulting matrix would be singular.");
invMat->Elem(i, i) = diag(i);
}
} else { // The matrix has not been allocated yet so that needs to be done. Sparse Matrix has a constructor that can build from the diagonals
invMat = std::make_unique<mfem::SparseMatrix>(diag);
}
}
PolytropeOperator::PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> M,
std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets) :
mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector
m_blockOffsets(blockOffsets),
m_jacobian(nullptr) {
m_M = std::move(M);
m_Q = std::move(Q);
m_D = std::move(D);
m_f = std::move(f);
}
void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
if (m_isFinalized) {
return; // do nothing if already finalized
}
m_negM_op = std::make_unique<mfem::ScaledOperator>(m_M.get(), -1.0);
m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Q.get(), -1.0);
// Set up the constant parts of the jacobian now
m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant)
m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant)
m_jacobian->SetBlock(1, 1, m_D.get()); //<- D (constant)
// Build the initial preconditioner based on some initial guess
const auto &grad = m_f->GetGradient(initTheta);
updatePreconditioner(grad);
m_isFinalized = true;
}
const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const {
if (m_jacobian == nullptr) {
MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call.");
}
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetJacobianOperator().");
}
return *m_jacobian;
}
mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const {
if (m_schurPreconditioner == nullptr) {
MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call.");
}
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner().");
}
return *m_schurPreconditioner;
}
void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Mult called before finalize");
}
// -- Create BlockVector views for input x and output y
mfem::BlockVector x_block(const_cast<mfem::Vector&>(x), m_blockOffsets);
mfem::BlockVector y_block(y, m_blockOffsets);
// -- Get Vector views for individual blocks
const mfem::Vector &x_theta = x_block.GetBlock(0);
const mfem::Vector &x_phi = x_block.GetBlock(1);
mfem::Vector &y_R0 = y_block.GetBlock(0); // Residual Block 0 (theta)
mfem::Vector &y_R1 = y_block.GetBlock(1); // Residual Block 1 (phi)
int theta_size = m_blockOffsets[1] - m_blockOffsets[0];
int phi_size = m_blockOffsets[2] - m_blockOffsets[1];
mfem::Vector f_term(theta_size);
mfem::Vector Mphi_term(theta_size);
mfem::Vector Dphi_term(phi_size);
mfem::Vector Qtheta_term(phi_size);
// Calculate R0 and R1 terms
// R0 = f(θ) - Mɸ
// R1 = Dɸ - Qθ
MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult");
m_f->Mult(x_theta, f_term); // fixme: this may be the wrong way to assemble m_f?
m_M->Mult(x_phi, Mphi_term);
m_D->Mult(x_phi, Dphi_term);
m_Q->Mult(x_theta, Qtheta_term);
subtract(f_term, Mphi_term, y_R0);
subtract(Dphi_term, Qtheta_term, y_R1);
// -- Apply essential boundary conditions --
for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
if (int idx = m_theta_ess_tdofs.first[i]; idx >= 0 && idx < y_R0.Size()) {
const double &targetValue = m_theta_ess_tdofs.second[i];
// y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form
}
}
// TODO: Are the residuals for φ being calculated correctly?
std::ofstream outfileθ("PolyOperatorMultCallTheta_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
outfileθ << "dof,R_θ\n";
for (int i = 0; i < y_R0.Size(); i++) {
outfileθ << i << "," << y_R0(i) << "\n";
}
outfileθ.close();
std::ofstream outfileφ("PolyOperatorMultCallPhi_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
outfileφ << "dof,R_ɸ\n";
for (int i = 0; i < y_R1.Size(); i++) {
outfileφ << i << "," << y_R1(i) << "\n";
}
outfileφ.close();
++s_PolyOperatorMultCalledCount;
}
void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const {
if (const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&grad); sparse_mat != nullptr) {
approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian");
} else {
MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not dynamically castable to a sparse matrix");
}
}
void PolytropeOperator::updateInverseSchurCompliment() const {
// TODO Add a flag in to make sure this tracks in parallel (i.e. every time the non linear jacobian inverse is updated set the flag to true and then check if the flag is true here and if so do work (if not throw error). then at the end of this function set it to false.
if (m_invNonlinearJacobian == nullptr) {
MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
}
mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D
// Calculate S = D - Q df^{-1} M
mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(m_Q->SpMat(), *m_invNonlinearJacobian); // Q * df^{-1}
const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, m_M->SpMat()); // Q * df^{-1} * M
// PERF: Could potentially add some caching in here to only update the preconditioner when some condition has been met
schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
approxJacobiInvert(*schurCompliment, m_invSchurCompliment, "Schur Compliment");
if (m_schurPreconditioner == nullptr) {
m_schurPreconditioner = std::make_unique<mfem::BlockDiagonalPreconditioner>(m_blockOffsets);
}
// ⎡ḟ(θ)^-1 0⎤
// ⎣0 S^-1⎦
m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get());
m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get());
}
void PolytropeOperator::updatePreconditioner(const mfem::Operator &grad) const {
updateInverseNonlinearJacobian(grad);
updateInverseSchurCompliment();
}
mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::GetGradient called before finalize");
}
// --- Get the gradient of f ---
mfem::BlockVector x_block(const_cast<mfem::Vector&>(x), m_blockOffsets);
const mfem::Vector& x_theta = x_block.GetBlock(0);
auto &grad = m_f->GetGradient(x_theta);
updatePreconditioner(grad);
m_jacobian->SetBlock(0, 0, &grad);
// The other blocks are set up in finalize since they are constant. Only J00 depends on the current state of theta
return *m_jacobian;
}
void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) {
m_isFinalized = false;
m_theta_ess_tdofs = theta_ess_tdofs;
m_phi_ess_tdofs = phi_ess_tdofs;
if (m_f) {
m_f->SetEssentialTrueDofs(theta_ess_tdofs.first);
} else {
MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
}
if (m_M) {
m_M->EliminateTestEssentialBC(theta_ess_tdofs.first);
m_M->EliminateTrialEssentialBC(phi_ess_tdofs.first);
} else {
MFEM_ABORT("m_M is null in PolytropeOperator::SetEssentialTrueDofs");
}
if (m_Q) {
m_Q->EliminateTrialEssentialBC(theta_ess_tdofs.first);
m_Q->EliminateTestEssentialBC(phi_ess_tdofs.first);
} else {
MFEM_ABORT("m_Q is null in PolytropeOperator::SetEssentialTrueDofs");
}
if (m_D) {
m_D->EliminateEssentialBC(phi_ess_tdofs.first);
} else {
MFEM_ABORT("m_D is null in PolytropeOperator::SetEssentialTrueDofs");
}
}
void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) {
SetEssentialTrueDofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second);
}
SSE::MFEMArrayPairSet PolytropeOperator::GetEssentialTrueDofs() const {
return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs);
}

98
src/poly/utils/public/operator.h
View File

@@ -1,98 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include "mfem.hpp"
#include "4DSTARTypes.h"
#include <memory>
#include "probe.h"
class PolytropeOperator final : public mfem::Operator {
public:
PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> M,
std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets);
~PolytropeOperator() override = default;
void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
mfem::Operator& GetGradient(const mfem::Vector &x) const override;
void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs);
void SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set);
SSE::MFEMArrayPairSet GetEssentialTrueDofs() const;
bool isFinalized() const { return m_isFinalized; }
void finalize(const mfem::Vector &initTheta);
const mfem::Array<int>& GetBlockOffsets() const { return m_blockOffsets; }
const mfem::BlockOperator &GetJacobianOperator() const;
mfem::BlockDiagonalPreconditioner &GetPreconditioner() const;
private:
Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
std::unique_ptr<mfem::MixedBilinearForm> m_M;
std::unique_ptr<mfem::MixedBilinearForm> m_Q;
std::unique_ptr<mfem::BilinearForm> m_D;
std::unique_ptr<mfem::NonlinearForm> m_f;
const mfem::Array<int> m_blockOffsets;
SSE::MFEMArrayPair m_theta_ess_tdofs;
SSE::MFEMArrayPair m_phi_ess_tdofs;
std::unique_ptr<mfem::ScaledOperator> m_negM_op;
std::unique_ptr<mfem::ScaledOperator> m_negQ_op;
mutable std::unique_ptr<mfem::BlockOperator> m_jacobian;
mutable std::unique_ptr<mfem::SparseMatrix> m_invSchurCompliment;
mutable std::unique_ptr<mfem::SparseMatrix> m_invNonlinearJacobian;
/*
* The schur preconditioner has the form
*
* ⎡ḟ(θ)^-1 0 ⎤
* ⎣ 0 S^-1 ⎦
*
* Where S is the Schur compliment of the system
*
*/
mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner;
bool m_isFinalized = false;
private:
void updateInverseNonlinearJacobian(const mfem::Operator &grad) const;
void updateInverseSchurCompliment() const;
void updatePreconditioner(const mfem::Operator &grad) const;
};

23
src/polytrope/coeff/meson.build Normal file
View File

@@ -0,0 +1,23 @@
polyCoeff_sources = files(
'private/polyCoeff.cpp'
)
polyCoeff_headers = files(
'public/polyCoeff.h'
)
libPolyCoeff = shared_library('polyCoeff',
polyCoeff_sources,
include_directories : include_directories('./public'),
cpp_args: ['-fvisibility=default'],
dependencies: [mfem_dep],
install: true
)
polycoeff_dep = declare_dependency(
include_directories : include_directories('./public'),
link_with : libPolyCoeff,
sources : polyCoeff_sources,
dependencies : [mfem_dep]
)

64
src/polytrope/coeff/private/polyCoeff.cpp Normal file
View File

@@ -0,0 +1,64 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 19, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "mfem.hpp"
#include "polyCoeff.h"
namespace serif {
namespace polytrope {
namespace polycoeff{
double nonlinearSourceCoeff(const mfem::Vector &x)
{
return 1;
// double r = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
// return std::pow(r, 2);
}
void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v)
{
v.SetSize(3);
v = -1;
// double r = x(0)*x(0) + x(1)*x(1) + x(2)*x(2);
// v = -std::pow(r, 2);
}
//PERF: One area of future optimization might be turning these into matrix operations since they are
//PERF: fundamentally just the dot product of [a b c ...] * [n^0 n^1 n^2 ...]
double x1(const double n)
{
double r = 0;
for (int i = 0; i < x1InterpCoeff::numTerms; i++) {
r += x1InterpCoeff::coeff[i] * std::pow(n, x1InterpCoeff::power[i]);
}
return r;
}
double thetaSurfaceFlux(const double n) {
double surfaceFlux = 0;
for (int i = 0; i < dThetaInterpCoeff::numTerms; i++) {
surfaceFlux += dThetaInterpCoeff::coeff[i] * std::pow(n, dThetaInterpCoeff::power[i]);
}
return surfaceFlux;
}
} // namespace polycoeff
} // namespace polytrope
} // namespace serif

70
src/polytrope/coeff/public/polyCoeff.h Normal file
View File

@@ -0,0 +1,70 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include "mfem.hpp"
#include <cmath>
namespace serif {
namespace polytrope {
namespace polycoeff
{
/**
* @brief Computes the xi coefficient function.
*
* @param x Input vector.
* @return double The computed xi coefficient.
*/
double nonlinearSourceCoeff(const mfem::Vector &x);
/**
* @brief Computes the vector xi coefficient function.
*
* @param x Input vector.
* @param v Output vector to store the computed xi coefficient.
*/
void diffusionCoeff(const mfem::Vector &x, mfem::Vector &v);
double x1(const double n);
double thetaSurfaceFlux(const double n);
/**
* @brief Coefficients for the interpolations of the surface location of a polytrope
* @param numTerms Number of terms in the polynomial interpolator
* @param power Array of the powers of the polynomial interpolator
* @param coeff Array of the coefficients of the polynomial interpolator
*/
struct x1InterpCoeff {
constexpr static int numTerms = 51;
constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
constexpr static double coeff[51] = {2.449260049071003e+00, 4.340987403815444e-01, 2.592462347475860e+00, -2.283794512491552e+01, 1.135535208066129e+02, -3.460093673587010e+02, 6.948139716983651e+02, -9.564368645158952e+02, 9.184546798891624e+02, -6.130308987803503e+02, 2.747228735318193e+02, -7.416795903196430e+01, 7.318638538580859e+00, 1.749441306797260e+00, -4.240148582456829e-01, -5.818809544982156e-02, 1.514877217199105e-02, 3.228634707578998e-03, -2.862524323980516e-04, -1.622486968261819e-04, -1.253644717104076e-05, 4.334945141292894e-06, 1.296452565229763e-06, 8.634802209209870e-08, -3.337511676486084e-08, -1.094796628367775e-08, -1.228178760540410e-09, 1.416744125622751e-10, 8.513777351265677e-11, 1.624582561364811e-11, 8.377207519041114e-13, -4.363812865112836e-13, -1.535862757816461e-13, -2.485085045037669e-14, -6.566281276491033e-16, 8.405047965853478e-16, 2.673804441025638e-16, 4.176337890285142e-17, 8.150073570140493e-19, -1.531673805257016e-18, -4.746996933716653e-19, -6.976825828390195e-20, 8.807513368604331e-22, 3.307739310180278e-21, 8.593260940093030e-22, 7.385969061093440e-23, -2.211130577977291e-23, -8.557291048388455e-24, -4.169359901215994e-25, 4.609379358875657e-25, -3.590870335035984e-26};
};
struct dThetaInterpCoeff {
constexpr static int numTerms = 51;
constexpr static int power[51] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50};
constexpr static double coeff[51] = {-8.164634353536671e-01, 8.651804355981230e-01, -6.391939109867470e-01, 4.882118355278944e-01, -3.898475553686603e-01, 3.102560043891329e-01, -2.365565113144737e-01, 1.639832755778444e-01, -9.645271947133248e-02, 4.450140636696046e-02, -1.460885194751961e-02, 2.816293786806504e-03, -8.178583546493117e-05, -8.584494556484958e-05, 8.476252127593868e-06, 2.923593421315422e-06, -2.206768214995963e-07, -1.203227957690371e-07, -3.381181730985542e-09, 4.022824706790907e-09, 7.041049107708875e-10, -3.562885681170365e-11, -3.281525407784209e-11, -5.031807464141896e-12, 2.034136401832885e-13, 2.361284283178230e-13, 4.602774507763180e-14, 1.809170850970874e-15, -1.333813332262995e-15, -4.045891156434286e-16, -5.197949512114809e-17, 2.222220713310119e-18, 2.625223897130583e-18, 6.226001466529447e-19, 6.571419077089260e-20, -6.672159423054950e-21, -4.476242224056620e-21, -9.790792477821165e-22, -9.222211318122281e-23, 1.427942034536028e-23, 7.759197090219954e-24, 1.546886518887300e-24, 9.585471984274525e-26, -4.005276449706623e-26, -1.459299762834743e-26, -1.870491354620814e-27, 2.271573838802745e-28, 1.360979028415734e-28, 1.102172718357361e-29, -6.919347646474293e-30, 4.875282352118995e-31};
};
} // namespace polycoeff
} // namespace polytrope
} // namespace serif

0
src/poly/meson.build → src/polytrope/meson.build
View File

55
src/polytrope/solver/meson.build Normal file
View File

@@ -0,0 +1,55 @@
# ***********************************************************************
#
# Copyright (C) 2025 -- The 4D-STAR Collaboration
# File Author: Emily Boudreaux
# Last Modified: March 19, 2025
#
# 4DSSE is free software; you can use it and/or modify
# it under the terms and restrictions the GNU General Library Public
# License version 3 (GPLv3) as published by the Free Software Foundation.
#
# 4DSSE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General Public License
# along with this software; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# *********************************************************************** #
polySolver_sources = files(
'private/polySolver.cpp'
)
polySolver_headers = files(
'public/polySolver.h'
)
dependencies = [
mfem_dep,
meshio_dep,
polycoeff_dep,
polyutils_dep,
macros_dep,
probe_dep,
quill_dep,
config_dep,
resourceManager_dep,
types_dep,
]
libPolySolver = shared_library('polySolver',
polySolver_sources,
include_directories : include_directories('./public'),
cpp_args: ['-fvisibility=default'],
dependencies: dependencies,
install: true
)
polysolver_dep = declare_dependency(
include_directories : include_directories('./public'),
link_with : libPolySolver,
sources : polySolver_sources,
dependencies : dependencies
)

501
src/polytrope/solver/private/polySolver.cpp Normal file
View File

@@ -0,0 +1,501 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "polySolver.h"
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>
#include "mfem.hpp"
#include "4DSTARTypes.h"
#include "config.h"
#include "integrators.h"
#include "mfem.hpp"
#include "polytropeOperator.h"
#include "polyCoeff.h"
#include "probe.h"
#include "resourceManager.h"
#include "resourceManagerTypes.h"
#include "utilities.h"
#include "quill/LogMacros.h"
namespace serif {
namespace polytrope {
namespace laneEmden {
double a (const int k, const double n) { // NOLINT(*-no-recursion)
if ( k == 0 ) { return 1; }
if ( k == 1 ) { return 0; }
else { return -(c(k-2, n)/(std::pow(k, 2)+k)); }
}
double c(const int m, const double n) { // NOLINT(*-no-recursion)
if ( m == 0 ) { return std::pow(a(0, n), n); }
else {
double termOne = 1.0/(m*a(0, n));
double acc = 0;
for (int k = 1; k <= m; k++) {
acc += (k*n-m+k)*a(k, n)*c(m-k, n);
}
return termOne*acc;
}
}
double thetaSeriesExpansion(const double xi, const double n, const int order) {
double acc = 0;
for (int k = 0; k < order; k++) {
acc += a(k, n) * std::pow(xi, k);
}
return acc;
}
} // namespace laneEmden
PolySolver::PolySolver(mfem::Mesh& mesh, const double n, const double order)
: m_config(serif::config::Config::getInstance()), // Updated
m_logManager(serif::probe::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")),
m_polytropicIndex(n),
m_feOrder(order),
m_mesh(mesh) {
// Use feOrder - 1 for the RT space to satisfy Brezzi-Babuska condition
// for the H1 and RT [H(div)] spaces
m_fecH1 = std::make_unique<mfem::H1_FECollection>(m_feOrder, m_mesh.SpaceDimension());
m_fecRT = std::make_unique<mfem::RT_FECollection>(m_feOrder - 1, m_mesh.SpaceDimension());
m_feTheta = std::make_unique<mfem::FiniteElementSpace>(&m_mesh, m_fecH1.get());
m_fePhi = std::make_unique<mfem::FiniteElementSpace>(&m_mesh, m_fecRT.get());
m_theta = std::make_unique<mfem::GridFunction>(m_feTheta.get());
m_phi = std::make_unique<mfem::GridFunction>(m_fePhi.get());
assembleBlockSystem();
}
PolySolver::PolySolver(const double n, const double order)
: PolySolver(prepareMesh(n), n, order){}
mfem::Mesh& PolySolver::prepareMesh(const double n) {
if (n > 4.99 || n < 0.0) {
throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(n));
}
const serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
const serif::resource::types::Resource& genericResource = rm.getResource("mesh:polySphere");
const auto &meshIO = std::get<std::unique_ptr<serif::mesh::MeshIO>>(genericResource);
meshIO->LinearRescale(polycoeff::x1(n)); // Assumes polycoeff is now serif::polytrope::polycoeff
return meshIO->GetMesh();
}
PolySolver::~PolySolver() = default;
void PolySolver::assembleBlockSystem() {
mfem::Array<int> blockOffsets = computeBlockOffsets();
const std::unique_ptr<formBundle> forms = buildIndividualForms(blockOffsets);
// const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1.0);
// mfem::Array<int> thetaCenterDofs, phiCenterDofs;
// std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
// mfem::SparseMatrix& D_mat = forms->D->SpMat();
//
// for (int i = 0; i < phiCenterDofs.Size(); ++i)
// {
// const int dof_idx = phiCenterDofs[i];
// if (dof_idx >= 0 && dof_idx < D_mat.Height()) {
// D_mat(dof_idx, dof_idx) += penalty_param;
// }
// }
// --- Build the BlockOperator ---
m_polytropOperator = std::make_unique<PolytropeOperator>(
std::move(forms->M),
std::move(forms->Q),
std::move(forms->D),
std::move(forms->S),
std::move(forms->f),
blockOffsets);
}
mfem::Array<int> PolySolver::computeBlockOffsets() const {
mfem::Array<int> blockOffsets;
blockOffsets.SetSize(3);
blockOffsets[0] = 0;
blockOffsets[1] = m_feTheta->GetVSize(); // Get actual number of dofs *before* applying BCs
blockOffsets[2] = m_fePhi->GetVSize();
blockOffsets.PartialSum(); // Cumulative sum to get the offsets
return blockOffsets;
}
std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<int> &blockOffsets) {
// --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) ---
auto forms = std::make_unique<formBundle>(
std::make_unique<mfem::MixedBilinearForm>(m_fePhi.get(), m_feTheta.get()),
std::make_unique<mfem::MixedBilinearForm>(m_feTheta.get(), m_fePhi.get()),
std::make_unique<mfem::BilinearForm>(m_fePhi.get()),
std::make_unique<mfem::BilinearForm>(m_feTheta.get()),
std::make_unique<mfem::NonlinearForm>(m_feTheta.get())
);
// --- Add the integrators to the forms ---
forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); // M ∫∇ψ^θ·N^φ dV
forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); // Q ∫ψ^φ·∇N^θ dV
forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator()); // D ∫ψ^φ·N^φ dV
forms->S->AddDomainIntegrator(new mfem::DiffusionIntegrator()); // S ∫∇ψ^θ·∇N^θ dV
forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex)); // Assumes polyMFEMUtils is now serif::polytrope::polyMFEMUtils
// --- Assemble and Finalize the forms ---
assembleAndFinalizeForm(forms->M);
assembleAndFinalizeForm(forms->Q);
assembleAndFinalizeForm(forms->D);
assembleAndFinalizeForm(forms->S);
// Note: The NonlinearForm does not need to / cannot be finalized, as it is not a matrix form. Rather, the operator
// will evaluate the nonlinear form during the solve phase.
return forms;
}
void PolySolver::assembleAndFinalizeForm(auto &f) {
// This constructs / ensures the matrix representation for each form
// Assemble => Computes the local element matrices across the domain. Adds these to the global matrix . Adds these to the global matrix.
// Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted (CSR).
f->Assemble();
f->Finalize();
}
void PolySolver::solve() const {
// --- Set the initial guess for the solution ---
setInitialGuess();
setOperatorEssentialTrueDofs();
// --- Cast the GridFunctions to mfem::Vector ---
const auto thetaVec = static_cast<mfem::Vector>(*m_theta); // NOLINT(*-slicing)
const auto phiVec = static_cast<mfem::Vector>(*m_phi); // NOLINT(*-slicing)
// --- Finalize the operator ---
// Finalize with the initial state of theta for the initial jacobian calculation
m_polytropOperator->finalize(thetaVec);
// --- Broadcast initial condition to the full state vector ---
const mfem::Array<int>& full_block_offsets = m_polytropOperator->get_block_offsets();
mfem::Vector x_full(full_block_offsets.Last());
mfem::BlockVector x_full_block(x_full, full_block_offsets);
x_full_block.GetBlock(0) = thetaVec; // NOLINT(*-slicing)
x_full_block.GetBlock(1) = phiVec; // NOLINT(*-slicing)
// --- Extract only the free DOFs from the full state vector ---
const mfem::Array<int>& freeDofs = m_polytropOperator->get_free_dofs();
mfem::Vector x_free(m_polytropOperator->get_reduced_system_size());
x_full.GetSubVector(freeDofs, x_free); // Extract the free DOFs from the full vector
// --- Initialize RHS ---
mfem::Vector zero_rhs(m_polytropOperator->get_reduced_system_size());
zero_rhs = 0.0;
// --- Setup and run the Newton solver ---
const solverBundle sb = setupNewtonSolver();
sb.newton.Mult(zero_rhs, x_free);
// --- Reconstruct the full state vector from the reduced solution ---
mfem::BlockVector solution = m_polytropOperator->reconstruct_full_block_state_vector(x_free);
// --- Save and view an approximate 1D solution ---
saveAndViewSolution(solution);
}
serif::types::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
mfem::Array<int> theta_ess_tdof_list;
mfem::Array<int> phi_ess_tdof_list;
mfem::Array<int> thetaCenterDofs, phiCenterDofs; // phiCenterDofs are not used
mfem::Array<double> thetaCenterVals, phiCenterVals;
std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
thetaCenterVals.SetSize(thetaCenterDofs.Size());
// phiCenterVals.SetSize(phiCenterDofs.Size());
//
// phiCenterVals = 0.0;
thetaCenterVals = 1.0;
mfem::Array<int> ess_brd(m_mesh.bdr_attributes.Max());
ess_brd = 1;
mfem::Array<double> thetaSurfaceVals, phiSurfaceVals;
m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list);
m_fePhi->GetEssentialTrueDofs(ess_brd, phi_ess_tdof_list);
thetaSurfaceVals.SetSize(theta_ess_tdof_list.Size());
thetaSurfaceVals = 0.0;
phiSurfaceVals.SetSize(phi_ess_tdof_list.Size());
phiSurfaceVals = polycoeff::thetaSurfaceFlux(m_polytropicIndex); // Assumes polycoeff is now serif::polytrope::polycoeff
// combine the essential dofs with the center dofs
theta_ess_tdof_list.Append(thetaCenterDofs);
thetaSurfaceVals.Append(thetaCenterVals);
// phi_ess_tdof_list.Append(phiCenterDofs);
// phiSurfaceVals.Append(phiCenterVals);
serif::types::MFEMArrayPair thetaPair = std::make_pair(theta_ess_tdof_list, thetaSurfaceVals);
serif::types::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiSurfaceVals);
serif::types::MFEMArrayPairSet pairSet = std::make_pair(thetaPair, phiPair);
return pairSet;
}
std::pair<mfem::Array<int>, mfem::Array<int>> PolySolver::findCenterElement() const {
mfem::Array<int> thetaCenterDofs;
mfem::Array<int> phiCenterDofs;
// --- 1. Find the index of the single mesh vertex at the origin ---
int center_vertex_idx = -1;
constexpr double tol = 1e-9; // A small tolerance for floating point comparison
for (int i = 0; i < m_mesh.GetNV(); ++i) {
const double* vertex_coords = m_mesh.GetVertex(i);
if (std::abs(vertex_coords[0]) < tol &&
std::abs(vertex_coords[1]) < tol &&
std::abs(vertex_coords[2]) < tol) {
center_vertex_idx = i;
break; // Found it, assume there's only one.
}
}
if (center_vertex_idx == -1) {
MFEM_ABORT("Could not find the center vertex at [0,0,0]. Check mesh construction.");
}
// --- 2. Get Theta (H1) DoFs associated ONLY with that vertex ---
// CORRECTED: Use GetVertexDofs, not GetVDofs.
m_feTheta->GetVertexDofs(center_vertex_idx, thetaCenterDofs);
mfem::Array<int> central_element_ids;
// PERF: could probably move this to a member variable and populate during construction
mfem::Table* vertex_to_elements_table = m_mesh.GetVertexToElementTable();
vertex_to_elements_table->Finalize();
mfem::Array<int> element_ids;
vertex_to_elements_table->GetRow(center_vertex_idx, element_ids);
delete vertex_to_elements_table;
for (int i = 0; i < element_ids.Size(); ++i) {
int element_id = element_ids[i];
mfem::Array<int> tempDofs;
m_fePhi->GetElementDofs(element_id, tempDofs);
// decode negative dofs to their true, physical, dof indices
for (int j = 0; j < tempDofs.Size(); ++j) {
int dof = tempDofs[j];
if (dof < 0) {
dof = -dof - 1; // Convert to positive index
}
phiCenterDofs.Append(dof);
}
}
phiCenterDofs.Sort();
phiCenterDofs.Unique();
return std::make_pair(thetaCenterDofs, phiCenterDofs);
}
void PolySolver::setInitialGuess() const {
// --- Set the initial guess for the solution ---
mfem::FunctionCoefficient thetaInitGuess (
[this](const mfem::Vector &x) {
const double r = x.Norml2();
// const double radius = Probe::getMeshRadius(*m_mesh);
// const double u = 1/radius;
// return (-1.0/radius) * r + 1;
// return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not
return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 10);
}
);
mfem::VectorFunctionCoefficient phiSurfaceVectors (m_mesh.SpaceDimension(),
[this](const mfem::Vector &x, mfem::Vector &y) {
const double r = x.Norml2();
mfem::Vector xh(x);
xh /= r; // Normalize the vector
y.SetSize(m_mesh.SpaceDimension());
y = xh;
y *= polycoeff::thetaSurfaceFlux(m_polytropicIndex); // Assumes polycoeff is now serif::polytrope::polycoeff
}
);
// We want to apply specific boundary conditions to the surface
mfem::Array<int> ess_brd(m_mesh.bdr_attributes.Max());
ess_brd = 1;
// θ = 0 at surface
mfem::ConstantCoefficient surfacePotential(0);
m_theta->ProjectCoefficient(thetaInitGuess);
m_theta->ProjectBdrCoefficient(surfacePotential, ess_brd);
mfem::GradientGridFunctionCoefficient phiInitGuess (m_theta.get());
m_phi->ProjectCoefficient(phiInitGuess);
// Note that this will not result in perfect boundary conditions
// because it must maintain H(div) continuity, this is
// why inhomogenous boundary conditions enforcement is needed for φ
// This manifests in PolytropeOperator::Mult where we do not
// just zero out the essential dof elements in the residuals vector
// for φ; rather, we need to set this to something which will push the
// solver towards a more consistent answer (x_φ - target)
m_phi->ProjectBdrCoefficientNormal(phiSurfaceVectors, ess_brd);
auto [thetaCenterDofs, phiCenterDofs] = findCenterElement();
for (int i = 0; i < phiCenterDofs.Size(); ++i)
{
(*m_phi)(phiCenterDofs[i]) = 0.0;
}
if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
serif::probe::glVisView(*m_theta, m_mesh, "θ init");
serif::probe::glVisView(*m_phi, m_mesh, "φ init");
}
std::cout << "HERE" << std::endl;
}
void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) const {
mfem::BlockVector x_block(const_cast<mfem::BlockVector&>(state_vector), m_polytropOperator->get_block_offsets());
mfem::Vector& x_theta = x_block.GetBlock(0);
mfem::Vector& x_phi = x_block.GetBlock(1);
if (m_config.get<bool>("Poly:Output:View", false)) {
serif::probe::glVisView(x_theta, *m_feTheta, "θ Solution");
serif::probe::glVisView(x_phi, *m_fePhi, "ɸ Solution");
}
// --- Extract the Solution ---
if (m_config.get<bool>("Poly:Output:1D:Save", true)) {
const auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
auto derivSolPath = "d" + solutionPath;
const auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
const auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
const auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
const std::vector rayDirection = {rayCoLatitude, rayLongitude};
serif::probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath);
// Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath);
}
}
void PolySolver::setOperatorEssentialTrueDofs() const {
const serif::types::MFEMArrayPairSet ess_tdof_pair_set = getEssentialTrueDof();
m_polytropOperator->set_essential_true_dofs(ess_tdof_pair_set);
}
void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const {
newtonRelTol = m_config.get<double>("Poly:Solver:Newton:RelTol", 1.e-4);
newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1.e-6);
newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 10);
newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 3);
gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1.e-12);
gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1.e-12);
gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 200);
gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", -1);
LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);
}
void PolySolver::GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename) {
mfem::Mesh *mesh = fes.GetMesh();
double r = serif::probe::getMeshRadius(*mesh);
std::ofstream outputFile(filename, std::ios::out | std::ios::trunc);
outputFile << "dof,R,r,x,y,z" << '\n';
const int nElements = mesh->GetNE();
mfem::Vector coords;
mfem::IntegrationPoint ipZero;
double p[3] = {0.0, 0.0, 0.0};
int actual_idx;
ipZero.Set3(p);
for (int i = 0; i < nElements; i++) {
mfem::Array<int> elemDofs;
fes.GetElementDofs(i, elemDofs);
mfem::ElementTransformation* T = mesh->GetElementTransformation(i);
mfem::Vector physCoord(3);
T->Transform(ipZero, physCoord);
for (int dofID = 0; dofID < elemDofs.Size(); dofID++) {
if (elemDofs[dofID] < 0) {
actual_idx = -elemDofs[dofID] - 1;
} else {
actual_idx = elemDofs[dofID];
}
outputFile << actual_idx;
if (dofID != elemDofs.Size() - 1) {
outputFile << "|";
} else {
outputFile << ",";
}
}
outputFile << r << "," << physCoord.Norml2() << "," << physCoord[0] << "," << physCoord[1] << "," << physCoord[2] << '\n';
}
outputFile.close();
}
solverBundle PolySolver::setupNewtonSolver() const {
// --- Load configuration parameters ---
double newtonRelTol, newtonAbsTol, gmresRelTol, gmresAbsTol;
int newtonMaxIter, newtonPrintLevel, gmresMaxIter, gmresPrintLevel;
LoadSolverUserParams(newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel, gmresRelTol, gmresAbsTol,
gmresMaxIter, gmresPrintLevel);
solverBundle solver; // Use this solver bundle to ensure lifetime safety
solver.solver.SetRelTol(gmresRelTol);
solver.solver.SetAbsTol(gmresAbsTol);
solver.solver.SetMaxIter(gmresMaxIter);
solver.solver.SetPrintLevel(gmresPrintLevel);
// solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner());
// --- Set up the Newton solver ---
solver.newton.SetRelTol(newtonRelTol);
solver.newton.SetAbsTol(newtonAbsTol);
solver.newton.SetMaxIter(newtonMaxIter);
solver.newton.SetPrintLevel(newtonPrintLevel);
solver.newton.SetOperator(*m_polytropOperator);
// --- Created the linear solver which is used to invert the jacobian ---
solver.newton.SetSolver(solver.solver);
return solver;
}
} // namespace polytrope
} // namespace serif

615
src/polytrope/solver/public/polySolver.h Normal file
View File

@@ -0,0 +1,615 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include "mfem.hpp"
#include <memory>
#include <utility>
#include "integrators.h"
#include "4DSTARTypes.h"
#include "polytropeOperator.h"
#include "config.h"
#include "meshIO.h"
#include "probe.h"
#include "quill/Logger.h"
namespace serif {
namespace polytrope {
/**
* @brief Namespace for Lane-Emden equation related utility functions.
*
* Provides functions to compute coefficients and evaluate the series expansion
* solution to the Lane-Emden equation, which describes the structure of a
* spherically symmetric polytropic star.
* The Lane-Emden equation is given by:
* \f[
* \frac{1}{\xi^2} \frac{d}{d\xi} \left( \xi^2 \frac{d\theta}{d\xi} \right) = -\theta^n
* \f]
* where \f$\xi\f$ is a dimensionless radius and \f$\theta\f$ is related to the density,
* and \f$n\f$ is the polytropic index.
*/
namespace laneEmden {
/**
* @brief Computes the coefficient \f$a_k\f$ for the Lane-Emden series expansion.
*
* The series solution for \f$\theta(\xi)\f$ is given by \f$\theta(\xi) = \sum_{k=0}^{\infty} a_k \xi^k\f$.
* The coefficients \f$a_k\f$ are determined by substituting the series into the Lane-Emden equation.
* Specifically, \f$a_0 = 1\f$, \f$a_1 = 0\f$, and for \f$k \ge 2\f$,
* \f$a_k = -\frac{c_{k-2,n}}{k(k+1)}\f$.
*
* @param k The index of the coefficient.
* @param n The polytropic index.
* @return The value of the coefficient \f$a_k\f$.
* @see c(const int m, const double n)
*/
double a (const int k, const double n);
/**
* @brief Computes the auxiliary coefficient \f$c_{m,n}\f$ used in determining \f$a_k\f$.
*
* The term \f$\theta^n\f$ in the Lane-Emden equation can also be expanded as a series
* \f$\theta^n(\xi) = \sum_{m=0}^{\infty} c_{m,n} \xi^m\f$.
* The coefficients \f$c_{m,n}\f$ are related to \f$a_k\f$ by:
* \f$c_{0,n} = a_0^n\f$
* \f$c_{m,n} = \frac{1}{m a_0} \sum_{j=1}^{m} (j n - m + j) a_j c_{m-j,n}\f$ for \f$m > 0\f$.
*
* @param m The index of the coefficient.
* @param n The polytropic index.
* @return The value of the coefficient \f$c_{m,n}\f$.
* @see a(const int k, const double n)
*/
double c(const int m, const double n);
/**
* @brief Computes the Lane-Emden function \f$\theta(\xi)\f$ using a series expansion.
*
* Evaluates the series \f$\theta(\xi) = \sum_{k=0}^{\text{order}-1} a_k \xi^k\f$ up to a specified order.
* This provides an approximate solution to the Lane-Emden equation, particularly accurate for small \f$\xi\f$.
*
* @param xi The dimensionless radius \f$\xi\f$.
* @param n The polytropic index.
* @param order The number of terms to include in the series expansion.
* @return The approximate value of \f$\theta(\xi)\f$.
*
* @example
* @code
* double xi = 0.5;
* double n = 1.5;
* int series_order = 10;
* double theta_val = laneEmden::thetaSeriesExpansion(xi, n, series_order);
* // theta_val will be an approximation of the Lane-Emden function at xi=0.5 for n=1.5
* @endcode
*/
double thetaSeriesExpansion(const double xi, const double n, const int order);
} // namespace laneEmden
/**
* @brief Structure to manage the lifetime of MFEM solver objects.
*
* This structure ensures that the `mfem::GMRESSolver` outlives the `mfem::NewtonSolver`
* that might use it. The `mfem::NewtonSolver` often takes a pointer or reference to a
* linear solver, and if the linear solver is destroyed first, it can lead to dangling
* pointers and crashes.
*
* @note The order of declaration of members is crucial: `solver` must be declared
* before `newton` to ensure it is constructed first and destroyed last.
*/
struct solverBundle {
mfem::GMRESSolver solver; ///< The linear solver (e.g., GMRES). Must be declared first.
mfem::NewtonSolver newton; ///< The nonlinear solver (e.g., Newton). Must be declared second.
};
/**
* @brief Structure to hold the various bilinear and nonlinear forms for the polytrope problem.
*
* This structure bundles the MFEM forms that represent the discretized weak formulation
* of the mixed variable polytropic equations.
* The system being solved is typically:
* \f[
* \begin{aligned}
* \boldsymbol{\phi} + \nabla \theta &= \mathbf{0} \\
* \nabla \cdot \boldsymbol{\phi} - \theta^n &= 0
* \end{aligned}
* \f]
* After integration by parts and discretization, these lead to matrix equations involving
* the forms stored here.
*/
struct formBundle {
/**
* @brief Mixed bilinear form \f$ M(\psi^\theta, N^\phi) = \int_\Omega \nabla\psi^\theta \cdot N^\phi \,dV \f$.
* This form arises from the term \f$\nabla \theta\f$ in the first equation, tested with a vector test function \f$N^\phi\f$.
* It couples the scalar field \f$\theta\f$ (potential) with the vector field \f$\boldsymbol{\phi}\f$ (flux).
*/
std::unique_ptr<mfem::MixedBilinearForm> M;
/**
* @brief Mixed bilinear form \f$ Q(\psi^\phi, N^\theta) = \int_\Omega \psi^\phi \cdot \nabla N^\theta \,dV \f$.
* This form arises from the term \f$\boldsymbol{\phi}\f$ in the first equation, tested with a scalar test function \f$N^\theta\f$,
* after integration by parts of the \f$\nabla \theta\f$ term.
* It can also be seen as related to the transpose of a discrete gradient operator.
*/
std::unique_ptr<mfem::MixedBilinearForm> Q;
/**
* @brief Bilinear form \f$ D(\psi^\phi, N^\phi) = \int_\Omega \psi^\phi \cdot N^\phi \,dV \f$.
* This is a mass matrix for the vector field \f$\boldsymbol{\phi}\f$. It arises from the \f$\boldsymbol{\phi}\f$ term
* in the first equation when tested with a vector test function \f$N^\phi\f$.
*/
std::unique_ptr<mfem::BilinearForm> D;
/**
* @brief Bilinear form \f$ S(\psi^\theta, N^\theta) = \int_\Omega \nabla\psi^\theta \cdot \nabla N^\theta \,dV \f$.
* This is a stiffness matrix (Laplacian) for the scalar field \f$\theta\f$. It is used for stabilization terms
* or alternative formulations.
*/
std::unique_ptr<mfem::BilinearForm> S;
/**
* @brief Nonlinear form \f$ f(\theta)(\psi^\theta) = \int_\Omega \psi^\theta \cdot \theta^n \,dV \f$.
* This form arises from the nonlinear source term \f$\theta^n\f$ in the second equation, tested with a scalar
* test function \f$\psi^\theta\f$.
*/
std::unique_ptr<mfem::NonlinearForm> f;
};
/**
* @brief Solves the Lane-Emden equation for a polytropic star using a mixed finite element method.
*
* This class sets up and solves the system of equations describing a polytropic
* stellar model. The Lane-Emden equation, in its second-order form, is:
* \f[
* \frac{1}{\xi^2} \frac{d}{d\xi} \left( \xi^2 \frac{d\theta}{d\xi} \right) = -\theta^n
* \f]
* where \f$\xi\f$ is a dimensionless radius, \f$\theta\f$ is a dimensionless temperature/potential
* related to density (\f$\rho = \rho_c \theta^n\f$), and \f$n\f$ is the polytropic index.
*
* This solver uses a mixed formulation by introducing \f$\boldsymbol{\phi} = -\nabla \theta\f$:
* \f[
* \begin{aligned}
* \boldsymbol{\phi} + \nabla \theta &= \mathbf{0} \\
* \nabla \cdot \boldsymbol{\phi} - \theta^n &= 0
* \end{aligned}
* \f]
* These equations are discretized using H1 finite elements for \f$\theta\f$ and
* RT (Raviart-Thomas) finite elements for \f$\boldsymbol{\phi}\f$. The resulting nonlinear
* system is solved using a Newton method, with the `PolytropeOperator` class
* defining the system residual and Jacobian.
*
* Boundary conditions typically include \f$\theta(\xi_c) = 1\f$ at the center (\f$\xi_c \approx 0\f$),
* \f$\theta(\xi_1) = 0\f$ at the surface (\f$\xi_1\f$ is the first root of \f$\theta\f$), and
* \f$\boldsymbol{\phi} \cdot \mathbf{n} = -\frac{d\theta}{d\xi}\f$ at the surface, which is related to the
* surface gravity. At the center, \f$\boldsymbol{\phi}(\xi_c) = \mathbf{0}\f$ due to symmetry.
*/
class PolySolver final{
public: // Public methods
/**
* @brief Constructs a PolySolver instance.
*
* Initializes the solver for a given polytropic index \f$n\f$ and finite element order.
* This constructor internally calls `prepareMesh` to load and scale the mesh.
*
* @param n The polytropic index. Must be \f$0 \le n < 5\f$.
* @param order The polynomial order of the H1 finite elements for \f$\theta\f$.
* The RT elements for \f$\boldsymbol{\phi}\f$ will be of order `order-1`.
* @throw std::runtime_error If \f$n \ge 5\f$ or \f$n < 0\f$ (via `prepareMesh`).
*
* @example
* @code
* try {
* PolySolver solver(1.5, 2); // Polytropic index n=1.5, FE order 2 for theta
* solver.solve();
* } catch (const std::exception& e) {
* std::cerr << "Error: " << e.what() << std::endl;
* }
* @endcode
*/
PolySolver(const double n, const double order);
/**
* @brief Destructor for PolySolver.
*
* Defaulted, handles cleanup of owned resources like `std::unique_ptr` members.
*/
~PolySolver();
/**
* @brief Solves the polytropic system.
*
* This method orchestrates the solution process:
* 1. Sets the initial guess for \f$\theta\f$ and \f$\boldsymbol{\phi}\f$.
* 2. Applies essential boundary conditions to the `PolytropeOperator`.
* 3. Finalizes the `PolytropeOperator` with the initial state.
* 4. Sets up the Newton solver (including the GMRES linear solver for Jacobian systems).
* 5. Runs the Newton iteration to find the solution for the free DOFs.
* 6. Reconstructs the full solution vector from the free DOFs.
* 7. Saves and/or views the solution based on configuration.
*
* @throw std::runtime_error or `mfem::ErrorException` if the solver fails to converge or encounters numerical issues.
*/
void solve() const;
/**
* @brief Gets the polytropic index \f$n\f$.
* @return The polytropic index.
*/
double getN() const { return m_polytropicIndex; }
/**
* @brief Gets the finite element order used for the \f$\theta\f$ variable.
* @return The polynomial order of the H1 finite elements.
*/
double getOrder() const { return m_feOrder; }
/**
* @brief Gets a reference to the computational mesh.
* @return A reference to the `mfem::Mesh` object.
*/
mfem::Mesh& getMesh() const { return m_mesh; }
/**
* @brief Gets a reference to the solution grid function for \f$\theta\f$.
* @return A reference to the `mfem::GridFunction` storing the \f$\theta\f$ solution.
* @note The solution is populated after `solve()` has been successfully called.
*/
mfem::GridFunction& getSolution() const { return *m_theta; }
private: // Private Attributes
// --- Configuration and Logging ---
serif::config::Config& m_config; ///< Reference to the global configuration manager instance.
serif::probe::LogManager& m_logManager; ///< Reference to the global log manager instance.
quill::Logger* m_logger; ///< Pointer to the specific logger instance for this class.
// --- Physical and Discretization Parameters ---
double m_polytropicIndex; ///< The polytropic index \f$n\f$.
double m_feOrder; ///< The polynomial order for H1 elements (\f$\theta\f$). RT elements for \f$\boldsymbol{\phi}\f$ are of order `m_feOrder - 1`.
// --- MFEM Core Objects ---
mfem::Mesh& m_mesh; ///< Reference to the computational mesh (owned by ResourceManager).
std::unique_ptr<mfem::H1_FECollection> m_fecH1; ///< Finite Element Collection for \f$\theta\f$ (H1 elements).
std::unique_ptr<mfem::RT_FECollection> m_fecRT; ///< Finite Element Collection for \f$\boldsymbol{\phi}\f$ (Raviart-Thomas elements).
std::unique_ptr<mfem::FiniteElementSpace> m_feTheta; ///< Finite Element Space for \f$\theta\f$.
std::unique_ptr<mfem::FiniteElementSpace> m_fePhi; ///< Finite Element Space for \f$\boldsymbol{\phi}\f$.
// --- Solution Vectors ---
std::unique_ptr<mfem::GridFunction> m_theta; ///< Grid function for the scalar potential \f$\theta\f$.
std::unique_ptr<mfem::GridFunction> m_phi; ///< Grid function for the vector flux \f$\boldsymbol{\phi}\f$.
// --- Operator and Solver Components ---
std::unique_ptr<PolytropeOperator> m_polytropOperator; ///< The main nonlinear operator for the mixed system.
std::unique_ptr<mfem::OperatorJacobiSmoother> m_prec; ///< Preconditioner (currently seems unused in `polySolver.cpp`).
private: // Private methods
/**
* @brief Private constructor that takes an existing mesh.
*
* Initializes FE collections, spaces, grid functions, and assembles the block system.
* This is called by the public constructor after `prepareMesh`.
*
* @param mesh A reference to an initialized `mfem::Mesh`.
* @param n The polytropic index.
* @param order The polynomial order for H1 finite elements.
*/
PolySolver(mfem::Mesh& mesh, double n, double order);
/**
* @brief Prepares the mesh for the simulation.
*
* Loads a generic sphere mesh from the `ResourceManager` and scales it
* to the dimensionless radius \f$\xi_1(n)\f$ (the first zero of the Lane-Emden function
* for the given polytropic index \f$n\f$).
*
* @param n The polytropic index.
* @return A reference to the prepared `mfem::Mesh` (owned by `ResourceManager`).
* @throw std::runtime_error If \f$n \ge 5.0\f$ or \f$n < 0.0\f$, as \f$\xi_1(n)\f$ is typically
* undefined or infinite outside this range for physical polytropes.
* @throw std::runtime_error If the mesh resource "mesh:polySphere" cannot be found or loaded.
*
* @details The scaling factor \f$\xi_1(n)\f$ is obtained from `polycoeff::x1(n)`.
* The mesh is expected to be a unit sphere initially.
*/
static mfem::Mesh& prepareMesh(double n);
/**
* @brief Assembles the block system operator `m_polytropOperator`.
*
* This involves:
* 1. Computing block offsets for \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ variables.
* 2. Building individual bilinear and nonlinear forms (M, Q, D, S, f) using `buildIndividualForms`.
* 3. Constructing the `PolytropeOperator` with these forms and offsets.
*/
void assembleBlockSystem();
/**
* @brief Compute the block offsets for the operator. These are the offsets that define which dofs belong to which variable.
*
* @details
* Create the block offsets. These define the start of each block in the combined vector.
* Block offsets will be \f$[0, N_\theta, N_\theta + N_\phi]\f$, where \f$N_\theta\f$ is the number of
* degrees of freedom for \f$\theta\f$ and \f$N_\phi\f$ is for \f$\boldsymbol{\phi}\f$.
* The interpretation of this is that each block tells the operator where in the flattened (1D) vector
* the degrees of freedom or coefficients for that free parameter start and end. I.e.
* we know that in any flattened vector will have a size \f$N_\theta + N_\phi\f$. The \f$\theta\f$ dofs will span
* from `blockOffsets[0]` to `blockOffsets[1]` and the \f$\boldsymbol{\phi}\f$ dofs will span from `blockOffsets[1]` to `blockOffsets[2]`.
*
* This is the same for matrices only in 2D (rows and columns).
*
* The key point here is that this is fundamentally an accounting structure, it is here to keep track of what
* parts of vectors and matrices belong to which variable.
*
* Also note that we use `VSize` rather than `Size`. `Size` refers to the number of true dofs (after eliminating
* boundary conditions). `VSize` refers to the total number of dofs before BC elimination, which is needed here.
*
* @return `mfem::Array<int>` The offsets for the blocks in the operator.
*
* @pre `m_feTheta` and `m_fePhi` must be valid, populated `mfem::FiniteElementSpace` objects.
*/
mfem::Array<int> computeBlockOffsets() const;
/**
* @brief Build the individual forms for the block operator (M, Q, D, S, and f).
*
* @param blockOffsets The offsets for the blocks in the operator, computed by `computeBlockOffsets`.
* @return A `std::unique_ptr<formBundle>` containing the assembled forms.
*
* @note These forms are built exactly how they are defined in the derivation.
* For example, `Mform` corresponds to \f$+M\f$, not \f$-M\f$.
* The `PolytropeOperator` handles any necessary sign changes for the final system assembly.
*
* @details
* This method initializes and assembles the discrete forms corresponding to the weak formulation
* of the mixed polytropic system:
* \f[
* M = \int_\Omega \nabla\psi^\theta \cdot N^\phi \,dV \quad (\text{from } \nabla\theta \text{ term})
* \f]
* \f[
* Q = \int_\Omega \psi^\phi \cdot \nabla N^\theta \,dV \quad (\text{from } \boldsymbol{\phi} \text{ term, related to } -M^T)
* \f]
* \f[
* D = \int_\Omega \psi^\phi \cdot N^\phi \,dV \quad (\text{mass matrix for } \boldsymbol{\phi})
* \f]
* \f[
* S = \int_\Omega \nabla\psi^\theta \cdot \nabla N^\theta \,dV \quad (\text{stiffness matrix for } \theta \text{, for stabilization})
* \f]
* \f[
* f(\theta)(\psi^\theta) = \int_\Omega \psi^\theta \cdot \theta^n \,dV \quad (\text{from nonlinear term } \theta^n)
* \f]
* The `PolytropeOperator` will then use these to form a system like:
* \f[
* R(X) = \begin{pmatrix}
* \text{nonlin_op}(\theta) + M\,\boldsymbol{\phi} \\
* D\,\boldsymbol{\phi} - Q\,\theta
* \end{pmatrix}
* = \mathbf{0}
* \f]
* (The exact structure depends on the `PolytropeOperator`'s internal assembly, which might involve stabilization terms using S).
*
* @pre `m_feTheta` and `m_fePhi` must be valid, populated `mfem::FiniteElementSpace` objects.
* `m_polytropicIndex` must be set.
* @post The returned `formBundle` contains unique pointers to assembled (and finalized where appropriate) MFEM forms.
*/
std::unique_ptr<formBundle> buildIndividualForms(const mfem::Array<int>& blockOffsets);
/**
* @brief Assemble and finalize a given MFEM form.
*
* This template function calls `Assemble()` and `Finalize()` on the provided form.
* `Assemble()` computes local element matrices and adds them to the global matrix.
* `Finalize()` builds the sparsity pattern and allows the `SparseMatrix` representation to be extracted.
*
* @tparam FormType The type of the MFEM form (e.g., `mfem::BilinearForm`, `mfem::MixedBilinearForm`).
* @param f A reference to the form to be assembled and finalized.
*
* @pre `f` must be a valid form object that supports `Assemble()` and `Finalize()` methods.
* @post `f` is assembled and finalized, and its sparse matrix representation is available.
*
* @note Nonlinear forms like `mfem::NonlinearForm` are typically not "finalized" in this sense,
* as they don't directly produce a sparse matrix but are evaluated. This function is
* intended for linear forms.
*/
static void assembleAndFinalizeForm(auto &f);
/**
* @brief Computes the essential true degrees of freedom (DOFs) and their values for boundary conditions.
*
* This method determines the DOFs that correspond to:
* 1. \f$\theta = 1\f$ at the center of the star.
* 2. \f$\theta = 0\f$ on the surface of the star (boundary attribute 1).
* 3. \f$\boldsymbol{\phi} \cdot \mathbf{n} = \text{surface\_flux}\f$ on the surface, where the flux is derived from \f$\theta'(\xi_1)\f$.
* (This is applied to the normal component of \f$\boldsymbol{\phi}\f$).
* 4. Potentially \f$\boldsymbol{\phi} = \mathbf{0}\f$ at the center (though the current implementation in `polySolver.cpp`
* seems to set components of \f$\boldsymbol{\phi}\f$ to zero at DOFs associated with the center element(s)).
*
* @return An `serif::types::MFEMArrayPairSet` containing two pairs:
* - The first pair is for \f$\theta\f$: (`mfem::Array<int>` of essential TDOF indices, `mfem::Array<double>` of corresponding values).
* - The second pair is for \f$\boldsymbol{\phi}\f$: (`mfem::Array<int>` of essential TDOF indices, `mfem::Array<double>` of corresponding values).
*
* @details "True DOFs" (tdof) in MFEM refer to the actual degrees of freedom in the global system,
* as opposed to local DOFs within an element. Essential boundary conditions fix the values of these DOFs.
* The center condition for \f$\theta\f$ is applied to DOFs identified by `findCenterElement()`.
* Surface conditions are applied to boundary attributes marked as essential (typically attribute 1).
* The surface flux for \f$\boldsymbol{\phi}\f$ is obtained from `polycoeff::thetaSurfaceFlux(m_polytropicIndex)`.
*
* @pre `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized. `m_polytropicIndex` must be set.
* @throw `mfem::ErrorException` (via `MFEM_ABORT`) if `findCenterElement()` fails to locate the center vertex.
*/
serif::types::MFEMArrayPairSet getEssentialTrueDof() const;
/**
* @brief Finds the degrees of freedom (DOFs) associated with the geometric center (origin) of the mesh.
*
* This method identifies:
* 1. The mesh vertex located at coordinates (0,0,0).
* 2. The H1 DOFs (\f$\theta\f$) associated with this center vertex.
* 3. The RT DOFs (\f$\boldsymbol{\phi}\f$) associated with mesh elements that share this center vertex.
*
* These DOFs are used to apply boundary conditions at the center of the polytrope,
* such as \f$\theta(\xi_c)=1\f$ and \f$\boldsymbol{\phi}(\xi_c)=\mathbf{0}\f$.
*
* @return A `std::pair` of `mfem::Array<int>`:
* - `first`: Array of \f$\theta\f$ (H1) DOF indices at the center.
* - `second`: Array of \f$\boldsymbol{\phi}\f$ (RT) DOF indices associated with central elements.
*
* @throw `mfem::ErrorException` (via `MFEM_ABORT`) if no vertex is found at the origin (within tolerance).
*
* @details For RT elements, DOFs are typically associated with faces (or edges in 2D). This method collects
* DOFs from all elements connected to the center vertex. For H1 elements, DOFs can be directly
* associated with vertices.
*
* @pre `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized.
*/
std::pair<mfem::Array<int>, mfem::Array<int>> findCenterElement() const;
/**
* @brief Sets the initial guess for the solution variables \f$\theta\f$ and \f$\boldsymbol{\phi}\f$.
*
* - For \f$\theta\f$:
* - The interior initial guess is based on the Lane-Emden series expansion:
* \f$\theta(\xi) = \sum a_k \xi^k\f$ (using `laneEmden::thetaSeriesExpansion`).
* - The boundary condition \f$\theta(\xi_1) = 0\f$ is projected onto the surface.
* - For \f$\boldsymbol{\phi}\f$:
* - The initial guess is \f$\boldsymbol{\phi} = -\nabla \theta_{\text{init}}\f$, where \f$\theta_{\text{init}}\f$ is the initial guess for \f$\theta\f$.
* - The boundary condition for the normal component \f$\boldsymbol{\phi} \cdot \mathbf{n}\f$ at the surface is projected.
* This value is \f$\theta'(\xi_1)\f$, obtained from `polycoeff::thetaSurfaceFlux`.
* - \f$\boldsymbol{\phi}\f$ components corresponding to central DOFs (from `findCenterElement`) are set to 0.
*
* The method uses `mfem::GridFunction::ProjectCoefficient` and `ProjectBdrCoefficient` (or `ProjectBdrCoefficientNormal`)
* for these projections.
*
* @note The projection of boundary conditions for \f$\boldsymbol{\phi}\f$ (RT elements) might not be exact due to
* H(div) continuity requirements. Inhomogeneous BC enforcement in `PolytropeOperator` handles this more robustly.
*
* @pre `m_theta`, `m_phi`, `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized.
* `m_polytropicIndex` must be set.
* Configuration settings for viewing initial guess should be loaded if desired.
* @post `m_theta` and `m_phi` grid functions are populated with the initial guess.
*/
void setInitialGuess() const;
/**
* @brief Saves the 1D radial solution and optionally displays the 2D/3D solution using GLVis.
*
* Extracts the \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ components from the converged `state_vector`.
* - If configured (`Poly:Output:1D:Save`), it extracts a 1D solution along a specified ray
* (defined by co-latitude and longitude) using `Probe::getRaySolution` and saves it to a CSV file.
* - If configured (`Poly:Output:View`), it displays \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ using `Probe::glVisView`.
*
* @param state_vector The full solution vector (block vector containing \f$\theta\f$ and \f$\boldsymbol{\phi}\f$)
* obtained from the Newton solver, typically after reconstruction by `PolytropeOperator`.
*
* @pre `m_polytropOperator`, `m_feTheta`, `m_fePhi` must be initialized.
* Configuration settings for output and viewing must be loaded.
* @post Solution data is saved and/or visualized according to configuration.
*/
void saveAndViewSolution(const mfem::BlockVector& state_vector) const;
/**
* @brief Sets up the Newton solver and its associated linear solver (GMRES).
*
* 1. Loads solver parameters (tolerances, max iterations, print levels) for both Newton
* and GMRES from the configuration file (via `LoadSolverUserParams`).
* 2. Creates a `solverBundle` to manage the lifetimes of `mfem::GMRESSolver` and `mfem::NewtonSolver`.
* 3. Configures the GMRES solver with its parameters.
* 4. Configures the Newton solver with its parameters, sets the `PolytropeOperator` as the
* nonlinear system operator, and sets the configured GMRES solver as the linear solver
* for inverting Jacobians.
*
* @return A `solverBundle` struct containing the configured Newton and GMRES solvers.
* The ownership of the solvers within the bundle is managed by the bundle itself.
*
* @pre `m_polytropOperator` must be initialized and finalized.
* Configuration settings for solver parameters must be available.
* @post A fully configured `solverBundle` is returned, ready for `newton.Mult()`.
*/
solverBundle setupNewtonSolver() const;
/**
* @brief Sets the essential true DOFs on the `PolytropeOperator`.
*
* Calls `getEssentialTrueDof()` to compute the boundary condition DOFs and values,
* and then passes them to `m_polytropOperator->set_essential_true_dofs()`.
* This step is crucial for the `PolytropeOperator` to correctly handle
* boundary conditions when forming reduced systems and applying residuals/Jacobians.
*
* @pre `m_polytropOperator` must be initialized.
* `m_mesh`, `m_feTheta`, `m_fePhi`, `m_polytropicIndex` must be set (for `getEssentialTrueDof`).
* @post Essential boundary conditions are registered with the `m_polytropOperator`.
* The `PolytropeOperator` will likely be marked as not finalized after this call.
*/
void setOperatorEssentialTrueDofs() const;
/**
* @brief Loads Newton and GMRES solver parameters from the configuration.
*
* Retrieves relative tolerance, absolute tolerance, maximum iterations, and print level
* for both the Newton solver and the GMRES linear solver from the `Config` instance.
* Default values are used if specific configuration keys are not found.
*
* Configuration keys are typically prefixed with `Poly:Solver:Newton:` and `Poly:Solver:GMRES:`.
* Example keys: `Poly:Solver:Newton:RelTol`, `Poly:Solver:GMRES:MaxIter`.
*
* @param[out] newtonRelTol Relative tolerance for Newton solver.
* @param[out] newtonAbsTol Absolute tolerance for Newton solver.
* @param[out] newtonMaxIter Maximum iterations for Newton solver.
* @param[out] newtonPrintLevel Print level for Newton solver.
* @param[out] gmresRelTol Relative tolerance for GMRES solver.
* @param[out] gmresAbsTol Absolute tolerance for GMRES solver.
* @param[out] gmresMaxIter Maximum iterations for GMRES solver.
* @param[out] gmresPrintLevel Print level for GMRES solver.
*
* @pre `m_config` must be a valid reference to a `Config` object.
* `m_logger` should be initialized if debug logging is enabled.
* @post Output parameters are populated with values from configuration or defaults.
* Solver parameters are logged at DEBUG level.
*/
void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const;
/**
* @brief Utility function to get and save the coordinates of degrees of freedom for a finite element space.
*
* For each element in the mesh:
* 1. Gets the DOFs associated with that element from the given `fes`.
* 2. Gets the physical coordinates of the element's center (or a reference point).
* 3. Writes the DOF indices, mesh radius, element center radius, and element center x,y,z coordinates
* to the specified output file in CSV format.
*
* This can be useful for debugging or analyzing the distribution of DOFs.
*
* @param fes The `mfem::FiniteElementSpace` for which to get DOF coordinates.
* @param filename The name of the output CSV file.
*
* @pre `fes` must be a valid, initialized `mfem::FiniteElementSpace`.
* The mesh associated with `fes` must be valid.
* @post A CSV file named `filename` is created/truncated and populated with DOF coordinate information.
*
* @note For DOFs shared by multiple elements, this function might list them multiple times,
* associated with each element they belong to. The output format lists all DOFs for an element
* on one line, pipe-separated, followed by coordinate data for that element.
*/
static void GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename);
};
} // namespace polytrope
} // namespace serif

6
src/poly/utils/meson.build → src/polytrope/utils/meson.build
View File

@@ -20,7 +20,8 @@
# *********************************************************************** #
polyutils_sources = files(
'private/integrators.cpp',
'private/operator.cpp'
'private/polytropeOperator.cpp',
'private/utilities.cpp',
)
dependencies = [
@@ -30,9 +31,10 @@ dependencies = [
quill_dep,
config_dep,
types_dep,
mfemanalysis_dep,
]
libpolyutils = library('polyutils',
libpolyutils = shared_library('polyutils',
polyutils_sources,
include_directories : include_directories('./public'),
cpp_args: ['-fvisibility=default'],

139
src/polytrope/utils/private/integrators.cpp Normal file
View File

@@ -0,0 +1,139 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "mfem.hpp"
#include <cmath>
#include "integrators.h"
#include "config.h"
#include <string>
namespace serif::polytrope::polyMFEMUtils {
NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
m_polytropicIndex(n),
m_epsilon(serif::config::Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
if (m_polytropicIndex < 0.0) {
throw std::invalid_argument("Polytropic index must be non-negative.");
}
if (m_polytropicIndex > 5.0) {
throw std::invalid_argument("Polytropic index must be less than 5.0.");
}
if (m_epsilon <= 0.0) {
throw std::invalid_argument("Epsilon (Poly:Solver:Epsilon) must be positive non-zero.");
}
}
void NonlinearPowerIntegrator::AssembleElementVector(
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::Vector &elvect) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elvect.SetSize(dof);
elvect = 0.0;
mfem::Vector shape(dof);
mfem::Vector physCoord;
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
mfem::IntegrationPoint ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
double u_nl;
Trans.Transform(ip, physCoord);
const double r = physCoord.Norml2();
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
} else {
u_nl = std::pow(u_val, m_polytropicIndex);
}
} else {
u_nl = 1.0 - m_polytropicIndex * m_regularizationCoeff * std::pow(r, 2);
}
for (int i = 0; i < dof; i++){
elvect(i) += shape(i) * u_nl * weight;
}
}
}
void NonlinearPowerIntegrator::AssembleElementGrad (
const mfem::FiniteElement &el,
mfem::ElementTransformation &Trans,
const mfem::Vector &elfun,
mfem::DenseMatrix &elmat) {
const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
const int dof = el.GetDof();
elmat.SetSize(dof);
elmat = 0.0;
mfem::Vector shape(dof);
mfem::DenseMatrix dshape(dof, 3);
mfem::DenseMatrix invJ(3, 3);
mfem::Vector physCoord;
for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
Trans.SetIntPoint(&ip);
const double weight = ip.weight * Trans.Weight();
Trans.Transform(ip, physCoord);
double r = physCoord.Norml2();
el.CalcShape(ip, shape);
double u_val = 0.0;
for (int j = 0; j < dof; j++) {
u_val += elfun(j) * shape(j);
}
double d_u_nl;
if (r > m_regularizationRadius) {
if (u_val < m_epsilon) {
d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
} else {
d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
}
} else {
d_u_nl = 0.0;
}
for (int i = 0; i < dof; i++) {
for (int j = 0; j < dof; j++) {
elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
}
}
}
}
}

493
src/polytrope/utils/private/polytropeOperator.cpp Normal file
View File

@@ -0,0 +1,493 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "polytropeOperator.h"
#include "4DSTARTypes.h"
#include "utilities.h"
#include "mfem.hpp"
#include "mfem_smout.h"
#include <memory>
#include "config.h"
namespace serif {
namespace polytrope {
// --- SchurCompliment Class Implementation ---
// SchurCompliment: Constructors
SchurCompliment::SchurCompliment(
const mfem::SparseMatrix &QOp,
const mfem::SparseMatrix &DOp,
const mfem::SparseMatrix &MOp,
const mfem::Solver &GradInvOp ) :
mfem::Operator(DOp.Height(), DOp.Width()) {
// Initialize sizes
SetOperator(QOp, DOp, MOp, GradInvOp);
m_nPhi = m_DOp->Height();
m_nTheta = m_MOp->Height();
}
SchurCompliment::SchurCompliment(
const mfem::SparseMatrix &QOp,
const mfem::SparseMatrix &DOp,
const mfem::SparseMatrix &MOp) :
mfem::Operator(DOp.Height(), DOp.Width())
{
updateConstantTerms(QOp, DOp, MOp);
m_nPhi = m_DOp->Height();
m_nTheta = m_MOp->Height();
}
// SchurCompliment: Public Interface Methods
void SchurCompliment::Mult(const mfem::Vector &x, mfem::Vector &y) const {
// Check that the input vector is the correct size
if (x.Size() != m_nPhi) {
MFEM_ABORT("Input vector x has size " + std::to_string(x.Size()) + ", expected " + std::to_string(m_nPhi));
}
if (y.Size() != m_nPhi) {
MFEM_ABORT("Output vector y has size " + std::to_string(y.Size()) + ", expected " + std::to_string(m_nPhi));
}
// Check that the operators are set
if (m_QOp == nullptr) {
MFEM_ABORT("QOp is null in SchurCompliment::Mult");
}
if (m_DOp == nullptr) {
MFEM_ABORT("DOp is null in SchurCompliment::Mult");
}
if (m_MOp == nullptr) {
MFEM_ABORT("MOp is null in SchurCompliment::Mult");
}
if (m_GradInvOp == nullptr) {
MFEM_ABORT("GradInvOp is null in SchurCompliment::Mult");
}
mfem::Vector v1(m_nTheta); // M * x
m_MOp -> Mult(x, v1); // M * x
mfem::Vector v2(m_nTheta); // GradInv * M * x
m_GradInvOp -> Mult(v1, v2); // GradInv * M * x
mfem::Vector v3(m_nPhi); // Q * GradInv * M * x
m_QOp -> Mult(v2, v3); // Q * GradInv * M * x
mfem::Vector v4(m_nPhi); // D * x
m_DOp -> Mult(x, v4); // D * x
subtract(v4, v3, y); // (D - Q * GradInv * M) * x
}
void SchurCompliment::SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp) {
updateConstantTerms(QOp, DOp, MOp);
updateInverseNonlinearJacobian(GradInvOp);
}
void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv) {
m_GradInvOp = &gradInv;
}
// SchurCompliment: Private Helper Methods
void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) {
m_QOp = &QOp;
m_DOp = &DOp;
m_MOp = &MOp;
}
// --- GMRESInverter Class Implementation ---
// GMRESInverter: Constructor
GMRESInverter::GMRESInverter(const SchurCompliment &op) :
mfem::Operator(op.Height(), op.Width()),
m_op(op) {
m_solver.SetOperator(m_op);
m_solver.SetMaxIter(100);
m_solver.SetRelTol(1e-1);
m_solver.SetAbsTol(1e-1);
}
// GMRESInverter: Public Interface Methods
void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const {
m_solver.Mult(x, y); // Approximates m_op^-1 * x
}
// --- PolytropeOperator Class Implementation ---
// PolytropeOperator: Constructor
PolytropeOperator::PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> M,
std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::BilinearForm> S,
std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets) :
// TODO: Need to update this so that the size is that of the reduced system operator
mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector
m_blockOffsets(blockOffsets),
m_jacobian(nullptr) {
m_M = std::move(M);
m_Q = std::move(Q);
m_D = std::move(D);
m_S = std::move(S);
m_f = std::move(f);
// Use Gauss-Seidel smoother to approximate the inverse of the matrix
// t = 0 (symmetric Gauss-Seidel), 1 (forward Gauss-Seidel), 2 (backward Gauss-Seidel)
// iterations = 3
m_invNonlinearJacobian = std::make_unique<mfem::GSSmoother>(0, 3);
}
// PolytropeOperator: Core Operator Overrides
void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::Mult called before finalize");
}
// TODO: confirm that the vectors xFree and m_freeDofs are always parallel
m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector
mfem::Vector y;
y.SetSize(m_blockOffsets.Last());
// -- Create BlockVector views for input x and output y
mfem::BlockVector x_block(const_cast<mfem::Vector&>(m_state), m_blockOffsets);
mfem::BlockVector y_block(y, m_blockOffsets);
// -- Get Vector views for individual blocks
const mfem::Vector &x_theta = x_block.GetBlock(0);
const mfem::Vector &x_phi = x_block.GetBlock(1);
mfem::Vector &y_R0 = y_block.GetBlock(0); // Residual Block 0 (theta)
mfem::Vector &y_R1 = y_block.GetBlock(1); // Residual Block 1 (phi)
int theta_size = m_blockOffsets[1] - m_blockOffsets[0];
int phi_size = m_blockOffsets[2] - m_blockOffsets[1];
mfem::Vector f_term(theta_size);
mfem::Vector Mphi_term(theta_size);
mfem::Vector Dphi_term(phi_size);
mfem::Vector Qtheta_term(phi_size);
mfem::Vector Stheta_term(theta_size);
// Calculate R0 and R1 terms
// R0 = f(θ) + (1+c)Mɸ + cSθ
// R1 = (1+c)Dɸ - (1+c)Qθ
MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult");
m_f->Mult(x_theta, f_term); // f(θ)
m_M->Mult(x_phi, Mphi_term); // Mɸ
m_D->Mult(x_phi, Dphi_term); // Dɸ
m_Q->Mult(x_theta, Qtheta_term); // Qθ
m_S->Mult(x_theta, Stheta_term); // Sθ
// PERF: these multiplications may be able to be refactored into the matrices so they only need to be done once.
Stheta_term *= m_stabilizationCoefficient; // cSθ
Qtheta_term *= m_IncrementedStabilizationCoefficient; // (1+c)Qθ
Mphi_term *= m_IncrementedStabilizationCoefficient; // (1+c)Mɸ
Dphi_term *= m_IncrementedStabilizationCoefficient; // (1+c)Dɸ
mfem::Vector y_R0_temp(theta_size);
add(f_term, Mphi_term, y_R0_temp); // R0 = f(θ) + (1+c)Mɸ
add(y_R0_temp, Stheta_term, y_R0); // R0 = f(θ) + (1+c)Mɸ + cSθ
subtract(Dphi_term, Qtheta_term, y_R1); // R1 = (1+c)Dɸ - (1+c)Qθ
// --- Scatter the residual vector y into the free dofs ---
yFree.SetSize(m_reducedBlockOffsets.Last());
MFEM_ASSERT(m_freeDofs.Size() == m_reducedBlockOffsets.Last(), "PolytropeOperator::Mult: Size of free dofs does not match reduced block offsets size.");
for (int i = 0, j = 0; i < y.Size(); ++i) {
if (m_freeDofs.Find(i) != -1) {
yFree[j] = y[i];
j++;
}
}
}
mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::GetGradient called before finalize");
}
m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector
// --- Get the gradient of f ---
mfem::BlockVector x_block(const_cast<mfem::Vector&>(m_state), m_blockOffsets);
const mfem::Vector& x_theta = x_block.GetBlock(0);
// PERF: There are a lot of copies and loops here, probably performance could be gained by flattering some of these.
auto &grad = m_f->GetGradient(x_theta);
// updatePreconditioner(grad);
const auto gradMatrix = dynamic_cast<mfem::SparseMatrix*>(&grad); // ∂f(θ)/∂θ
if (gradMatrix == nullptr) {
MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix.");
}
m_gradReduced = std::make_unique<mfem::SparseMatrix> (
serif::utilities::build_reduced_matrix(
*gradMatrix,
m_theta_ess_tdofs.first,
m_theta_ess_tdofs.first
)
); // ∂f(θ)/∂θ (now reduced to only the free DOFs)
// TODO: Confirm this actually does what I want it to do
*m_gradReduced += *m_ScaledSReduced; // ∂f(θ)/∂θ + cS (reduced to only the free DOFs)
m_jacobian->SetBlock(0, 0, m_gradReduced.get());
return *m_jacobian;
}
// PolytropeOperator: Setup and Finalization
void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
using serif::utilities::build_reduced_matrix;
if (m_isFinalized) {
return;
}
// These functions must be called in this order since they depend on each others post state
// TODO: Refactor this so that either there are explicit checks to make sure the order is correct or make
// them pure functions
construct_matrix_representations();
construct_reduced_block_offsets();
construct_jacobian_constant_terms();
scatter_boundary_conditions();
populate_free_dof_array();
// Override the size based on the reduced system
height = m_reducedBlockOffsets.Last();
width = m_reducedBlockOffsets.Last();
m_isFinalized = true;
}
// PolytropeOperator: Essential True DOF Management
void PolytropeOperator::set_essential_true_dofs(const serif::types::MFEMArrayPair& theta_ess_tdofs, const serif::types::MFEMArrayPair& phi_ess_tdofs) {
m_isFinalized = false;
m_theta_ess_tdofs = theta_ess_tdofs;
m_phi_ess_tdofs = phi_ess_tdofs;
if (m_f) {
m_f->SetEssentialTrueDofs(theta_ess_tdofs.first);
} else {
MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
}
}
void PolytropeOperator::set_essential_true_dofs(const serif::types::MFEMArrayPairSet& ess_tdof_pair_set) {
set_essential_true_dofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second);
}
serif::types::MFEMArrayPairSet PolytropeOperator::get_essential_true_dofs() const {
return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs);
}
// PolytropeOperator: Getter Methods
const mfem::BlockOperator &PolytropeOperator::get_jacobian_operator() const {
if (m_jacobian == nullptr) {
MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call.");
}
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetJacobianOperator().");
}
return *m_jacobian;
}
mfem::BlockDiagonalPreconditioner& PolytropeOperator::get_preconditioner() const {
if (m_schurPreconditioner == nullptr) {
MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call.");
}
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner().");
}
return *m_schurPreconditioner;
}
int PolytropeOperator::get_reduced_system_size() const {
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator not finalized prior to call to GetReducedSystemSize().");
}
return m_reducedBlockOffsets.Last();
}
// PolytropeOperator: State Reconstruction
const mfem::Vector &PolytropeOperator::reconstruct_full_state_vector(const mfem::Vector &reducedState) const {
m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector
return m_state;
}
const mfem::BlockVector PolytropeOperator::reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const {
m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector
mfem::BlockVector x_block(m_state, m_blockOffsets);
return x_block;
}
// PolytropeOperator: DOF Population
void PolytropeOperator::populate_free_dof_array() {
m_freeDofs.SetSize(0);
for (int i = 0; i < m_blockOffsets.Last(); i++) {
const int thetaSearchIndex = i;
const int phiSearchIndex = i - m_blockOffsets[1];
if (phiSearchIndex < 0){
if (m_theta_ess_tdofs.first.Find(thetaSearchIndex) == -1) {
m_freeDofs.Append(i);
}
} else {
if (m_phi_ess_tdofs.first.Find(phiSearchIndex) == -1) {
m_freeDofs.Append(i);
}
}
}
}
// PolytropeOperator: Private Helper Methods - Construction and Setup
void PolytropeOperator::construct_matrix_representations() {
// --- Construct the sparse matrix representations of M, Q, D, and S ---
m_Mmat = std::make_unique<mfem::SparseMatrix>(m_M->SpMat());
m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
m_Smat = std::make_unique<mfem::SparseMatrix>(m_S->SpMat());
// --- Reduce the matrices to only the free DOFs ---
m_MReduced = std::make_unique<mfem::SparseMatrix>(
serif::utilities::build_reduced_matrix(
*m_Mmat,
m_phi_ess_tdofs.first,
m_theta_ess_tdofs.first)
);
m_QReduced = std::make_unique<mfem::SparseMatrix>(
serif::utilities::build_reduced_matrix(
*m_Qmat,
m_theta_ess_tdofs.first,
m_phi_ess_tdofs.first)
);
m_DReduced = std::make_unique<mfem::SparseMatrix>(
serif::utilities::build_reduced_matrix(
*m_Dmat,
m_phi_ess_tdofs.first,
m_phi_ess_tdofs.first)
);
m_SReduced = std::make_unique<mfem::SparseMatrix>(
serif::utilities::build_reduced_matrix(
*m_Smat,
m_theta_ess_tdofs.first,
m_theta_ess_tdofs.first)
);
// --- Scale the reduced matrices by the stabilization coefficients ---
mfem::SparseMatrix MScaledTemp(*m_MReduced); // Create a temporary copy of the M matrix for scaling
MScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the M matrix by the incremented stabilization coefficient
m_MScaledReduced = std::make_unique<mfem::SparseMatrix>(MScaledTemp); // Store the scaled M matrix so that it persists
mfem::SparseMatrix QScaledTemp(*m_QReduced);
QScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the Q matrix by the incremented stabilization coefficient
m_QScaledReduced = std::make_unique<mfem::SparseMatrix>(QScaledTemp); // Store the scaled Q matrix so that it persists
mfem::SparseMatrix DRScaledTemp(*m_DReduced);
DRScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the D matrix by the incremented stabilization coefficient
m_DScaledReduced = std::make_unique<mfem::SparseMatrix>(DRScaledTemp); // Store the scaled D matrix so that it persists
// Scale the S matrix by the stabilization coefficient (so that the allocations only need to be done once)
mfem::SparseMatrix SScaledTemp(*m_SReduced); // Create a temporary copy of the S matrix for scaling
SScaledTemp *= m_stabilizationCoefficient; // Scale the S matrix by the stabilization coefficient
m_ScaledSReduced = std::make_unique<mfem::SparseMatrix>(SScaledTemp); // Store the scaled S matrix so that it persists
// --- Create the scaled operator for -(1+c)Q ---
m_negQ_mat = std::make_unique<mfem::ScaledOperator>(m_QScaledReduced.get(), -1.0);
}
void PolytropeOperator::construct_reduced_block_offsets() {
m_reducedBlockOffsets.SetSize(3);
m_reducedBlockOffsets[0] = 0; // R0 block (theta)
m_reducedBlockOffsets[1] = m_MReduced->Height(); // R1 block (theta)
m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi + theta)
}
void PolytropeOperator::construct_jacobian_constant_terms() {
m_jacobian = std::make_unique<mfem::BlockOperator>(m_reducedBlockOffsets);
m_jacobian->SetBlock(0, 1, m_MScaledReduced.get()); //<- (1+c)M (constant, reduced to only free DOFs)
m_jacobian->SetBlock(1, 0, m_negQ_mat.get()); //<- -(1+c)Q (constant, reduced to only free DOFs)
m_jacobian->SetBlock(1, 1, m_DScaledReduced.get()); //<- (1+c)D (constant, reduced to only free DOFs)
}
void PolytropeOperator::scatter_boundary_conditions() {
mfem::Vector thetaStateValues(m_theta_ess_tdofs.first.Size());
for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
thetaStateValues[i] = m_theta_ess_tdofs.second[i];
}
mfem::Vector phiStateValues(m_phi_ess_tdofs.first.Size());
for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
phiStateValues[i] = m_phi_ess_tdofs.second[i]; // TODO: figure out if this needs to be normalized
}
mfem::Array<int> phiDofIndices(m_phi_ess_tdofs.first.Size());
for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
phiDofIndices[i] = m_phi_ess_tdofs.first[i] + m_blockOffsets[1];
}
m_state.SetSize(m_blockOffsets.Last());
m_state = 0.0;
m_state.SetSubVector(m_theta_ess_tdofs.first, thetaStateValues);
m_state.SetSubVector(phiDofIndices, phiStateValues);
}
// PolytropeOperator: Private Helper Methods - Jacobian and Preconditioner Updates
void PolytropeOperator::update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const {
m_invNonlinearJacobian->SetOperator(grad);
}
void PolytropeOperator::update_inverse_schur_compliment() const {
// TODO: This entire function could probably be refactored out
if (!m_isFinalized) {
MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before finalize");
}
if (m_invNonlinearJacobian == nullptr) {
MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
}
if (m_schurCompliment == nullptr) {
MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseSchurCompliment");
}
m_schurCompliment->updateInverseNonlinearJacobian(*m_invNonlinearJacobian);
if (m_schurPreconditioner == nullptr) {
m_schurPreconditioner = std::make_unique<mfem::BlockDiagonalPreconditioner>(m_blockOffsets);
}
m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get());
m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get());
}
void PolytropeOperator::update_preconditioner(const mfem::Operator &grad) const {
update_inverse_nonlinear_jacobian(grad);
update_inverse_schur_compliment();
}
} // namespace polytrope
} // namespace serif

132
src/polytrope/utils/private/utilities.cpp Normal file
View File

@@ -0,0 +1,132 @@
#include "utilities.h"
#include "mfem.hpp"
#include <memory>
namespace serif::utilities {
mfem::SparseMatrix build_reduced_matrix(
const mfem::SparseMatrix& matrix,
const mfem::Array<int>& trialEssentialDofs,
const mfem::Array<int>& testEssentialDofs
) {
int M_orig = matrix.Height();
int N_orig = matrix.Width();
// 1. Create boolean lookup tables for eliminated rows/columns
// These tables help quickly check if an original row/column index is eliminated.
mfem::Array<bool> row_is_eliminated(M_orig);
row_is_eliminated = false; // Initialize all to false (no rows eliminated yet)
for (int i = 0; i < testEssentialDofs.Size(); ++i) {
int r_idx = testEssentialDofs[i];
if (r_idx >= 0 && r_idx < M_orig) { // Check for valid index
row_is_eliminated[r_idx] = true;
}
}
mfem::Array<bool> col_is_eliminated(N_orig);
col_is_eliminated = false; // Initialize all to false (no columns eliminated yet)
for (int i = 0; i < trialEssentialDofs.Size(); ++i) {
int c_idx = trialEssentialDofs[i];
if (c_idx >= 0 && c_idx < N_orig) { // Check for valid index
col_is_eliminated[c_idx] = true;
}
}
// 2. Create mappings from old (original) indices to new indices
// Also, count the number of rows and columns in the new, reduced matrix.
mfem::Array<int> old_row_to_new_row(M_orig);
int num_new_rows = 0;
for (int i = 0; i < M_orig; ++i) {
if (row_is_eliminated[i]) {
old_row_to_new_row[i] = -1; // Mark as eliminated (no corresponding new row)
} else {
old_row_to_new_row[i] = num_new_rows++; // Assign the next available new row index
}
}
mfem::Array<int> old_col_to_new_col(N_orig);
int num_new_cols = 0;
for (int i = 0; i < N_orig; ++i) {
if (col_is_eliminated[i]) {
old_col_to_new_col[i] = -1; // Mark as eliminated (no corresponding new column)
} else {
old_col_to_new_col[i] = num_new_cols++; // Assign the next available new column index
}
}
// 3. Create the new sparse matrix with the calculated reduced dimensions.
// It's initially empty (no non-zero entries).
mfem::SparseMatrix A_new(num_new_rows, num_new_cols);
// 4. Iterate through the non-zero entries of the original matrix (A_orig).
// A_orig is typically stored in Compressed Sparse Row (CSR) format.
// GetI() returns row pointers, GetJ() returns column indices, GetData() returns values.
const int* I_orig = matrix.GetI();
const int* J_orig = matrix.GetJ();
const double* V_orig = matrix.GetData(); // Assuming double, can be templated if needed
for (int r_orig = 0; r_orig < M_orig; ++r_orig) {
// If the original row is marked for elimination, skip all its entries.
if (row_is_eliminated[r_orig]) {
continue;
}
// Get the new row index for the current original row.
int r_new = old_row_to_new_row[r_orig];
// Iterate through non-zero entries in the current original row r_orig.
// I_orig[r_orig] is the start index in J_orig and V_orig for row r_orig.
// I_orig[r_orig+1]-1 is the end index.
for (int k = I_orig[r_orig]; k < I_orig[r_orig + 1]; ++k) {
int c_orig = J_orig[k]; // Original column index of the non-zero entry
double val = V_orig[k]; // Value of the non-zero entry
if (col_is_eliminated[c_orig]) {
continue;
}
int c_new = old_col_to_new_col[c_orig];
A_new.Add(r_new, c_new, val);
}
}
A_new.Finalize();
return A_new;
}
mfem::Vector build_dof_identification_vector(const mfem::Array<int>& allDofs, const::mfem::Array<int>& highlightDofs) {
mfem::Vector v(allDofs.Size());
v = 0.0; // Initialize the vector to zero
v.SetSubVector(highlightDofs, 1.0); // Set the highlighted dofs to 1.0
return v;
}
mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf) {
mfem::Mesh* mesh = phi_gf.FESpace()->GetMesh();
const int dim = mesh->Dimension();
// Match the polynomial order of the original RT space for consistency.
const int order = phi_gf.FESpace()->GetOrder(0);
mfem::Vector curl_mag_vec;
for (int ne = 0; ne < mesh->GetNE(); ++ne) {
if (mesh->GetElementType(ne) != mfem::Element::TRIANGLE &&
mesh->GetElementType(ne) != mfem::Element::QUADRILATERAL &&
mesh->GetElementType(ne) != mfem::Element::TETRAHEDRON &&
mesh->GetElementType(ne) != mfem::Element::HEXAHEDRON) {
throw std::invalid_argument("Mesh element type not supported for curl computation.");
}
mfem::IsoparametricTransformation T;
mesh->GetElementTransformation(ne, &T);
phi_gf.GetCurl(T, curl_mag_vec);
std::cout << "HERE" << std::endl;
}
mfem::L2_FECollection fac(order, dim);
mfem::FiniteElementSpace fs(mesh, &fac);
mfem::GridFunction curl_gf(&fs);
curl_gf = 0.0;
return curl_gf;
}
}

38
src/poly/utils/public/integrators.h → src/polytrope/utils/public/integrators.h
View File

@@ -31,7 +31,8 @@
* @brief A collection of utilities for working with MFEM and solving the lane-emden equation.
*/
namespace serif {
namespace polytrope {
/**
* @namespace polyMFEMUtils
* @brief A namespace for utilities for working with MFEM and solving the lane-emden equation.
@@ -69,9 +70,38 @@ namespace polyMFEMUtils {
*/
virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
private:
Config& m_config = Config::getInstance();
Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
double m_polytropicIndex;
double m_epsilon;
static constexpr double m_regularizationRadius = 0.15; ///< Regularization radius for the epsilon function, used to avoid singularities in the power law.
static constexpr double m_regularizationCoeff = 1.0/6.0; ///< Coefficient for the regularization term, used to ensure smoothness in the power law.
};
} // namespace polyMFEMUtils
inline double dfmod(const double epsilon, const double n) {
if (n == 0.0) {
return 0.0;
}
if (n == 1.0) {
return 1.0;
}
return n * std::pow(epsilon, n - 1.0);
}
inline double fmod(const double theta, const double n, const double epsilon) {
if (n == 0.0) {
return 1.0;
}
// For n != 0
const double y_prime_at_epsilon = dfmod(epsilon, n); // Uses the robust dfmod
const double y_at_epsilon = std::pow(epsilon, n); // epsilon^n
// f_mod(theta) = y_at_epsilon + y_prime_at_epsilon * (theta - epsilon)
return y_at_epsilon + y_prime_at_epsilon * (theta - epsilon);
}
} // namespace polyMFEMUtils
} // namespace polytrope
} // namespace serif

399
src/polytrope/utils/public/polytropeOperator.h Normal file
View File

@@ -0,0 +1,399 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: April 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include "mfem.hpp"
#include "4DSTARTypes.h"
#include <memory>
#include "probe.h"
namespace serif {
namespace polytrope {
/**
* @brief Represents the Schur complement operator used in the solution process.
*
* This class computes S = D - Q * GradInv * M.
*/
class SchurCompliment final : public mfem::Operator {
public:
/**
* @brief Constructs a SchurCompliment operator.
* @param QOp The Q matrix operator.
* @param DOp The D matrix operator.
* @param MOp The M matrix operator.
* @param GradInvOp The inverse of the gradient operator.
*/
SchurCompliment(
const mfem::SparseMatrix &QOp,
const mfem::SparseMatrix &DOp,
const mfem::SparseMatrix &MOp,
const mfem::Solver &GradInvOp
);
/**
* @brief Constructs a SchurCompliment operator without the inverse gradient initially.
* The inverse gradient must be set later using updateInverseNonlinearJacobian.
* @param QOp The Q matrix operator.
* @param DOp The D matrix operator.
* @param MOp The M matrix operator.
*/
SchurCompliment(
const mfem::SparseMatrix &QOp,
const mfem::SparseMatrix &DOp,
const mfem::SparseMatrix &MOp
);
/**
* @brief Destructor.
*/
~SchurCompliment() override = default;
/**
* @brief Applies the Schur complement operator: y = (D - Q * GradInv * M) * x.
* @param x The input vector.
* @param y The output vector.
*/
void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
/**
* @brief Sets all operators for the Schur complement.
* @param QOp The Q matrix operator.
* @param DOp The D matrix operator.
* @param MOp The M matrix operator.
* @param GradInvOp The inverse of the gradient operator.
*/
void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp);
/**
* @brief Updates the inverse of the nonlinear Jacobian (GradInv).
* @param gradInv The new inverse gradient solver.
*/
void updateInverseNonlinearJacobian(const mfem::Solver &gradInv);
private:
/**
* @brief Updates the constant matrix terms (Q, D, M).
* @param QOp The Q matrix operator.
* @param DOp The D matrix operator.
* @param MOp The M matrix operator.
*/
void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp);
private:
// Pointers to external operators (not owned by this class)
const mfem::SparseMatrix* m_QOp = nullptr; ///< Pointer to the Q matrix operator.
const mfem::SparseMatrix* m_DOp = nullptr; ///< Pointer to the D matrix operator.
const mfem::SparseMatrix* m_MOp = nullptr; ///< Pointer to the M matrix operator.
const mfem::Solver* m_GradInvOp = nullptr; ///< Pointer to the inverse of the gradient operator.
// Dimensions
int m_nPhi = 0; ///< Dimension related to the phi variable (typically number of rows/cols of D).
int m_nTheta = 0; ///< Dimension related to the theta variable (typically number of rows of M).
// Owned resources
mutable std::unique_ptr<mfem::SparseMatrix> m_matrixForm; ///< Optional: if a matrix representation is ever explicitly formed.
};
/**
* @brief Provides an approximate inverse of the SchurCompliment operator using GMRES.
*/
class GMRESInverter final : public mfem::Operator {
public:
/**
* @brief Constructs a GMRESInverter.
* @param op The SchurCompliment operator to invert.
*/
explicit GMRESInverter(const SchurCompliment& op);
/**
* @brief Destructor.
*/
~GMRESInverter() override = default;
/**
* @brief Applies the GMRES solver to approximate op^-1 * x.
* @param x The input vector.
* @param y The output vector (approximation of op^-1 * x).
*/
void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
private:
const SchurCompliment& m_op; ///< Reference to the SchurCompliment operator to be inverted.
mfem::GMRESSolver m_solver; ///< GMRES solver instance.
};
/**
* @brief Represents the coupled nonlinear operator for the polytropic system.
*
* This operator defines the system F(X) = 0, where X = [θ, φ]^T,
* and F(X) = [ f(θ) + Mφ ]
* [ Dφ - Qθ ].
* It handles essential boundary conditions and the construction of the Jacobian.
*/
class PolytropeOperator final : public mfem::Operator {
public:
/**
* @brief Constructs the PolytropeOperator.
* @param M The M bilinear form (coupling θ and φ).
* @param Q The Q bilinear form (coupling φ and θ).
* @param D The D bilinear form (acting on φ).
* @param f The nonlinear form f(θ).
* @param blockOffsets Offsets defining the blocks for θ and φ variables.
*/
PolytropeOperator(
std::unique_ptr<mfem::MixedBilinearForm> M,
std::unique_ptr<mfem::MixedBilinearForm> Q,
std::unique_ptr<mfem::BilinearForm> D,
std::unique_ptr<mfem::BilinearForm> S,
std::unique_ptr<mfem::NonlinearForm> f,
const mfem::Array<int> &blockOffsets);
/**
* @brief Destructor.
*/
~PolytropeOperator() override = default;
/**
* @brief Applies the PolytropeOperator: y = F(x).
* This computes the residual of the nonlinear system.
* @param xFree The vector of free (non-essential) degrees of freedom.
* @param yFree The resulting residual vector corresponding to free DOFs.
*/
void Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const override;
/**
* @brief Computes the Jacobian of the PolytropeOperator at a given state xFree.
* The Jacobian is of the form:
* J = [ G M ]
* [ -Q D ]
* where G is the gradient of f(θ).
* @param xFree The vector of free DOFs at which to evaluate the gradient.
* @return A reference to the Jacobian operator.
*/
mfem::Operator& GetGradient(const mfem::Vector &xFree) const override;
/**
* @brief Finalizes the operator setup.
* This must be called after setting essential true DOFs and before using Mult or GetGradient.
* It constructs reduced matrices, block offsets, and populates free DOFs.
* @param initTheta Initial guess for θ, used to evaluate the nonlinear form if necessary during setup.
*/
void finalize(const mfem::Vector &initTheta);
/**
* @brief Checks if the operator has been finalized.
* @return True if finalize() has been successfully called, false otherwise.
*/
bool isFinalized() const { return m_isFinalized; }
/**
* @brief Sets the essential true degrees of freedom for both θ and φ variables.
* Marks the operator as not finalized.
* @param theta_ess_tdofs Pair of arrays: (indices of essential DOFs for θ, values at these DOFs).
* @param phi_ess_tdofs Pair of arrays: (indices of essential DOFs for φ, values at these DOFs).
*/
void set_essential_true_dofs(const serif::types::MFEMArrayPair& theta_ess_tdofs, const serif::types::MFEMArrayPair& phi_ess_tdofs);
/**
* @brief Sets the essential true degrees of freedom using a pair of pairs.
* Marks the operator as not finalized.
* @param ess_tdof_pair_set A pair containing the essential TDOF pairs for theta and phi.
*/
void set_essential_true_dofs(const serif::types::MFEMArrayPairSet& ess_tdof_pair_set);
/**
* @brief Reconstructs the full state vector (including essential DOFs) from a reduced state vector (free DOFs).
* @param reducedState The vector containing only the free degrees of freedom.
* @return Constant reference to the internal full state vector, updated with the reducedState.
*/
[[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const;
/**
* @breif Reconstruct the full state vector (including essential DOFs) from a reduced state vector (free DOFs) as well as the block offsets.
* @param reducedState The vector containing only the free degrees of freedom.
* @return Constant reference to the internal full state vector, updated with the reducedState as a block vector.
*/
[[nodiscard]] const mfem::BlockVector reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const;
/**
* @brief Populates the internal array of free (non-essential) degree of freedom indices.
* This is called during finalize().
*/
void populate_free_dof_array();
/// --- Getters for key internal state and operators ---
/**
* @brief Gets the Jacobian operator.
* Asserts that the operator is finalized and the Jacobian has been computed.
* @return Constant reference to the block Jacobian operator.
*/
const mfem::BlockOperator &get_jacobian_operator() const;
/**
* @brief Gets the block diagonal preconditioner for the Schur complement system.
* Asserts that the operator is finalized and the preconditioner has been computed.
* @return Reference to the block diagonal preconditioner.
*/
mfem::BlockDiagonalPreconditioner &get_preconditioner() const;
/// --- Getters for information on internal state ---
/**
* @brief Gets the full state vector, including essential DOFs.
* @return Constant reference to the internal state vector.
*/
const mfem::Array<int>& get_free_dofs() const { return m_freeDofs; } ///< Getter for the free DOFs array.
/**
* @brief Gets the size of the reduced system (number of free DOFs).
* Asserts that the operator is finalized.
* @return The total number of free degrees of freedom.
*/
int get_reduced_system_size() const;
/**
* @brief Gets the currently set essential true degrees of freedom.
* @return A pair containing the essential TDOF pairs for theta and phi.
*/
serif::types::MFEMArrayPairSet get_essential_true_dofs() const;
/**
* @brief Gets the block offsets for the full (unreduced) system.
* @return Constant reference to the array of block offsets.
*/
const mfem::Array<int>& get_block_offsets() const { return m_blockOffsets; }
/**
* @brief Gets the block offsets for the reduced system (after eliminating essential DOFs).
* @return Constant reference to the array of reduced block offsets.
*/
const mfem::Array<int>& get_reduced_block_offsets() const {return m_reducedBlockOffsets; }
private:
// --- Logging ---
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance(); ///< Reference to the global log manager.
quill::Logger* m_logger = m_logManager.getLogger("log"); ///< Pointer to the specific logger instance.
// --- Input Bilinear/Nonlinear Forms (owned) ---
std::unique_ptr<mfem::MixedBilinearForm> m_M; ///< Bilinear form M, coupling θ and φ.
std::unique_ptr<mfem::MixedBilinearForm> m_Q; ///< Bilinear form Q, coupling φ and θ.
std::unique_ptr<mfem::BilinearForm> m_D; ///< Bilinear form D, acting on φ.
std::unique_ptr<mfem::BilinearForm> m_S;
std::unique_ptr<mfem::NonlinearForm> m_f; ///< Nonlinear form f, acting on θ.
// --- Full Matrix Representations (owned, derived from forms) ---
std::unique_ptr<mfem::SparseMatrix> m_Mmat; ///< Sparse matrix representation of M.
std::unique_ptr<mfem::SparseMatrix> m_Qmat; ///< Sparse matrix representation of Q.
std::unique_ptr<mfem::SparseMatrix> m_Dmat; ///< Sparse matrix representation of D.
std::unique_ptr<mfem::SparseMatrix> m_Smat;
// --- Reduced Matrix Representations (owned, after eliminating essential DOFs) ---
std::unique_ptr<mfem::SparseMatrix> m_MReduced; ///< Reduced M matrix (free DOFs only).
std::unique_ptr<mfem::SparseMatrix> m_QReduced; ///< Reduced Q matrix (free DOFs only).
std::unique_ptr<mfem::SparseMatrix> m_DReduced; ///< Reduced D matrix (free DOFs only).
std::unique_ptr<mfem::SparseMatrix> m_SReduced; ///< Reduced S matrix (free DOFs only, used for least squares stabilization).
mutable std::unique_ptr<mfem::SparseMatrix> m_gradReduced; ///< Reduced gradient operator (G) for the nonlinear part f(θ).
// --- Scaled Reduced Matrix Representations (owned, after eliminating essential DOFs and scaling by stabilization coefficients) ---
std::unique_ptr<mfem::SparseMatrix> m_MScaledReduced; ///< Scaled M matrix (free DOFs only, scaled by stabilization coefficient).
std::unique_ptr<mfem::SparseMatrix> m_QScaledReduced; ///< Scaled Q matrix (free DOFs only, scaled by stabilization coefficient).
std::unique_ptr<mfem::SparseMatrix> m_DScaledReduced; ///< Scaled D matrix (free DOFs only, scaled by stabilization coefficient).
std::unique_ptr<mfem::SparseMatrix> m_ScaledSReduced; ///< Scaled S matrix (free DOFs only, scaled by stabilization coefficient).
// --- Stabilization Coefficient --- (Perhapses this should be a user parameter...) // TODO: Sort out why this is negative (sign error?)
static constexpr double m_stabilizationCoefficient = -2.0; ///< Stabilization coefficient for the system, used to more tightly couple ∇θ and φ.
static constexpr double m_IncrementedStabilizationCoefficient = 1.0 + m_stabilizationCoefficient; ///< 1 + Stabilization coefficient for the system, used to more tightly couple ∇θ and φ.
// --- State Vectors and DOF Management ---
mutable mfem::Vector m_state; ///< Full state vector [θ, φ]^T, including essential DOFs.
mfem::Array<int> m_freeDofs; ///< Array of indices for free (non-essential) DOFs.
// --- Block Offsets ---
const mfem::Array<int> m_blockOffsets; ///< Block offsets for the full system [0, size(θ), size(θ)+size(φ)].
mfem::Array<int> m_reducedBlockOffsets; ///< Block offsets for the reduced system (free DOFs).
// --- Essential Boundary Conditions ---
serif::types::MFEMArrayPair m_theta_ess_tdofs; ///< Essential true DOFs for θ (indices and values).
serif::types::MFEMArrayPair m_phi_ess_tdofs; ///< Essential true DOFs for φ (indices and values).
// --- Jacobian and Preconditioner Components (owned, mutable) ---
std::unique_ptr<mfem::ScaledOperator> m_negQ_mat; ///< Scaled operator for -Q_reduced.
mutable std::unique_ptr<mfem::BlockOperator> m_jacobian; ///< Jacobian operator J = [G M; -Q D]_reduced.
mutable std::unique_ptr<SchurCompliment> m_schurCompliment; ///< Schur complement S = D_reduced - Q_reduced * G_inv_reduced * M_reduced.
mutable std::unique_ptr<GMRESInverter> m_invSchurCompliment; ///< Approximate inverse of the Schur complement.
mutable std::unique_ptr<mfem::Solver> m_invNonlinearJacobian; ///< Solver for the inverse of the G block (gradient of f(θ)_reduced).
mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner; ///< Block diagonal preconditioner for the system.
// --- State Flags ---
bool m_isFinalized = false; ///< Flag indicating if finalize() has been called.
private:
/**
* @brief Constructs the sparse matrix representations of M, Q, and D, and their reduced forms.
* Called during finalize().
*/
void construct_matrix_representations();
/**
* @brief Constructs the block offsets for the reduced system.
* Called during finalize().
*/
void construct_reduced_block_offsets();
/**
* @brief Constructs the constant terms of the Jacobian operator (M, -Q, D).
* The (0,0) block (gradient of f) is set in GetGradient.
* Called during finalize().
*/
void construct_jacobian_constant_terms();
/**
* @brief Scatters the values of essential boundary conditions into the full state vector.
* Called during finalize().
*/
void scatter_boundary_conditions();
/**
* @brief Updates the solver for the inverse of the nonlinear Jacobian block (G_00).
* @param grad The gradient operator (G_00) of the nonlinear part f(θ).
*/
void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const;
/**
* @brief Updates the inverse Schur complement operator and its components.
* This is typically called after the nonlinear Jacobian part has been updated.
*/
void update_inverse_schur_compliment() const;
/**
* @brief Updates the preconditioner components.
* This involves updating the inverse nonlinear Jacobian and then the inverse Schur complement.
* @param grad The gradient operator (G_00) of the nonlinear part f(θ).
*/
void update_preconditioner(const mfem::Operator &grad) const;
};
} // namespace polytrope
} // namespace serif

62
src/polytrope/utils/public/utilities.h Normal file
View File

@@ -0,0 +1,62 @@
#pragma once
#include "mfem.hpp"
namespace serif::utilities {
[[nodiscard]] mfem::SparseMatrix build_reduced_matrix(
const mfem::SparseMatrix& matrix,
const mfem::Array<int>& trialEssentialDofs,
const mfem::Array<int>& testEssentialDofs
);
/**
* @brief Generate a vector of 1s and 0s where 1 elemetns cooresponds to queried dofs. Useful for degugging
* @param allDofs array, counding from 0, of all dofs in the system
* @param highlightDofs the dofs that you want to identify
* @return
*
* *Example Usage:*
* One could use this to identify, for example, which dofs are being identified as the central dofs
* @code
* ...
* mfem::Array<int> phiDofs, thetaDofs;
* phiDofs.SetSize(m_fePhi->GetNDofs());
* thetaDofs.SetSize(m_feTheta->GetNDofs());
* const mfem::Vector phiHighlightVector = serif::utilities::build_dof_identification_vector(phiDofs, phiCenterDofs);
* const mfem::Vector thetaHighlightVector = serif::utilities::build_dof_identification_vector(thetaDofs, thetaCenterDofs);
* Probe::glVisView(
* const_cast<mfem::Vector&>(phiHighlightVector),
* *m_fePhi,
* "Phi Center Dofs"
* );
* Probe::glVisView(
* const_cast<mfem::Vector&>(thetaHighlightVector),
* *m_feTheta,
* "Theta Center Dofs"
* );
* @endcode
*/
mfem::Vector build_dof_identification_vector(
const mfem::Array<int>& allDofs,
const::mfem::Array<int>& highlightDofs
);
/**
* @brief Computes the curl of a given H(div) grid function (e.g., from an RT space).
*
* This function is crucial for diagnosing spurious, non-physical modes in mixed FEM
* formulations where the curl of a gradient field is expected to be zero.
*
* @param phi_gf The GridFunction representing the vector field (e.g., φ). It is expected
* to be in an H(div)-conforming space like Raviart-Thomas (RT).
* @return A std::pair containing two new grid functions:
* - pair.first: A unique_ptr to the vector curl field (∇ × φ). This field will
* be in an H(curl)-conforming Nedelec (ND) space.
* - pair.second: A unique_ptr to the scalar magnitude of the curl (||∇ × φ||).
* This field will be in an L2 space.
*
* @note The returned unique_ptrs manage the lifetime of the new GridFunctions and their
* associated FiniteElementSpaces and FECollections, preventing memory leaks.
*/
mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf);
}

27
src/probe/private/probe.cpp
View File

@@ -37,12 +37,12 @@
#include "mfem.hpp"
#include "config.h"
#include "probe.h"
#include "probe.h"
#include "warning_control.h"
namespace Probe {
namespace serif::probe {
void pause() {
std::cout << "Execution paused. Please press enter to continue..."
@@ -56,22 +56,20 @@ void wait(int seconds) {
void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
const std::string& windowTitle, const std::string& keyset) {
Config& config = Config::getInstance();
serif::config::Config& config = serif::config::Config::getInstance();
quill::Logger* logger = LogManager::getInstance().getLogger("log");
std::string usedKeyset;
if (config.get<bool>("Probe:GLVis:Visualization", true)) {
std::string usedKeyset;
LOG_INFO(logger, "Visualizing solution using GLVis...");
LOG_INFO(logger, "Window title: {}", windowTitle);
if (keyset == "") {
if (keyset.empty()) {
usedKeyset = config.get<std::string>("Probe:GLVis:DefaultKeyset", "");
} else {
usedKeyset = keyset;
}
LOG_INFO(logger, "Keyset: {}", usedKeyset);
std::string vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
int visport = config.get<int>("Probe:GLVis:Port", 19916);
std::cout << "GLVis visualization enabled. Opening GLVis window... " << std::endl;
std::cout << "Using host: " << vishost << " and port: " << visport << std::endl;
const auto vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
const auto visport = config.get<int>("Probe:GLVis:Port", 19916);
mfem::socketstream sol_sock(vishost.c_str(), visport);
sol_sock.precision(8);
sol_sock << "solution\n" << mesh << u
@@ -83,7 +81,7 @@ void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle, const std::string& keyset) {
mfem::GridFunction gf(&fes);
DEPRECATION_WARNING_OFF
gf.SetData(vec);
DEPRECATION_WARNING_ON
@@ -110,8 +108,8 @@ double getMeshRadius(mfem::Mesh& mesh) {
std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
const std::vector<double>& rayDirection,
int numSamples, std::string filename) {
Config& config = Config::getInstance();
Probe::LogManager& logManager = Probe::LogManager::getInstance();
serif::config::Config& config = serif::config::Config::getInstance();
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
LOG_INFO(logger, "Getting ray solution...");
// Check if the directory to write to exists
@@ -172,7 +170,8 @@ std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::GridFun
LOG_DEBUG(logger, "Probe::getRaySolution() : Ray point {} found in element {} with r={:0.2f} and theta={:0.2f}", i, elementId, r, sampleValue);
samples.push_back(sampleValue);
} else { // If the point was not found in an element
samples.push_back(0.0);
LOG_INFO(logger, "Probe::getRaySolution() : Ray point {} not found", i);
samples.push_back(0.0);
}
}
std::pair<std::vector<double>, std::vector<double>> samplesPair(radialPoints, samples);
@@ -205,7 +204,7 @@ std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::Vector
}
LogManager::LogManager() {
Config& config = Config::getInstance();
serif::config::Config& config = serif::config::Config::getInstance();
quill::Backend::start();
auto CLILogger = quill::Frontend::create_or_get_logger(
"root",

14
src/probe/public/probe.h
View File

@@ -19,8 +19,7 @@
//
// *********************************************************************** */
//=== Probe.h ===
#ifndef PROBE_H
#define PROBE_H
#pragma once
#include <string>
#include <map>
@@ -33,7 +32,7 @@
/**
* @brief The Probe namespace contains utility functions for debugging and logging.
*/
namespace Probe {
namespace serif::probe {
/**
* @brief Pause the execution and wait for user input.
*/
@@ -54,7 +53,7 @@ namespace Probe {
*/
void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
const std::string& windowTitle = "grid function", const std::string& keyset="");
/**
* @brief Visualize a vector using GLVis.
* @param vec The vector to visualize.
@@ -63,8 +62,8 @@ namespace Probe {
* @param keyset The keyset to use for visualization.
*/
void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes,
const std::string &windowTitle = "vector", const std::string& keyset="");
const std::string &windowTitle = "vector", const std::string& keyset="");
double getMeshRadius(mfem::Mesh& mesh);
std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
@@ -135,5 +134,4 @@ namespace Probe {
const std::string& loggerName);
};
} // namespace Probe
#endif
} // namespace Probe

4
src/python/bindings.cpp
View File

@@ -7,8 +7,8 @@
#include "config/bindings.h"
#include "eos/bindings.h"
PYBIND11_MODULE(fourdsse_bindings, m) {
m.doc() = "Python bindings for the 4DSSE project";
PYBIND11_MODULE(serif, m) {
m.doc() = "Python bindings for the SERiF project";
auto compMod = m.def_submodule("composition", "Composition-module bindings");
register_comp_bindings(compMod);

70
src/python/composition/bindings.cpp
View File

@@ -15,7 +15,7 @@ std::string sv_to_string(std::string_view sv) {
return std::string(sv);
}
std::string get_ostream_str(const composition::Composition& comp) {
std::string get_ostream_str(const serif::composition::Composition& comp) {
std::ostringstream oss;
oss << comp;
return oss.str();
@@ -24,38 +24,38 @@ std::string get_ostream_str(const composition::Composition& comp) {
void register_comp_bindings(pybind11::module &comp_submodule) {
// --- Bindings for composition and species module ---
py::class_<composition::GlobalComposition>(comp_submodule, "GlobalComposition")
.def_readonly("specificNumberDensity", &composition::GlobalComposition::specificNumberDensity)
.def_readonly("meanParticleMass", &composition::GlobalComposition::meanParticleMass)
py::class_<serif::composition::GlobalComposition>(comp_submodule, "GlobalComposition")
.def_readonly("specificNumberDensity", &serif::composition::GlobalComposition::specificNumberDensity)
.def_readonly("meanParticleMass", &serif::composition::GlobalComposition::meanParticleMass)
.def("__repr__", // Add a string representation for easy printing in Python
[](const composition::GlobalComposition &gc) {
[](const serif::composition::GlobalComposition &gc) {
return "<GlobalComposition(specNumDens=" + std::to_string(gc.specificNumberDensity) +
", meanMass=" + std::to_string(gc.meanParticleMass) + ")>";
});
py::class_<composition::CompositionEntry>(comp_submodule, "CompositionEntry")
.def("symbol", &composition::CompositionEntry::symbol)
py::class_<serif::composition::CompositionEntry>(comp_submodule, "CompositionEntry")
.def("symbol", &serif::composition::CompositionEntry::symbol)
.def("mass_fraction",
py::overload_cast<>(&composition::CompositionEntry::mass_fraction, py::const_),
py::overload_cast<>(&serif::composition::CompositionEntry::mass_fraction, py::const_),
"Gets the mass fraction of the species.")
.def("mass_fraction",
py::overload_cast<double>(&composition::CompositionEntry::mass_fraction, py::const_),
py::overload_cast<double>(&serif::composition::CompositionEntry::mass_fraction, py::const_),
py::arg("meanMolarMass"), // Name the argument in Python
"Gets the mass fraction of the species given the mean molar mass.")
.def("number_fraction",
py::overload_cast<>(&composition::CompositionEntry::number_fraction, py::const_),
py::overload_cast<>(&serif::composition::CompositionEntry::number_fraction, py::const_),
"Gets the number fraction of the species.")
.def("number_fraction",
py::overload_cast<double>(&composition::CompositionEntry::number_fraction, py::const_),
py::overload_cast<double>(&serif::composition::CompositionEntry::number_fraction, py::const_),
py::arg("totalMoles"),
"Gets the number fraction of the species given the total moles.")
.def("rel_abundance", &composition::CompositionEntry::rel_abundance)
.def("isotope", &composition::CompositionEntry::isotope) // Assuming Species is bound or convertible
.def("getMassFracMode", &composition::CompositionEntry::getMassFracMode)
.def("rel_abundance", &serif::composition::CompositionEntry::rel_abundance)
.def("isotope", &serif::composition::CompositionEntry::isotope) // Assuming Species is bound or convertible
.def("getMassFracMode", &serif::composition::CompositionEntry::getMassFracMode)
.def("__repr__", // Optional: nice string representation
[](const composition::CompositionEntry &ce) {
[](const serif::composition::CompositionEntry &ce) {
// You might want to include more info here now
return "<CompositionEntry(symbol='" + ce.symbol() + "', " +
"mass_frac=" + std::to_string(ce.mass_fraction()) + ", " +
@@ -63,7 +63,7 @@ void register_comp_bindings(pybind11::module &comp_submodule) {
});
// --- Binding for the main Composition class ---
py::class_<composition::Composition>(comp_submodule, "Composition")
py::class_<serif::composition::Composition>(comp_submodule, "Composition")
// Constructors
.def(py::init<>(), "Default constructor")
.def(py::init<const std::vector<std::string>&>(),
@@ -75,61 +75,61 @@ void register_comp_bindings(pybind11::module &comp_submodule) {
"Constructor taking symbols, fractions, and mode (True=Mass, False=Number)")
// Methods
.def("finalize", &composition::Composition::finalize, py::arg("norm") = false,
.def("finalize", &serif::composition::Composition::finalize, py::arg("norm") = false,
"Finalize the composition, optionally normalizing fractions to sum to 1.")
.def("registerSymbol", py::overload_cast<const std::string&, bool>(&composition::Composition::registerSymbol),
.def("registerSymbol", py::overload_cast<const std::string&, bool>(&serif::composition::Composition::registerSymbol),
py::arg("symbol"), py::arg("massFracMode") = true, "Register a single symbol.")
.def("registerSymbol", py::overload_cast<const std::vector<std::string>&, bool>(&composition::Composition::registerSymbol),
.def("registerSymbol", py::overload_cast<const std::vector<std::string>&, bool>(&serif::composition::Composition::registerSymbol),
py::arg("symbols"), py::arg("massFracMode") = true, "Register multiple symbols.")
.def("getRegisteredSymbols", &composition::Composition::getRegisteredSymbols,
.def("getRegisteredSymbols", &serif::composition::Composition::getRegisteredSymbols,
"Get the set of registered symbols.")
.def("setMassFraction", py::overload_cast<const std::string&, const double&>(&composition::Composition::setMassFraction),
.def("setMassFraction", py::overload_cast<const std::string&, const double&>(&serif::composition::Composition::setMassFraction),
py::arg("symbol"), py::arg("mass_fraction"), "Set mass fraction for a single symbol (requires massFracMode). Returns old value.")
.def("setMassFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&composition::Composition::setMassFraction),
.def("setMassFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&serif::composition::Composition::setMassFraction),
py::arg("symbols"), py::arg("mass_fractions"), "Set mass fractions for multiple symbols (requires massFracMode). Returns list of old values.")
.def("setNumberFraction", py::overload_cast<const std::string&, const double&>(&composition::Composition::setNumberFraction),
.def("setNumberFraction", py::overload_cast<const std::string&, const double&>(&serif::composition::Composition::setNumberFraction),
py::arg("symbol"), py::arg("number_fraction"), "Set number fraction for a single symbol (requires !massFracMode). Returns old value.")
.def("setNumberFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&composition::Composition::setNumberFraction),
.def("setNumberFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&serif::composition::Composition::setNumberFraction),
py::arg("symbols"), py::arg("number_fractions"), "Set number fractions for multiple symbols (requires !massFracMode). Returns list of old values.")
.def("mix", &composition::Composition::mix, py::arg("other"), py::arg("fraction"),
.def("mix", &serif::composition::Composition::mix, py::arg("other"), py::arg("fraction"),
"Mix with another composition. Returns new Composition.")
.def("getMassFraction", py::overload_cast<const std::string&>(&composition::Composition::getMassFraction, py::const_),
.def("getMassFraction", py::overload_cast<const std::string&>(&serif::composition::Composition::getMassFraction, py::const_),
py::arg("symbol"), "Get mass fraction for a symbol (calculates if needed). Requires finalization.")
.def("getMassFraction", py::overload_cast<>(&composition::Composition::getMassFraction, py::const_),
.def("getMassFraction", py::overload_cast<>(&serif::composition::Composition::getMassFraction, py::const_),
"Get dictionary of all mass fractions. Requires finalization.")
.def("getNumberFraction", py::overload_cast<const std::string&>(&composition::Composition::getNumberFraction, py::const_),
.def("getNumberFraction", py::overload_cast<const std::string&>(&serif::composition::Composition::getNumberFraction, py::const_),
py::arg("symbol"), "Get number fraction for a symbol (calculates if needed). Requires finalization.")
.def("getNumberFraction", py::overload_cast<>(&composition::Composition::getNumberFraction, py::const_),
.def("getNumberFraction", py::overload_cast<>(&serif::composition::Composition::getNumberFraction, py::const_),
"Get dictionary of all number fractions. Requires finalization.")
// Note: pybind11 automatically converts std::pair to a Python tuple
.def("getComposition", py::overload_cast<const std::string&>(&composition::Composition::getComposition, py::const_),
.def("getComposition", py::overload_cast<const std::string&>(&serif::composition::Composition::getComposition, py::const_),
py::arg("symbol"), "Returns a tuple (CompositionEntry, GlobalComposition) for the symbol. Requires finalization.")
// Binding the version returning map<string, Entry> requires a bit more care or helper function
// to convert the map to a Python dict if needed directly. Let's bind the pair version for now.
.def("getComposition", py::overload_cast<>(&composition::Composition::getComposition, py::const_),
.def("getComposition", py::overload_cast<>(&serif::composition::Composition::getComposition, py::const_),
"Returns a tuple (dict[str, CompositionEntry], GlobalComposition) for all symbols. Requires finalization.")
.def("subset", &composition::Composition::subset, py::arg("symbols"), py::arg("method") = "norm",
.def("subset", &serif::composition::Composition::subset, py::arg("symbols"), py::arg("method") = "norm",
"Create a new Composition containing only the specified symbols.")
.def("hasSymbol", &composition::Composition::hasSymbol, py::arg("symbol"),
.def("hasSymbol", &serif::composition::Composition::hasSymbol, py::arg("symbol"),
"Check if a symbol is registered.")
.def("setCompositionMode", &composition::Composition::setCompositionMode, py::arg("massFracMode"),
.def("setCompositionMode", &serif::composition::Composition::setCompositionMode, py::arg("massFracMode"),
"Set the mode (True=Mass, False=Number). Requires finalization before switching.")
// Operator overload
.def(py::self + py::self, "Mix equally with another composition.") // Binds operator+
// Add __repr__ or __str__
.def("__repr__", [](const composition::Composition &comp) {
.def("__repr__", [](const serif::composition::Composition &comp) {
return get_ostream_str(comp); // Use helper for C++ operator<<
});

10
src/python/config/bindings.cpp
View File

@@ -14,7 +14,7 @@ template <typename T>
void def_config_get(py::module &m) {
m.def("get",
[](const std::string &key, T defaultValue) {
return Config::getInstance().get<T>(key, defaultValue);
return serif::config::Config::getInstance().get<T>(key, defaultValue);
},
py::arg("key"), py::arg("defaultValue"),
"Get configuration value (type inferred from default)");
@@ -28,28 +28,28 @@ void register_config_bindings(pybind11::module &config_submodule) {
config_submodule.def("loadConfig",
[](const std::string& configFilePath) {
return Config::getInstance().loadConfig(configFilePath);
return serif::config::Config::getInstance().loadConfig(configFilePath);
},
py::arg("configFilePath"),
"Load configuration from a YAML file.");
config_submodule.def("has",
[](const std::string &key) {
return Config::getInstance().has(key);
return serif::config::Config::getInstance().has(key);
},
py::arg("key"),
"Check if a key exists in the configuration.");
config_submodule.def("keys",
[]() {
return py::cast(Config::getInstance().keys());
return py::cast(serif::config::Config::getInstance().keys());
},
"Get a list of all configuration keys.");
config_submodule.def("__repr__",
[]() {
std::ostringstream oss;
oss << Config::getInstance(); // Use the existing operator<<
oss << serif::config::Config::getInstance(); // Use the existing operator<<
return std::string("<fourdsse_bindings.config module accessing C++ Singleton>\n") + oss.str();
});
}

26
src/python/const/bindings.cpp
View File

@@ -10,44 +10,44 @@ namespace py = pybind11;
void register_const_bindings(pybind11::module &const_submodule) {
py::class_<Constant>(const_submodule, "Constant")
.def_readonly("name", &Constant::name)
.def_readonly("value", &Constant::value)
.def_readonly("uncertainty", &Constant::uncertainty)
.def_readonly("unit", &Constant::unit)
.def_readonly("reference", &Constant::reference)
.def("__repr__", [](const Constant &c) {
py::class_<serif::constant::Constant>(const_submodule, "Constant")
.def_readonly("name", &serif::constant::Constant::name)
.def_readonly("value", &serif::constant::Constant::value)
.def_readonly("uncertainty", &serif::constant::Constant::uncertainty)
.def_readonly("unit", &serif::constant::Constant::unit)
.def_readonly("reference", &serif::constant::Constant::reference)
.def("__repr__", [](const serif::constant::Constant &c) {
return "<Constant(name='" + c.name + "', value=" + std::to_string(c.value) +
", uncertainty=" + std::to_string(c.uncertainty) +
", unit='" + c.unit + "')>";
});
py::class_<Constants>(const_submodule, "Constants")
.def_property_readonly("loaded", &Constants::isLoaded)
py::class_<serif::constant::Constants>(const_submodule, "Constants")
.def_property_readonly("loaded", &serif::constant::Constants::isLoaded)
.def_static("get",
[](const std::string &name) {
return py::cast(
Constants::getInstance().get(name)
serif::constant::Constants::getInstance().get(name)
);
},
"Get a constant by name. Returns None if not found."
)
.def_static("has",
[](const std::string &name) {
return Constants::getInstance().has(name);
return serif::constant::Constants::getInstance().has(name);
},
"Check if a constant exists by name.")
.def_static("keys",
[]() {
return py::cast(
Constants::getInstance().keys()
serif::constant::Constants::getInstance().keys()
);
},
"Get a list of all constant names.")
.def_static("__class_getitem__",
[](const std::string &name) {
return py::cast(
Constants::getInstance().get(name)
serif::constant::Constants::getInstance().get(name)
);
});

83
src/resource/private/resourceManager.cpp
View File

@@ -34,55 +34,56 @@
#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)
ResourceManager::ResourceManager() {
std::string defaultDataDir = TOSTRING(DATA_DIR);
m_dataDir = m_config.get<std::string>("Data:Dir", defaultDataDir);
// -- Get the index file path using filesytem to make it a system safe path
std::string indexFilePath = m_dataDir + "/index.yaml";
std::filesystem::path indexFile(indexFilePath);
namespace serif::resource {
ResourceManager::ResourceManager() {
const std::string defaultDataDir = TOSTRING(DATA_DIR);
m_dataDir = m_config.get<std::string>("Data:Dir", defaultDataDir);
// -- Get the index file path using filesystem to make it a system safe path
const std::string indexFilePath = m_dataDir + "/index.yaml";
// TODO Add checks to make sure data dir exists and index.yaml exists
const std::filesystem::path indexFile(indexFilePath);
m_resourceConfig.loadConfig(indexFile.string());
std::vector<std::string> assets = m_resourceConfig.keys();
for (auto key : assets ) {
load(key);
m_resourceConfig.loadConfig(indexFile.string());
for (const auto& key : m_resourceConfig.keys() ) {
load(key);
}
}
}
std::vector<std::string> ResourceManager::getAvaliableResources() {
std::vector<std::string> resources;
resources = m_resourceConfig.keys();
return resources;
}
const Resource& ResourceManager::getResource(const std::string &name) const {
auto it = m_resources.find(name);
if (it != m_resources.end()) {
return it->second;
std::vector<std::string> ResourceManager::getAvailableResources() const {
const std::vector<std::string> resources = m_resourceConfig.keys();
return resources;
}
throw std::runtime_error("Resource " + name + " not found");
}
bool ResourceManager::loadResource(std::string& name) {
return load(name);
}
bool ResourceManager::load(const std::string& name) {
const std::string resourcePath = m_dataDir + "/" + m_resourceConfig.get<std::string>(name);
std::filesystem::path resourceFile(resourcePath);
if (!std::filesystem::exists(resourceFile)) {
LOG_ERROR(m_logger, "Resource file not found: {}", resourceFile.string());
return false;
const serif::resource::types::Resource& ResourceManager::getResource(const std::string &name) const {
if (const auto it = m_resources.find(name); it != m_resources.end()) {
return it->second;
}
throw std::runtime_error("Resource " + name + " not found");
}
LOG_INFO(m_logger, "Loading resource: {}", resourceFile.string());
if (m_resources.find(name) != m_resources.end()) {
LOG_INFO(m_logger, "Resource already loaded: {}", name);
bool ResourceManager::loadResource(std::string& name) {
return load(name);
}
bool ResourceManager::load(const std::string& name) {
// TODO : make this path system more system agnostic
const std::string resourcePath = m_dataDir + "/" + m_resourceConfig.get<std::string>(name);
const std::filesystem::path resourceFile(resourcePath);
if (!std::filesystem::exists(resourceFile)) {
LOG_ERROR(m_logger, "Resource file not found: {}", resourceFile.string());
return false;
}
LOG_INFO(m_logger, "Loading resource: {}", resourceFile.string());
if (m_resources.contains(name)) {
LOG_INFO(m_logger, "Resource already loaded: {}", name);
return true;
}
serif::resource::types::Resource resource = serif::resource::types::createResource(name, resourcePath);
m_resources[name] = std::move(resource);
// -- Check if the resource is already in the map
return true;
}
Resource resource = createResource(name, resourcePath);
m_resources[name] = std::move(resource);
// -- Check if the resource is already in the map
return true;
}

56
src/resource/private/resourceManagerTypes.cpp
View File

@@ -27,35 +27,37 @@
#include "debug.h"
std::string getFirstSegment(const std::string& input) {
size_t pos = input.find(':');
if (pos == std::string::npos) {
// No colon found, return the entire string
return input;
} else {
// Return substring from start to the position of the first colon
return input.substr(0, pos);
namespace serif::resource::types {
std::string getFirstSegment(const std::string& input) {
size_t pos = input.find(':');
if (pos == std::string::npos) {
// No colon found, return the entire string
return input;
} else {
// Return substring from start to the position of the first colon
return input.substr(0, pos);
}
}
}
Resource createResource(const std::string& type, const std::string& path) {
static const std::unordered_map<std::string, std::function<Resource(const std::string&)>> factoryMap = {
{"opac", [](const std::string& p) { return Resource(
std::make_unique<opat::OPAT>(opat::readOPAT(p)));
}},
{"mesh", [](const std::string& p) { return Resource(
std::make_unique<MeshIO>(p));
}},
{"eos", [](const std::string& p) { return Resource(
std::make_unique<EOSio>(p));
}}
// Add more mappings as needed
};
auto it = factoryMap.find(getFirstSegment(type));
if (it != factoryMap.end()) {
return it->second(path);
} else {
throw std::invalid_argument("Unknown resource type: " + type);
Resource createResource(const std::string& type, const std::string& path) {
static const std::unordered_map<std::string, std::function<Resource(const std::string&)>> factoryMap = {
{"opac", [](const std::string& p) { return Resource(
std::make_unique<opat::OPAT>(opat::readOPAT(p)));
}},
{"mesh", [](const std::string& p) { return Resource(
std::make_unique<serif::mesh::MeshIO>(p));
}},
{"eos", [](const std::string& p) { return Resource(
std::make_unique<serif::eos::EOSio>(p));
}}
// Add more mappings as needed
};
auto it = factoryMap.find(getFirstSegment(type));
if (it != factoryMap.end()) {
return it->second(path);
} else {
throw std::invalid_argument("Unknown resource type: " + type);
}
}
}

172
src/resource/public/resourceManager.h
View File

@@ -22,11 +22,13 @@
#include <vector>
#include <string>
#include <stdexcept>
#include <unordered_map>
#include "resourceManagerTypes.h"
#include "config.h"
#include "probe.h"
#include "quill/LogMacros.h"
/**
* @class ResourceManager
@@ -35,96 +37,98 @@
* The ResourceManager class is responsible for loading, storing, and providing access to resources.
* It follows the Singleton design pattern to ensure only one instance of the manager exists.
*/
class ResourceManager {
private:
/**
* @brief Private constructor to prevent instantiation.
*/
ResourceManager();
namespace serif::resource {
class ResourceManager {
private:
/**
* @brief Private constructor to prevent instantiation.
*/
ResourceManager();
/**
* @brief Deleted copy constructor to prevent copying.
*/
ResourceManager(const ResourceManager&) = delete;
/**
* @brief Deleted copy constructor to prevent copying.
*/
ResourceManager(const ResourceManager&) = delete;
/**
* @brief Deleted assignment operator to prevent assignment.
*/
ResourceManager& operator=(const ResourceManager&) = delete;
/**
* @brief Deleted assignment operator to prevent assignment.
*/
ResourceManager& operator=(const ResourceManager&) = delete;
Config& m_config = Config::getInstance();
Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
Config m_resourceConfig;
std::string m_dataDir;
std::unordered_map<std::string, Resource> m_resources;
serif::config::Config m_resourceConfig;
std::string m_dataDir;
std::unordered_map<std::string, types::Resource> m_resources;
/**
* @brief Loads a resource by name.
* @param name The name of the resource to load.
* @return True if the resource was loaded successfully, false otherwise.
*/
bool load(const std::string& name);
/**
* @brief Loads a resource by name.
* @param name The name of the resource to load.
* @return True if the resource was loaded successfully, false otherwise.
*/
bool load(const std::string& name);
public:
/**
* @brief Gets the singleton instance of the ResourceManager.
* @return The singleton instance of the ResourceManager.
*/
static ResourceManager& getInstance() {
static ResourceManager instance;
return instance;
}
public:
/**
* @brief Gets the singleton instance of the ResourceManager.
* @return The singleton instance of the ResourceManager.
*/
static ResourceManager& getInstance() {
static ResourceManager instance;
return instance;
}
/**
* @brief Gets a list of available resources.
* @return A vector of strings containing the names of available resources.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* std::vector<std::string> resources = manager.getAvaliableResources();
* @endcode
*/
std::vector<std::string> getAvaliableResources();
/**
* @brief Gets a list of available resources.
* @return A vector of strings containing the names of available resources.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* std::vector<std::string> resources = manager.getAvailableResources();
* @endcode
*/
std::vector<std::string> getAvailableResources() const;
/**
* @brief Gets a resource by name.
* @param name The name of the resource to retrieve.
* @return A constant reference to the requested resource.
* @throws std::runtime_error if the resource is not found.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* const Resource& resource = manager.getResource("exampleResource");
* @endcode
*/
const Resource& getResource(const std::string &name) const;
/**
* @brief Loads a resource by name.
* @param name The name of the resource to load.
* @return True if the resource was loaded successfully, false otherwise.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* bool success = manager.loadResource("exampleResource");
* @endcode
*/
bool loadResource(std::string& name);
/**
* @brief Gets a resource by name.
* @param name The name of the resource to retrieve.
* @return A constant reference to the requested resource.
* @throws std::runtime_error if the resource is not found.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* const Resource& resource = manager.getResource("exampleResource");
* @endcode
*/
const types::Resource& getResource(const std::string &name) const;
/**
* @brief Loads all resources.
* @return An unordered map with resource names as keys and load success as values.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* std::unordered_map<std::string, bool> results = manager.loadAllResources();
* @endcode
*/
std::unordered_map<std::string, bool> loadAllResources();
};
/**
* @brief Loads a resource by name.
* @param name The name of the resource to load.
* @return True if the resource was loaded successfully, false otherwise.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* bool success = manager.loadResource("exampleResource");
* @endcode
*/
bool loadResource(std::string& name);
/**
* @brief Loads all resources.
* @return An unordered map with resource names as keys and load success as values.
*
* Example usage:
* @code
* ResourceManager& manager = ResourceManager::getInstance();
* std::unordered_map<std::string, bool> results = manager.loadAllResources();
* @endcode
*/
std::unordered_map<std::string, bool> loadAllResources();
};
}

116
src/resource/public/resourceManagerTypes.h
View File

@@ -18,8 +18,6 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef RESOURCE_MANAGER_TYPES_H
#define RESOURCE_MANAGER_TYPES_H
#include <memory>
#include <variant>
@@ -29,63 +27,63 @@
#include "meshIO.h"
#include "EOSio.h"
/**
* @file resourceManagerTypes.h
* @brief Defines types and functions for managing resources.
*
* This file provides type definitions and functions for handling different
* types of resources in a unified manner.
*/
namespace serif::resource::types {
/**
* @file resourceManagerTypes.h
* @brief Defines types and functions for managing resources.
*
* This file provides type definitions and functions for handling different
* types of resources in a unified manner.
*/
// -- Valid resource types
/**
* @brief A variant type that can hold different types of resources.
*
* The Resource type is a std::variant that can hold a unique pointer to
* an OpatIO, MeshIO, or EOSio object.
*
* Example usage:
* @code
* Resource resource = std::make_unique<OpatIO>(...);
* @endcode
*/
using Resource = std::variant<
std::unique_ptr<opat::OPAT>,
std::unique_ptr<MeshIO>,
std::unique_ptr<EOSio>>;
// -- Valid resource types
/**
* @brief A variant type that can hold different types of resources.
*
* The Resource type is a std::variant that can hold a unique pointer to
* an OpatIO, MeshIO, or EOSio object.
*
* Example usage:
* @code
* Resource resource = std::make_unique<OpatIO>(...);
* @endcode
*/
using Resource = std::variant<
std::unique_ptr<opat::OPAT>,
std::unique_ptr<serif::mesh::MeshIO>,
std::unique_ptr<serif::eos::EOSio>>;
/**
* @brief Extracts the first segment of a given string.
*
* This function takes a string input and returns the first segment
* separated by a delimiter (default is '/').
*
* @param input The input string to be processed.
* @return The first segment of the input string.
*
* Example usage:
* @code
* std::string segment = getFirstSegment("path/to/resource");
* // segment == "path"
* @endcode
*/
std::string getFirstSegment(const std::string& input);
/**
* @brief Extracts the first segment of a given string.
*
* This function takes a string input and returns the first segment
* separated by a delimiter (default is '/').
*
* @param input The input string to be processed.
* @return The first segment of the input string.
*
* Example usage:
* @code
* std::string segment = getFirstSegment("path/to/resource");
* // segment == "path"
* @endcode
*/
std::string getFirstSegment(const std::string& input);
/**
* @brief Creates a resource based on the specified type and path.
*
* This function creates a resource object based on the provided type
* and initializes it using the given path.
*
* @param type The type of the resource to be created (e.g., "OpatIO", "MeshIO", "EOSio").
* @param path The path to initialize the resource with.
* @return A Resource object initialized with the specified type and path.
*
* Example usage:
* @code
* Resource resource = createResource("OpatIO", "path/to/opat");
* @endcode
*/
Resource createResource(const std::string& type, const std::string& path);
#endif // RESOURCE_MANAGER_TYPES_H
/**
* @brief Creates a resource based on the specified type and path.
*
* This function creates a resource object based on the provided type
* and initializes it using the given path.
*
* @param type The type of the resource to be created (e.g., "OpatIO", "MeshIO", "EosIO").
* @param path The path to initialize the resource with.
* @return A Resource object initialized with the specified type and path.
*
* Example usage:
* @code
* Resource resource = createResource("OpatIO", "path/to/opat");
* @endcode
*/
Resource createResource(const std::string& type, const std::string& path);
}

8
src/types/public/4DSTARTypes.h
View File

@@ -18,16 +18,12 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#ifndef _4DSTAR_TYPES_H
#define _4DSTAR_TYPES_H
#pragma once
#include <utility>
#include "mfem.hpp"
// TODO : Need a better namespace name for these types
namespace SSE {
namespace serif::types {
typedef std::pair<mfem::Array<int>, mfem::Array<double>> MFEMArrayPair;
typedef std::pair<MFEMArrayPair, MFEMArrayPair> MFEMArrayPairSet;
}
#endif // _4DSTAR_TYPES_H