From 4f1c26044451dd4840735e10ecdc319948a70d05 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 29 Sep 2025 13:35:48 -0400 Subject: [PATCH] feat(trigger): added working robust repartitioning trigger system more work is needed to identify the most robust set of criteria to trigger on but the system is now very easy to exend, probe, and use. --- meson.build | 2 +- .../gridfire/engine/views/engine_multiscale.h | 14 +- .../triggers/engine_partitioning_trigger.h | 94 ++++ .../trigger/procedures/trigger_pprint.h | 17 + .../gridfire/trigger/trigger_abstract.h | 26 ++ .../gridfire/trigger/trigger_logical.h | 440 ++++++++++++++++++ src/include/gridfire/trigger/trigger_result.h | 13 + src/lib/engine/views/engine_multiscale.cpp | 130 +----- .../strategies/CVODE_solver_strategy.cpp | 173 ++++--- .../triggers/engine_partitioning_trigger.cpp | 263 +++++++++++ src/meson.build | 3 +- tests/graphnet_sandbox/main.cpp | 2 +- 12 files changed, 980 insertions(+), 197 deletions(-) create mode 100644 src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h create mode 100644 src/include/gridfire/trigger/procedures/trigger_pprint.h create mode 100644 src/include/gridfire/trigger/trigger_abstract.h create mode 100644 src/include/gridfire/trigger/trigger_logical.h create mode 100644 src/include/gridfire/trigger/trigger_result.h create mode 100644 src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp diff --git a/meson.build b/meson.build index 2dae91f6..ba0ed33d 100644 --- a/meson.build +++ b/meson.build @@ -81,7 +81,7 @@ subdir('src') if get_option('build-python') message('Configuring Python bindings...') - subdir('src-pybind') + subdir('build-python') else message('Skipping Python bindings...') endif diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index cf4ff91c..04179302 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -57,9 +57,9 @@ namespace gridfire { // TODO: We should probably sort out how to adjust these from absolute to relative tolerances. QSECacheConfig m_cacheConfig = { - 1e-3, // Default tolerance for T9 - 1e-1, // Default tolerance for rho - 1e-3 // Default tolerance for species abundances + 1e-10, // Default tolerance for T9 + 1e-10, // Default tolerance for rho + 1e-10 // Default tolerance for species abundances }; /** @@ -763,7 +763,7 @@ namespace gridfire { std::string toString() const; - std::string toString(DynamicEngine &engine) const; + std::string toString(const DynamicEngine &engine) const; }; /** @@ -989,6 +989,11 @@ namespace gridfire { * @brief Species that are treated as algebraic (in QSE) in the QSE groups. */ std::vector m_algebraic_species; + + /** + * @breif Stateful storage of the current algebraic species abundances. This is updated every time the update method is called. + */ + std::vector m_Y_algebraic; /** * @brief Indices of algebraic species in the full network. */ @@ -1003,7 +1008,6 @@ namespace gridfire { */ std::vector m_activeReactionIndices; - // TODO: Enhance the hashing for the cache to consider not just T and rho but also the current abundance in some careful way that automatically ignores small changes (i.e. network should only be repartitioned sometimes) /** * @brief Cache for QSE abundances based on T9, rho, and Y. * diff --git a/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h b/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h new file mode 100644 index 00000000..eb0d8fdb --- /dev/null +++ b/src/include/gridfire/solver/strategies/triggers/engine_partitioning_trigger.h @@ -0,0 +1,94 @@ +#pragma once + +#include "gridfire/trigger/trigger_abstract.h" +#include "gridfire/trigger/trigger_result.h" +#include "gridfire/solver/strategies/CVODE_solver_strategy.h" +#include "fourdst/logging/logging.h" + +#include +#include +#include + +namespace gridfire::trigger::solver::CVODE { + class SimulationTimeTrigger final : public Trigger { + public: + explicit SimulationTimeTrigger(double interval); + bool check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + void update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) override; + void reset() override; + + std::string name() const override; + TriggerResult why(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + private: + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + + double m_interval; + + mutable double m_last_trigger_time = 0.0; + mutable double m_last_trigger_time_delta = 0.0; + }; + + class OffDiagonalTrigger final : public Trigger { + public: + explicit OffDiagonalTrigger(double threshold); + bool check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + void update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) override; + void reset() override; + + std::string name() const override; + TriggerResult why(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + private: + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + + double m_threshold; + }; + + class TimestepCollapseTrigger final : public Trigger { + public: + explicit TimestepCollapseTrigger(double threshold, bool relative); + explicit TimestepCollapseTrigger(double threshold, bool relative, size_t windowSize); + bool check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + void update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) override; + void reset() override; + + std::string name() const override; + TriggerResult why(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + private: + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + + double m_threshold; + bool m_relative; + size_t m_windowSize; + + std::deque m_timestep_window; + }; + + std::unique_ptr> makeEnginePartitioningTrigger( + const double simulationTimeInterval, + const double offDiagonalThreshold, + const double timestepGrowthThreshold, + const bool timestepGrowthRelative, + const size_t timestepGrowthWindowSize + ); +} \ No newline at end of file diff --git a/src/include/gridfire/trigger/procedures/trigger_pprint.h b/src/include/gridfire/trigger/procedures/trigger_pprint.h new file mode 100644 index 00000000..c1729462 --- /dev/null +++ b/src/include/gridfire/trigger/procedures/trigger_pprint.h @@ -0,0 +1,17 @@ +#pragma once + +#include "gridfire/trigger/trigger_result.h" + +#include + +namespace gridfire::trigger { + inline void printWhy(const TriggerResult& result, const int indent = 0) { + const std::string prefix(indent * 2, ' '); + std::cout << prefix << "• [" << (result.value ? "TRUE" : "FALSE") + << "] " << result.name << ": " << result.description << std::endl; + + for (const auto& cause : result.causes) { + printWhy(cause, indent + 1); + } + } +} \ No newline at end of file diff --git a/src/include/gridfire/trigger/trigger_abstract.h b/src/include/gridfire/trigger/trigger_abstract.h new file mode 100644 index 00000000..f2fe5ceb --- /dev/null +++ b/src/include/gridfire/trigger/trigger_abstract.h @@ -0,0 +1,26 @@ +#pragma once + +#include "gridfire/trigger/trigger_result.h" + +#include +#include + +namespace gridfire::trigger { + template + class Trigger { + public: + virtual ~Trigger() = default; + + virtual bool check(const TriggerContextStruct& ctx) const = 0; + + virtual void update(const TriggerContextStruct& ctx) = 0; + virtual void reset() = 0; + + virtual std::string name() const = 0; + virtual std::string describe() const = 0; + virtual TriggerResult why(const TriggerContextStruct& ctx) const = 0; + + virtual size_t numTriggers() const = 0; + virtual size_t numMisses() const = 0; + }; +} \ No newline at end of file diff --git a/src/include/gridfire/trigger/trigger_logical.h b/src/include/gridfire/trigger/trigger_logical.h new file mode 100644 index 00000000..b86c90ab --- /dev/null +++ b/src/include/gridfire/trigger/trigger_logical.h @@ -0,0 +1,440 @@ +#pragma once + +#include "gridfire/trigger/trigger_abstract.h" +#include "gridfire/trigger/trigger_result.h" + +#include +#include +#include + +namespace gridfire::trigger { + template + class LogicalTrigger : public Trigger {}; + + template + class AndTrigger final : public LogicalTrigger { + public: + AndTrigger(std::unique_ptr> A, std::unique_ptr> B); + ~AndTrigger() override = default; + + bool check(const TriggerContextStruct& ctx) const override; + void update(const TriggerContextStruct& ctx) override; + void reset() override; + std::string name() const override; + TriggerResult why(const TriggerContextStruct& ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + private: + std::unique_ptr> m_A; + std::unique_ptr> m_B; + + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + }; + + template + class OrTrigger final : public LogicalTrigger { + public: + OrTrigger(std::unique_ptr> A, std::unique_ptr> B); + ~OrTrigger() override = default; + + bool check(const TriggerContextStruct& ctx) const override; + void update(const TriggerContextStruct& ctx) override; + void reset() override; + std::string name() const override; + TriggerResult why(const TriggerContextStruct& ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + + private: + std::unique_ptr> m_A; + std::unique_ptr> m_B; + + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + }; + + template + class NotTrigger final : public LogicalTrigger { + public: + explicit NotTrigger(std::unique_ptr> A); + ~NotTrigger() override = default; + + bool check(const TriggerContextStruct& ctx) const override; + void update(const TriggerContextStruct& ctx) override; + void reset() override; + + std::string name() const override; + TriggerResult why(const TriggerContextStruct& ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + + private: + std::unique_ptr> m_A; + + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + + }; + + template + class EveryNthTrigger final : public LogicalTrigger { + public: + explicit EveryNthTrigger(std::unique_ptr> A, size_t N); + ~EveryNthTrigger() override = default; + + bool check(const TriggerContextStruct& ctx) const override; + void update(const TriggerContextStruct& ctx) override; + void reset() override; + + std::string name() const override; + TriggerResult why(const TriggerContextStruct& ctx) const override; + std::string describe() const override; + size_t numTriggers() const override; + size_t numMisses() const override; + + private: + std::unique_ptr> m_A; + + size_t m_N; + mutable size_t m_counter = 0; + mutable size_t m_hits = 0; + mutable size_t m_misses = 0; + mutable size_t m_updates = 0; + mutable size_t m_resets = 0; + }; + + /////////////////////////////// + // Templated Implementations // + /////////////////////////////// + + + template + AndTrigger::AndTrigger( + std::unique_ptr> A, + std::unique_ptr> B + ) : m_A(std::move(A)), m_B(std::move(B)) {} + + template + bool AndTrigger::check(const TriggerContextStruct &ctx) const { + const bool valid = m_A->check(ctx) && m_B->check(ctx); + if (valid) { + m_hits++; + } else { + m_misses++; + } + return valid; + } + + template + void AndTrigger::update(const TriggerContextStruct &ctx) { + m_A->update(ctx); + m_B->update(ctx); + m_updates++; + } + + template + void AndTrigger::reset() { + m_A->reset(); + m_B->reset(); + m_resets++; + + m_hits = 0; + m_misses = 0; + m_updates = 0; + } + + template + std::string AndTrigger::name() const { + return "AND Trigger"; + } + + template + TriggerResult AndTrigger::why(const TriggerContextStruct &ctx) const { + TriggerResult result; + + result.name = name(); + + TriggerResult A_result = m_A->why(ctx); + result.causes.push_back(A_result); + + if (!A_result.value) { + // Short Circuit + result.value = false; + result.description = "Failed because A (" + A_result.name + ") is false."; + return result; + } + + TriggerResult B_result = m_B->why(ctx); + result.causes.push_back(B_result); + + if (!B_result.value) { + result.value = false; + result.description = "Failed because B (" + B_result.name + ") is false."; + return result; + } + + result.value = true; + result.description = "Succeeded because both A (" + A_result.name + ") and B (" + B_result.description + ") are true."; + return result; + } + + template + std::string AndTrigger::describe() const { + return "(" + m_A->describe() + ") AND (" + m_B->describe() + ")"; + } + + template + size_t AndTrigger::numTriggers() const { + return m_hits; + } + + template + size_t AndTrigger::numMisses() const { + return m_misses; + } + + template + OrTrigger::OrTrigger( + std::unique_ptr> A, + std::unique_ptr> B + ) : m_A(std::move(A)), m_B(std::move(B)) {} + + template + bool OrTrigger::check(const TriggerContextStruct &ctx) const { + const bool valid = m_A->check(ctx) || m_B->check(ctx); + if (valid) { + m_hits++; + } else { + m_misses++; + } + return valid; + } + + template + void OrTrigger::update(const TriggerContextStruct &ctx) { + m_A->update(ctx); + m_B->update(ctx); + m_updates++; + } + + template + void OrTrigger::reset() { + m_A->reset(); + m_B->reset(); + m_resets++; + + m_hits = 0; + m_misses = 0; + m_updates = 0; + } + + template + std::string OrTrigger::name() const { + return "OR Trigger"; + } + + template + TriggerResult OrTrigger::why(const TriggerContextStruct &ctx) const { + TriggerResult result; + result.name = name(); + + TriggerResult A_result = m_A->why(ctx); + result.causes.push_back(A_result); + + if (A_result.value) { + // Short Circuit + result.value = true; + result.description = "Succeeded because A (" + A_result.name + ") is true."; + return result; + } + + TriggerResult B_result = m_B->why(ctx); + result.causes.push_back(B_result); + + if (B_result.value) { + result.value = true; + result.description = "Succeeded because B (" + B_result.name + ") is true."; + return result; + } + + result.value = false; + result.description = "Failed because both A (" + A_result.name + ") and B (" + B_result.name + ") are false."; + return result; + } + + template + std::string OrTrigger::describe() const { + return "(" + m_A->describe() + ") OR (" + m_B->describe() + ")"; + } + + template + size_t OrTrigger::numTriggers() const { + return m_hits; + } + + template + size_t OrTrigger::numMisses() const { + return m_misses; + } + + template + NotTrigger::NotTrigger( + std::unique_ptr> A + ) : m_A(std::move(A)) {} + + template + bool NotTrigger::check(const TriggerContextStruct &ctx) const { + const bool valid = !m_A->check(ctx); + if (valid) { + m_hits++; + } else { + m_misses++; + } + return valid; + } + + template + void NotTrigger::update(const TriggerContextStruct &ctx) { + m_A->update(ctx); + m_updates++; + } + + template + void NotTrigger::reset() { + m_A->reset(); + m_resets++; + + m_hits = 0; + m_misses = 0; + m_updates = 0; + } + + template + std::string NotTrigger::name() const { + return "NOT Trigger"; + } + + template + TriggerResult NotTrigger::why(const TriggerContextStruct &ctx) const { + TriggerResult result; + result.name = name(); + + TriggerResult A_result = m_A->why(ctx); + result.causes.push_back(A_result); + if (A_result.value) { + result.value = false; + result.description = "Failed because A (" + A_result.name + ") is true."; + return result; + } + result.value = true; + result.description = "Succeeded because A (" + A_result.name + ") is false."; + return result; + } + + template + std::string NotTrigger::describe() const { + return "NOT (" + m_A->describe() + ")"; + } + + template + size_t NotTrigger::numTriggers() const { + return m_hits; + } + + template + size_t NotTrigger::numMisses() const { + return m_misses; + } + + template + EveryNthTrigger::EveryNthTrigger(std::unique_ptr> A, const size_t N) : m_A(std::move(A)), m_N(N) { + if (N == 0) { + throw std::invalid_argument("N must be greater than 0."); + } + } + + template + bool EveryNthTrigger::check(const TriggerContextStruct &ctx) const + { + if (m_A->check(ctx) && (m_counter % m_N == 0)) { + m_hits++; + return true; + } + m_misses++; + return false; + } + + template + void EveryNthTrigger::update(const TriggerContextStruct &ctx) { + if (m_A->check(ctx)) { + m_counter++; + } + m_A->update(ctx); + m_updates++; + } + + template + void EveryNthTrigger::reset() { + m_A->reset(); + m_resets++; + + m_counter = 0; + m_hits = 0; + m_misses = 0; + m_updates = 0; + } + + template + std::string EveryNthTrigger::name() const { + return "Every Nth Trigger"; + } + + template + TriggerResult EveryNthTrigger::why(const TriggerContextStruct &ctx) const { + TriggerResult result; + result.name = name(); + TriggerResult A_result = m_A->why(ctx); + result.causes.push_back(A_result); + if (!A_result.value) { + result.value = false; + result.description = "Failed because A (" + A_result.name + ") is false."; + return result; + } + if (m_counter % m_N == 0) { + result.value = true; + result.description = "Succeeded because A (" + A_result.name + ") is true and the counter (" + std::to_string(m_counter) + ") is a multiple of N (" + std::to_string(m_N) + ")."; + return result; + } + result.value = false; + result.description = "Failed because the counter (" + std::to_string(m_counter) + ") is not a multiple of N (" + std::to_string(m_N) + ")."; + return result; + } + + template + std::string EveryNthTrigger::describe() const { + return "Every " + std::to_string(m_N) + "th (" + m_A->describe() + ")"; + } + + template + size_t EveryNthTrigger::numTriggers() const { + return m_hits; + } + + template + size_t EveryNthTrigger::numMisses() const { + return m_misses; + } + + + +} diff --git a/src/include/gridfire/trigger/trigger_result.h b/src/include/gridfire/trigger/trigger_result.h new file mode 100644 index 00000000..9b865a91 --- /dev/null +++ b/src/include/gridfire/trigger/trigger_result.h @@ -0,0 +1,13 @@ +#pragma once + +#include +#include + +namespace gridfire::trigger { + struct TriggerResult { + std::string name; + std::string description; + bool value; + std::vector causes; + }; +} \ No newline at end of file diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 634aefe1..8f43c5a1 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -176,20 +176,6 @@ namespace gridfire { } // Check the cache to see if the network needs to be repartitioned. Note that the QSECacheKey manages binning of T9, rho, and Y_full to ensure that small changes (which would likely not result in a repartitioning) do not trigger a cache miss. - const QSECacheKey key(T9, rho, Y_full); - if (! m_qse_abundance_cache.contains(key)) { - m_cacheStats.miss(CacheStats::operators::CalculateRHSAndEnergy); - LOG_ERROR( - m_logger, - "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). calculateRHSAndEnergy does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", - T9, - rho, - m_cacheStats.misses(), - m_cacheStats.hits() - ); - return std::unexpected{expectations::StaleEngineError(expectations::StaleEngineErrorTypes::SYSTEM_RESIZED)}; - } - m_cacheStats.hit(CacheStats::operators::CalculateRHSAndEnergy); const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -215,21 +201,6 @@ namespace gridfire { const double T9, const double rho ) const { - const QSECacheKey key(T9, rho, Y_full); - if (!m_qse_abundance_cache.contains(key)) { - m_cacheStats.miss(CacheStats::operators::GenerateJacobianMatrix); - LOG_ERROR( - m_logger, - "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). generateJacobianMatrix does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", - T9, - rho, - m_cacheStats.misses(), - m_cacheStats.hits() - ); - throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); - } - m_cacheStats.hit(CacheStats::operators::GenerateJacobianMatrix); - // TODO: Add sparsity pattern to this to prevent base engine from doing unnecessary work. m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); } @@ -268,27 +239,11 @@ namespace gridfire { const double T9, const double rho ) const { - const auto key = QSECacheKey(T9, rho, Y_full); - if (!m_qse_abundance_cache.contains(key)) { - m_cacheStats.miss(CacheStats::operators::CalculateMolarReactionFlow); - LOG_ERROR( - m_logger, - "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). calculateMolarReactionFlow does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", - T9, - rho, - m_cacheStats.misses(), - m_cacheStats.hits() - ); - throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); - } - m_cacheStats.hit(CacheStats::operators::CalculateMolarReactionFlow); - std::vector Y_algebraic = m_qse_abundance_cache.at(key); - - assert(Y_algebraic.size() == m_algebraic_species_indices.size()); + assert(m_Y_algebraic.size() == m_algebraic_species_indices.size()); // Fix the algebraic species to the equilibrium abundances we calculate. std::vector Y_mutable = Y_full; - for (const auto& [index, Yi] : std::views::zip(m_algebraic_species_indices, Y_algebraic)) { + for (const auto& [index, Yi] : std::views::zip(m_algebraic_species_indices, m_Y_algebraic)) { Y_mutable[index] = Yi; } @@ -309,20 +264,6 @@ namespace gridfire { const double T9, const double rho ) const { - const auto key = QSECacheKey(T9, rho, Y); - if (!m_qse_abundance_cache.contains(key)) { - m_cacheStats.miss(CacheStats::operators::GetSpeciesTimescales); - LOG_ERROR( - m_logger, - "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). getSpeciesTimescales does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", - T9, - rho, - m_cacheStats.misses(), - m_cacheStats.hits() - ); - throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); - } - m_cacheStats.hit(CacheStats::operators::GetSpeciesTimescales); const auto result = m_baseEngine.getSpeciesTimescales(Y, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -337,23 +278,9 @@ namespace gridfire { std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::getSpeciesDestructionTimescales( const std::vector &Y, - double T9, - double rho + const double T9, + const double rho ) const { - const auto key = QSECacheKey(T9, rho, Y); - if (!m_qse_abundance_cache.contains(key)) { - m_cacheStats.miss(CacheStats::operators::GetSpeciesDestructionTimescales); - LOG_ERROR( - m_logger, - "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). getSpeciesDestructionTimescales does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", - T9, - rho, - m_cacheStats.misses(), - m_cacheStats.hits() - ); - throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); - } - m_cacheStats.hit(CacheStats::operators::GetSpeciesDestructionTimescales); const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -367,16 +294,7 @@ namespace gridfire { fourdst::composition::Composition MultiscalePartitioningEngineView::update(const NetIn &netIn) { const fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn); - double T9 = netIn.temperature / 1.0e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) - const auto preKey = QSECacheKey( - T9, - netIn.density, - packCompositionToVector(baseUpdatedComposition, m_baseEngine) - ); - if (m_qse_abundance_cache.contains(preKey)) { - return baseUpdatedComposition; - } NetIn baseUpdatedNetIn = netIn; baseUpdatedNetIn.composition = baseUpdatedComposition; const fourdst::composition::Composition equilibratedComposition = equilibrateNetwork(baseUpdatedNetIn); @@ -386,15 +304,7 @@ namespace gridfire { Y_algebraic[i] = equilibratedComposition.getMolarAbundance(m_baseEngine.getNetworkSpecies()[species_index]); } - // We store the algebraic abundances in the cache for both pre- and post-conditions to avoid recalculating them. - m_qse_abundance_cache[preKey] = Y_algebraic; - - const auto postKey = QSECacheKey( - T9, - netIn.density, - packCompositionToVector(equilibratedComposition, m_baseEngine) - ); - m_qse_abundance_cache[postKey] = Y_algebraic; + m_Y_algebraic = std::move(Y_algebraic); return equilibratedComposition; } @@ -594,11 +504,6 @@ namespace gridfire { m_qse_groups.size(), m_qse_groups.size() == 1 ? "" : "s" ); - - // throw std::runtime_error( - // "Partitioning complete. Throwing an error to end the program during debugging. This error should not be caught by the caller. " - // ); - } void MultiscalePartitioningEngineView::partitionNetwork( @@ -1129,29 +1034,6 @@ namespace gridfire { coupling_flux += flow * coupling_fraction; } - // if (leakage_flux < 1e-99) { - // LOG_TRACE_L1( - // m_logger, - // "Group containing {} is in equilibrium due to vanishing leakage: leakage flux = {}, coupling flux = {}, ratio = {}", - // [&]() -> std::string { - // std::stringstream ss; - // int count = 0; - // for (const auto& idx : group.algebraic_indices) { - // ss << m_baseEngine.getNetworkSpecies()[idx].name(); - // if (count < group.species_indices.size() - 1) { - // ss << ", "; - // } - // count++; - // } - // return ss.str(); - // }(), - // leakage_flux, - // coupling_flux, - // coupling_flux / leakage_flux - // ); - // validated_groups.emplace_back(group); - // validated_groups.back().is_in_equilibrium = true; - // } else if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) { if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) { LOG_TRACE_L1( m_logger, @@ -1703,7 +1585,7 @@ namespace gridfire { } - std::string MultiscalePartitioningEngineView::QSEGroup::toString(DynamicEngine &engine) const { + std::string MultiscalePartitioningEngineView::QSEGroup::toString(const DynamicEngine &engine) const { std::stringstream ss; ss << "QSEGroup(Algebraic: ["; size_t count = 0; diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index f8afb5d9..8777b2f4 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -17,6 +17,8 @@ #include #include "fourdst/composition/exceptions/exceptions_composition.h" +#include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h" +#include "gridfire/trigger/procedures/trigger_pprint.h" namespace { @@ -78,6 +80,43 @@ namespace { namespace gridfire::solver { + CVODESolverStrategy::TimestepContext::TimestepContext( + const double t, + const N_Vector &state, + const double dt, + const double last_step_time, + const double t9, + const double rho, + const int num_steps, + const DynamicEngine &engine, + const std::vector &networkSpecies + ) : + t(t), + state(state), + dt(dt), + last_step_time(last_step_time), + T9(t9), + rho(rho), + num_steps(num_steps), + engine(engine), + networkSpecies(networkSpecies) + {} + + std::vector> CVODESolverStrategy::TimestepContext::describe() const { + std::vector> description; + description.emplace_back("t", "Current Time"); + description.emplace_back("state", "Current State Vector (N_Vector)"); + description.emplace_back("dt", "Last Timestep Size"); + description.emplace_back("last_step_time", "Time at Last Step"); + description.emplace_back("T9", "Temperature in GK"); + description.emplace_back("rho", "Density in g/cm^3"); + description.emplace_back("num_steps", "Number of Steps Taken So Far"); + description.emplace_back("engine", "Reference to the DynamicEngine"); + description.emplace_back("networkSpecies", "Reference to the list of network species"); + return description; + } + + CVODESolverStrategy::CVODESolverStrategy(DynamicEngine &engine): NetworkSolverStrategy(engine) { // TODO: In order to support MPI this function must be changed const int flag = SUNContext_Create(SUN_COMM_NULL, &m_sun_ctx); @@ -96,6 +135,8 @@ namespace gridfire::solver { } NetOut CVODESolverStrategy::evaluate(const NetIn& netIn) { + auto trigger = trigger::solver::CVODE::makeEnginePartitioningTrigger(1e12, 1e10, 1, true, 10); + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const auto absTol = m_config.get("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8); @@ -121,74 +162,63 @@ namespace gridfire::solver { double accumulated_energy = 0.0; int total_update_stages_triggered = 0; while (current_time < netIn.tMax) { - try { - user_data.T9 = T9; - user_data.rho = netIn.density; - user_data.networkSpecies = &m_engine.getNetworkSpecies(); - user_data.captured_exception.reset(); + user_data.T9 = T9; + user_data.rho = netIn.density; + user_data.networkSpecies = &m_engine.getNetworkSpecies(); + user_data.captured_exception.reset(); - check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData"); + check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData"); - int flag = -1; - if (m_stdout_logging_enabled) { - flag = CVode(m_cvode_mem, netIn.tMax, m_Y, ¤t_time, CV_ONE_STEP); - } else { - flag = CVode(m_cvode_mem, netIn.tMax, m_Y, ¤t_time, CV_NORMAL); - } + int flag = -1; + if (m_stdout_logging_enabled) { + flag = CVode(m_cvode_mem, netIn.tMax, m_Y, ¤t_time, CV_ONE_STEP); + } else { + flag = CVode(m_cvode_mem, netIn.tMax, m_Y, ¤t_time, CV_NORMAL); + } - if (user_data.captured_exception){ - std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception)); - } + if (user_data.captured_exception){ + std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception)); + } - check_cvode_flag(flag, "CVode"); + check_cvode_flag(flag, "CVode"); - long int n_steps; - double last_step_size; - CVodeGetNumSteps(m_cvode_mem, &n_steps); - CVodeGetLastStep(m_cvode_mem, &last_step_size); - long int nliters, nlcfails; - CVodeGetNumNonlinSolvIters(m_cvode_mem, &nliters); - CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nlcfails); + long int n_steps; + double last_step_size; + CVodeGetNumSteps(m_cvode_mem, &n_steps); + CVodeGetLastStep(m_cvode_mem, &last_step_size); + long int nliters, nlcfails; + CVodeGetNumNonlinSolvIters(m_cvode_mem, &nliters); + CVodeGetNumNonlinSolvConvFails(m_cvode_mem, &nlcfails); - sunrealtype* y_data = N_VGetArrayPointer(m_Y); - const double current_energy = y_data[numSpecies]; // Specific energy rate - std::cout << std::scientific << std::setprecision(3) - << "Step: " << std::setw(6) << n_steps - << " | Time: " << current_time << " [s]" - << " | Last Step Size: " << last_step_size - << " | Accumulated Energy: " << current_energy << " [erg/g]" - << " | NonlinIters: " << std::setw(2) << nliters - << " | ConvFails: " << std::setw(2) << nlcfails - << std::endl; - if (n_steps % 300 == 0) { - std::cout << "Manually triggering engine update at step " << n_steps << "..." << std::endl; - exceptions::StaleEngineTrigger::state staleState { - T9, - netIn.density, - std::vector(y_data, y_data + numSpecies), - current_time, - static_cast(n_steps), - current_energy - }; - throw exceptions::StaleEngineTrigger(staleState); - } + sunrealtype* y_data = N_VGetArrayPointer(m_Y); + const double current_energy = y_data[numSpecies]; // Specific energy rate + std::cout << std::scientific << std::setprecision(3) + << "Step: " << std::setw(6) << n_steps + << " | Time: " << current_time << " [s]" + << " | Last Step Size: " << last_step_size + << " | Accumulated Energy: " << current_energy << " [erg/g]" + << " | NonLinIters: " << std::setw(2) << nliters + << " | ConvFails: " << std::setw(2) << nlcfails + << std::endl; - // if (n_steps % 50 == 0) { - // std::cout << "Logging step diagnostics at step " << n_steps << "..." << std::endl; - // log_step_diagnostics(user_data); - // } - // if (n_steps == 300) { - // log_step_diagnostics(user_data); - // exit(0); - // } - // log_step_diagnostics(user_data); + auto ctx = TimestepContext( + current_time, + reinterpret_cast(y_data), + last_step_size, + last_callback_time, + T9, + netIn.density, + n_steps, + m_engine, + m_engine.getNetworkSpecies()); - } catch (const exceptions::StaleEngineTrigger& e) { - exceptions::StaleEngineTrigger::state staleState = e.getState(); - accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy + if (trigger->check(ctx)) { + trigger::printWhy(trigger->why(ctx)); + trigger->update(ctx); + accumulated_energy += current_energy; // Add the specific energy rate to the accumulated energy LOG_INFO( m_logger, - "Engine Update Triggered due to StaleEngineTrigger exception at time {} ({} update{} triggered). Current total specific energy {} [erg/g]", + "Engine Update Triggered at time {} ({} update{} triggered). Current total specific energy {} [erg/g]", current_time, total_update_stages_triggered, total_update_stages_triggered == 1 ? "" : "s", @@ -197,21 +227,31 @@ namespace gridfire::solver { fourdst::composition::Composition temp_comp; std::vector mass_fractions; - size_t num_species_at_stop = e.numSpecies(); + size_t num_species_at_stop = m_engine.getNetworkSpecies().size(); + + if (num_species_at_stop > m_Y->ops->nvgetlength(m_Y) - 1) { + LOG_ERROR( + m_logger, + "Number of species at engine update ({}) exceeds the number of species in the CVODE solver ({}). This should never happen.", + num_species_at_stop, + m_Y->ops->nvgetlength(m_Y) - 1 // -1 due to energy in the last index + ); + throw std::runtime_error("Number of species at engine update exceeds the number of species in the CVODE solver. This should never happen."); + } mass_fractions.reserve(num_species_at_stop); for (size_t i = 0; i < num_species_at_stop; ++i) { const auto& species = m_engine.getNetworkSpecies()[i]; temp_comp.registerSpecies(species); - mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction + mass_fractions.push_back(y_data[i] * species.mass()); // Convert from molar abundance to mass fraction } temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions); temp_comp.finalize(true); NetIn netInTemp = netIn; - netInTemp.temperature = e.temperature(); - netInTemp.density = e.density(); + netInTemp.temperature = T9 * 1e9; // Convert back to Kelvin + netInTemp.density = netIn.density; netInTemp.composition = temp_comp; fourdst::composition::Composition currentComposition = m_engine.update(netInTemp); @@ -225,13 +265,16 @@ namespace gridfire::solver { numSpecies = m_engine.getNetworkSpecies().size(); N = numSpecies + 1; + cleanup_cvode_resources(true); + + m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx); + check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); + initialize_cvode_integration_resources(N, numSpecies, current_time, currentComposition, absTol, relTol, accumulated_energy); check_cvode_flag(CVodeReInit(m_cvode_mem, current_time, m_Y), "CVodeReInit"); - } catch (fourdst::composition::exceptions::InvalidCompositionError& e) { - log_step_diagnostics(user_data); - std::rethrow_exception(std::make_exception_ptr(e)); } + } sunrealtype* y_data = N_VGetArrayPointer(m_Y); diff --git a/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp b/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp new file mode 100644 index 00000000..6d8789f4 --- /dev/null +++ b/src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp @@ -0,0 +1,263 @@ +#include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h" +#include "gridfire/solver/strategies/CVODE_solver_strategy.h" + +#include "gridfire/trigger/trigger_logical.h" +#include "gridfire/trigger/trigger_abstract.h" + +#include "quill/LogMacros.h" + +#include +#include +#include + +namespace { + template + void push_to_fixed_deque(std::deque& dq, T value, size_t max_size) { + dq.push_back(value); + if (dq.size() > max_size) { + dq.pop_front(); + } + } +} + +namespace gridfire::trigger::solver::CVODE { + SimulationTimeTrigger::SimulationTimeTrigger(double interval) : m_interval(interval) { + if (interval <= 0.0) { + LOG_ERROR(m_logger, "Interval must be positive, currently it is {}", interval); + throw std::invalid_argument("Interval must be positive, currently it is " + std::to_string(interval)); + } + } + + bool SimulationTimeTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { + if (ctx.t - m_last_trigger_time >= m_interval) { + m_hits++; + LOG_TRACE_L2(m_logger, "SimulationTimeTrigger triggered at t = {}, last trigger time was {}, delta = {}", ctx.t, m_last_trigger_time, m_last_trigger_time_delta); + return true; + } + m_misses++; + return false; + } + + void SimulationTimeTrigger::update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) { + if (check(ctx)) { + m_last_trigger_time_delta = (ctx.t - m_last_trigger_time) - m_interval; + m_last_trigger_time = ctx.t; + m_updates++; + } + } + + void SimulationTimeTrigger::reset() { + m_misses = 0; + m_hits = 0; + m_updates = 0; + m_last_trigger_time = 0.0; + m_last_trigger_time_delta = 0.0; + m_resets++; + } + + std::string SimulationTimeTrigger::name() const { + return "Simulation Time Trigger"; + } + + TriggerResult SimulationTimeTrigger::why(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { + TriggerResult result; + result.name = name(); + if (check(ctx)) { + result.value = true; + result.description = "Triggered because current time " + std::to_string(ctx.t) + " - last trigger time " + std::to_string(m_last_trigger_time - m_last_trigger_time_delta) + " >= interval " + std::to_string(m_interval); + } else { + result.value = false; + result.description = "Not triggered because current time " + std::to_string(ctx.t) + " - last trigger time " + std::to_string(m_last_trigger_time) + " < interval " + std::to_string(m_interval); + } + return result; + } + + std::string SimulationTimeTrigger::describe() const { + return "SimulationTimeTrigger(interval=" + std::to_string(m_interval) + ")"; + } + + size_t SimulationTimeTrigger::numTriggers() const { + return m_hits; + } + + size_t SimulationTimeTrigger::numMisses() const { + return m_misses; + } + + OffDiagonalTrigger::OffDiagonalTrigger( + double threshold + ) : m_threshold(threshold) { + if (threshold < 0.0) { + LOG_ERROR(m_logger, "Threshold must be non-negative, currently it is {}", threshold); + throw std::invalid_argument("Threshold must be non-negative, currently it is " + std::to_string(threshold)); + } + } + + bool OffDiagonalTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { + const size_t numSpecies = ctx.engine.getNetworkSpecies().size(); + for (int row = 0; row < numSpecies; ++row) { + for (int col = 0; col < numSpecies; ++col) { + double DRowDCol = std::abs(ctx.engine.getJacobianMatrixEntry(row, col)); + if (row != col && DRowDCol > m_threshold) { + m_hits++; + LOG_TRACE_L2(m_logger, "OffDiagonalTrigger triggered at t = {} due to entry ({}, {}) = {}", ctx.t, row, col, DRowDCol); + return true; + } + } + } + m_misses++; + return false; + } + + + void OffDiagonalTrigger::update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) { + m_updates++; + } + + void OffDiagonalTrigger::reset() { + m_misses = 0; + m_hits = 0; + m_updates = 0; + m_resets++; + } + + std::string OffDiagonalTrigger::name() const { + return "Off-Diagonal Trigger"; + } + + TriggerResult OffDiagonalTrigger::why(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { + TriggerResult result; + result.name = name(); + + if (check(ctx)) { + result.value = true; + result.description = "Triggered because an off-diagonal Jacobian entry exceeded the threshold " + std::to_string(m_threshold); + } else { + result.value = false; + result.description = "Not triggered because no off-diagonal Jacobian entry exceeded the threshold " + std::to_string(m_threshold); + } + return result; + } + + std::string OffDiagonalTrigger::describe() const { + return "OffDiagonalTrigger(threshold=" + std::to_string(m_threshold) + ")"; + } + + size_t OffDiagonalTrigger::numTriggers() const { + return m_hits; + } + + size_t OffDiagonalTrigger::numMisses() const { + return m_misses; + } + + TimestepCollapseTrigger::TimestepCollapseTrigger( + const double threshold, + const bool relative + ) : TimestepCollapseTrigger(threshold, relative, 1){} + + + TimestepCollapseTrigger::TimestepCollapseTrigger( + double threshold, + const bool relative, + const size_t windowSize + ) : m_threshold(threshold), m_relative(relative), m_windowSize(windowSize) { + if (threshold < 0.0) { + LOG_ERROR(m_logger, "Threshold must be non-negative, currently it is {}", threshold); + throw std::invalid_argument("Threshold must be non-negative, currently it is " + std::to_string(threshold)); + } + if (relative && threshold > 1.0) { + LOG_ERROR(m_logger, "Relative threshold must be between 0 and 1, currently it is {}", threshold); + throw std::invalid_argument("Relative threshold must be between 0 and 1, currently it is " + std::to_string(threshold)); + } + } + + bool TimestepCollapseTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const { + if (m_timestep_window.size() < 1) { + m_misses++; + return false; + } + double averageTimestep = 0.0; + for (const auto& dt : m_timestep_window) { + averageTimestep += dt; + } + averageTimestep /= m_timestep_window.size(); + if (m_relative && (std::abs(ctx.dt - averageTimestep) / averageTimestep) >= m_threshold) { + m_hits++; + LOG_TRACE_L2(m_logger, "TimestepCollapseTrigger triggered at t = {} due to relative growth: dt = {}, average dt = {}, threshold = {}", ctx.t, ctx.dt, averageTimestep, m_threshold); + return true; + } else if (!m_relative && std::abs(ctx.dt - averageTimestep) >= m_threshold) { + m_hits++; + LOG_TRACE_L2(m_logger, "TimestepCollapseTrigger triggered at t = {} due to absolute growth: dt = {}, average dt = {}, threshold = {}", ctx.t, ctx.dt, averageTimestep, m_threshold); + return true; + } + m_misses++; + return false; + } + + void TimestepCollapseTrigger::update(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) { + push_to_fixed_deque(m_timestep_window, ctx.dt, m_windowSize); + m_updates++; + } + + void TimestepCollapseTrigger::reset() { + m_misses = 0; + m_hits = 0; + m_updates = 0; + m_resets++; + m_timestep_window.clear(); + } + + std::string TimestepCollapseTrigger::name() const { + return "TimestepCollapseTrigger"; + } + + TriggerResult TimestepCollapseTrigger::why( + const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx + ) const { + TriggerResult result; + result.name = name(); + + if (check(ctx)) { + result.value = true; + result.description = "Triggered because timestep change exceeded the threshold " + std::to_string(m_threshold); + } else { + result.value = false; + result.description = "Not triggered because timestep change did not exceed the threshold " + std::to_string(m_threshold); + } + return result; + } + + std::string TimestepCollapseTrigger::describe() const { + return "TimestepCollapseTrigger(threshold=" + std::to_string(m_threshold) + ", relative=" + (m_relative ? "true" : "false") + ", windowSize=" + std::to_string(m_windowSize) + ")"; + } + + size_t TimestepCollapseTrigger::numTriggers() const { + return m_hits; + } + + size_t TimestepCollapseTrigger::numMisses() const { + return m_misses; + } + + std::unique_ptr> makeEnginePartitioningTrigger( + const double simulationTimeInterval, + const double offDiagonalThreshold, + const double timestepGrowthThreshold, + const bool timestepGrowthRelative, + const size_t timestepGrowthWindowSize + ) { + using ctx_t = gridfire::solver::CVODESolverStrategy::TimestepContext; + + // Create the individual conditions that can trigger a repartitioning + auto simulationTimeTrigger = std::make_unique>(std::make_unique(simulationTimeInterval), 1000); + auto offDiagTrigger = std::make_unique(offDiagonalThreshold); + auto timestepGrowthTrigger = std::make_unique>(std::make_unique(timestepGrowthThreshold, timestepGrowthRelative, timestepGrowthWindowSize), 10); + + // Combine the triggers using logical OR + auto orTriggerA = std::make_unique>(std::move(simulationTimeTrigger), std::move(offDiagTrigger)); + auto orTriggerB = std::make_unique>(std::move(orTriggerA), std::move(timestepGrowthTrigger)); + + return orTriggerB; + } +} \ No newline at end of file diff --git a/src/meson.build b/src/meson.build index adff336f..c8ca70b3 100644 --- a/src/meson.build +++ b/src/meson.build @@ -15,6 +15,7 @@ gridfire_sources = files( 'lib/io/network_file.cpp', 'lib/solver/solver.cpp', 'lib/solver/strategies/CVODE_solver_strategy.cpp', + 'lib/solver/strategies/triggers/engine_partitioning_trigger.cpp', 'lib/screening/screening_types.cpp', 'lib/screening/screening_weak.cpp', 'lib/screening/screening_bare.cpp', @@ -56,7 +57,7 @@ install_subdir('include/gridfire', install_dir: get_option('includedir')) if get_option('build-python') message('Configuring Python bindings...') - subdir('src-pybind') + subdir('python') else message('Skipping Python bindings...') endif \ No newline at end of file diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 761332b8..a63cc543 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -122,7 +122,7 @@ int main(int argc, char* argv[]){ netIn.temperature = 1.5e7; netIn.density = 1.6e2; netIn.energy = 0; - netIn.tMax = 3e13; + netIn.tMax = 3e16; // netIn.tMax = 1e-14; netIn.dt0 = 1e-12;