reaction.h

Go to the documentation of this file.

#pragma once
 
#include <string_view>
 
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
#include <unordered_map>
#include <vector>
#include <unordered_set>
 
 
#include "cppad/cppad.hpp"
#include "xxhash64.h"

namespace gridfire::reaction {
    struct RateCoefficientSet {
        double a0; 
        double a1; 
        double a2; 
        double a3; 
        double a4; 
        double a5; 
        double a6; 

        friend std::ostream& operator<<(std::ostream& os, const RateCoefficientSet& r) {
            os << "[" << r.a0 << ", " << r.a1 << ", " << r.a2 << ", "
               << r.a3 << ", " << r.a4 << ", " << r.a5 << ", " << r.a6 << "]";
            return os;
        }
    };

    class Reaction {
    public:
        virtual ~Reaction() = default;

        Reaction(
            const std::string_view id,
            const std::string_view peName,
            const int chapter,
            const std::vector<fourdst::atomic::Species> &reactants,
            const std::vector<fourdst::atomic::Species> &products,
            const double qValue,
            const std::string_view label,
            const RateCoefficientSet &sets,
            const bool reverse = false);

        [[nodiscard]] virtual double calculate_rate(const double T9) const;

        [[nodiscard]] virtual CppAD::AD<double> calculate_rate(const CppAD::AD<double> T9) const;
 
        [[nodiscard]] virtual double calculate_forward_rate_log_derivative(const double T9) const;

        [[nodiscard]] virtual std::string_view peName() const { return m_peName; }

        [[nodiscard]] int chapter() const { return m_chapter; }

        [[nodiscard]] std::string_view sourceLabel() const { return m_sourceLabel; }

        [[nodiscard]] const RateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; }

        [[nodiscard]] bool contains(const fourdst::atomic::Species& species) const;

        [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const;

        [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const;

        [[nodiscard]] std::unordered_set<fourdst::atomic::Species> all_species() const;

        [[nodiscard]] std::unordered_set<fourdst::atomic::Species> reactant_species() const;

        [[nodiscard]] std::unordered_set<fourdst::atomic::Species> product_species() const;

        [[nodiscard]] size_t num_species() const;

        [[nodiscard]] int stoichiometry(const fourdst::atomic::Species& species) const;

        [[nodiscard]] std::unordered_map<fourdst::atomic::Species, int> stoichiometry() const;

        [[nodiscard]] std::string_view id() const { return m_id; }

        [[nodiscard]] double qValue() const { return m_qValue; }

        [[nodiscard]] const std::vector<fourdst::atomic::Species>& reactants() const { return m_reactants; }

        [[nodiscard]] const std::vector<fourdst::atomic::Species>& products() const { return m_products; }

        [[nodiscard]] bool is_reverse() const { return m_reverse; }

        [[nodiscard]] double excess_energy() const;

        bool operator==(const Reaction& other) const { return m_id == other.m_id; }

        bool operator!=(const Reaction& other) const { return !(*this == other); }

        [[nodiscard]] uint64_t hash(uint64_t seed = 0) const;
 
        friend std::ostream& operator<<(std::ostream& os, const Reaction& r) {
            return os << "(Reaction:" << r.m_id << ")";
        }
 
    protected:
        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
        std::string m_id; 
        std::string m_peName; 
        int m_chapter;    
        double m_qValue = 0.0; 
        std::vector<fourdst::atomic::Species> m_reactants; 
        std::vector<fourdst::atomic::Species> m_products; 
        std::string m_sourceLabel; 
        RateCoefficientSet m_rateCoefficients; 
        bool m_reverse = false; 
    private:
        template <typename T>
        [[nodiscard]] T calculate_rate(const T T9) const {
            const T T913 = CppAD::pow(T9, 1.0/3.0);
            const T T953 = CppAD::pow(T9, 5.0/3.0);
            const T logT9 = CppAD::log(T9);
            const T exponent = m_rateCoefficients.a0 +
                   m_rateCoefficients.a1 / T9 +
                   m_rateCoefficients.a2 / T913 +
                   m_rateCoefficients.a3 * T913 +
                   m_rateCoefficients.a4 * T9 +
                   m_rateCoefficients.a5 * T953 +
                   m_rateCoefficients.a6 * logT9;
            return CppAD::exp(exponent);
        }
    };
 
 

    class LogicalReaction final : public Reaction {
    public:
        explicit LogicalReaction(const std::vector<Reaction> &reactions);

        void add_reaction(const Reaction& reaction);

        [[nodiscard]] size_t size() const { return m_rates.size(); }

        [[nodiscard]] std::vector<std::string> sources() const { return m_sources; }

        [[nodiscard]] double calculate_rate(const double T9) const override;
 
        [[nodiscard]] virtual double calculate_forward_rate_log_derivative(const double T9) const override;

        [[nodiscard]] CppAD::AD<double> calculate_rate(const CppAD::AD<double> T9) const override;

        auto begin() { return m_rates.begin(); }
        [[nodiscard]] auto begin() const { return m_rates.cbegin(); }
        auto end() { return m_rates.end(); }
        [[nodiscard]] auto end() const { return m_rates.cend(); }
 
        friend std::ostream& operator<<(std::ostream& os, const LogicalReaction& r) {
            os << "(LogicalReaction: " << r.id() << ", reverse: " << r.is_reverse() << ")";
            return os;
        }
 
    private:
        std::vector<std::string> m_sources; 
        std::vector<RateCoefficientSet> m_rates; 
 
    private:
        template <typename T>
        [[nodiscard]] T calculate_rate(const T T9) const {
            T sum = static_cast<T>(0.0);
            const T T913 = CppAD::pow(T9, 1.0/3.0);
            const T T953 = CppAD::pow(T9, 5.0/3.0);
            const T logT9 = CppAD::log(T9);
            // ReSharper disable once CppUseStructuredBinding
            for (const auto& rate : m_rates) {
                const T exponent = rate.a0 +
                       rate.a1 / T9 +
                       rate.a2 / T913 +
                       rate.a3 * T913 +
                       rate.a4 * T9 +
                       rate.a5 * T953 +
                       rate.a6 * logT9;
                sum += CppAD::exp(exponent);
            }
            return sum;
        }
    };
 
    template <typename ReactionT>
    class TemplatedReactionSet final {
    public:
        explicit TemplatedReactionSet(std::vector<ReactionT> reactions);
 
        TemplatedReactionSet();

        TemplatedReactionSet(const TemplatedReactionSet<ReactionT>& other);

        TemplatedReactionSet<ReactionT>& operator=(const TemplatedReactionSet<ReactionT>& other);

        void add_reaction(ReactionT reaction);

        void remove_reaction(const ReactionT& reaction);

        [[nodiscard]] bool contains(const std::string_view& id) const;

        [[nodiscard]] bool contains(const Reaction& reaction) const;

        [[nodiscard]] size_t size() const { return m_reactions.size(); }

        void clear();

        [[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const;

        [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const;

        [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const;

        [[nodiscard]] const ReactionT& operator[](size_t index) const;

        [[nodiscard]] const ReactionT& operator[](const std::string_view& id) const;

        bool operator==(const TemplatedReactionSet& other) const;

        bool operator!=(const TemplatedReactionSet& other) const;

        [[nodiscard]] uint64_t hash(uint64_t seed = 0) const;

        auto begin() { return m_reactions.begin(); }
        [[nodiscard]] auto begin() const { return m_reactions.cbegin(); }
        auto end() { return m_reactions.end(); }
        [[nodiscard]] auto end() const { return m_reactions.cend(); }
        friend std::ostream& operator<<(std::ostream& os, const TemplatedReactionSet<ReactionT>& r) {
            os << "(ReactionSet: [";
            size_t counter = 0;
            for (const auto& reaction : r.m_reactions) {
                os << reaction;
                if (counter < r.m_reactions.size() - 2) {
                    os << ", ";
                } else if (counter == r.m_reactions.size() - 2) {
                    os << " and ";
                }
                ++counter;
            }
            os << "])";
            return os;
        }
 
        [[nodiscard]] std::unordered_set<fourdst::atomic::Species> getReactionSetSpecies() const;
    private:
        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
        std::vector<ReactionT> m_reactions;
        std::string m_id;
        std::unordered_map<std::string, ReactionT> m_reactionNameMap; 
 
    };
 
    using ReactionSet = TemplatedReactionSet<Reaction>; 
    using LogicalReactionSet = TemplatedReactionSet<LogicalReaction>; 
 
    LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet);
 
    template <typename ReactionT>
    TemplatedReactionSet<ReactionT>::TemplatedReactionSet(
        std::vector<ReactionT> reactions
    ) :
    m_reactions(std::move(reactions)) {
        if (m_reactions.empty()) {
            return; // Case where the reactions will be added later.
        }
        m_reactionNameMap.reserve(reactions.size());
        for (const auto& reaction : m_reactions) {
            m_id += reaction.id();
            m_reactionNameMap.emplace(reaction.id(), reaction);
        }
    }
 
    template<typename ReactionT>
    TemplatedReactionSet<ReactionT>::TemplatedReactionSet() {}
 
    template <typename ReactionT>
    TemplatedReactionSet<ReactionT>::TemplatedReactionSet(const TemplatedReactionSet<ReactionT> &other) {
        m_reactions.reserve(other.m_reactions.size());
        for (const auto& reaction_ptr: other.m_reactions) {
            m_reactions.push_back(reaction_ptr);
        }
 
        m_reactionNameMap.reserve(other.m_reactionNameMap.size());
        for (const auto& reaction_ptr : m_reactions) {
            m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr);
        }
    }
 
    template <typename ReactionT>
    TemplatedReactionSet<ReactionT>& TemplatedReactionSet<ReactionT>::operator=(const TemplatedReactionSet<ReactionT> &other) {
        if (this != &other) {
            TemplatedReactionSet temp(other);
            std::swap(m_reactions, temp.m_reactions);
            std::swap(m_reactionNameMap, temp.m_reactionNameMap);
        }
        return *this;
    }
 
    template <typename ReactionT>
    void TemplatedReactionSet<ReactionT>::add_reaction(ReactionT reaction) {
        m_reactions.emplace_back(reaction);
        m_id += m_reactions.back().id();
        m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back());
    }
 
    template <typename ReactionT>
    void TemplatedReactionSet<ReactionT>::remove_reaction(const ReactionT& reaction) {
        if (!m_reactionNameMap.contains(std::string(reaction.id()))) {
            return;
        }
 
        m_reactionNameMap.erase(std::string(reaction.id()));
 
        std::erase_if(m_reactions, [&reaction](const Reaction& r) {
            return r == reaction;
        });
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::contains(const std::string_view& id) const {
        for (const auto& reaction : m_reactions) {
            if (reaction.id() == id) {
                return true;
            }
        }
        return false;
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::contains(const Reaction& reaction) const {
        for (const auto& r : m_reactions) {
            if (r == reaction) {
                return true;
            }
        }
        return false;
    }
 
    template <typename ReactionT>
    void TemplatedReactionSet<ReactionT>::clear() {
        m_reactions.clear();
        m_reactionNameMap.clear();
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::contains_species(const fourdst::atomic::Species& species) const {
        for (const auto& reaction : m_reactions) {
            if (reaction.contains(species)) {
                return true;
            }
        }
        return false;
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::contains_reactant(const fourdst::atomic::Species& species) const {
        for (const auto& r : m_reactions) {
            if (r.contains_reactant(species)) {
                return true;
            }
        }
        return false;
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::contains_product(const fourdst::atomic::Species& species) const {
        for (const auto& r : m_reactions) {
            if (r.contains_product(species)) {
                return true;
            }
        }
        return false;
    }
 
    template <typename ReactionT>
    const ReactionT& TemplatedReactionSet<ReactionT>::operator[](const size_t index) const {
        if (index >= m_reactions.size()) {
            m_logger -> flush_log();
            throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + ".");
        }
        return m_reactions[index];
    }
 
    template <typename ReactionT>
    const ReactionT& TemplatedReactionSet<ReactionT>::operator[](const std::string_view& id) const {
        if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) {
            return it->second;
        }
        m_logger -> flush_log();
        throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet.");
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::operator==(const TemplatedReactionSet<ReactionT>& other) const {
        if (size() != other.size()) {
            return false;
        }
        return hash() == other.hash();
    }
 
    template <typename ReactionT>
    bool TemplatedReactionSet<ReactionT>::operator!=(const TemplatedReactionSet<ReactionT>& other) const {
        return !(*this == other);
    }
 
    template <typename ReactionT>
    uint64_t TemplatedReactionSet<ReactionT>::hash(uint64_t seed) const {
        if (m_reactions.empty()) {
            return XXHash64::hash(nullptr, 0, seed);
        }
        std::vector<uint64_t> individualReactionHashes;
        individualReactionHashes.reserve(m_reactions.size());
        for (const auto& reaction : m_reactions) {
            individualReactionHashes.push_back(reaction.hash(seed));
        }
 
        std::ranges::sort(individualReactionHashes);
 
        const auto data = static_cast<const void*>(individualReactionHashes.data());
        const size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t);
        return XXHash64::hash(data, sizeInBytes, seed);
    }
 
    template<typename ReactionT>
    std::unordered_set<fourdst::atomic::Species> TemplatedReactionSet<ReactionT>::getReactionSetSpecies() const {
        std::unordered_set<fourdst::atomic::Species> species;
        for (const auto& reaction : m_reactions) {
            const auto reactionSpecies = reaction.all_species();
            species.insert(reactionSpecies.begin(), reactionSpecies.end());
        }
        return species;
    }
}