system.h

Go to the documentation of this file.

#ifndef CHEMMISOL_SYSTEM_H
#define CHEMMISOL_SYSTEM_H
 
#include "chemmisol/math/linear.h"
#include "reaction.h"
 
namespace chemmisol {
    namespace solver {
        class Solver;
    }
 
    struct ComponentReagent {
        double coefficient;
        ChemicalComponent* component;
 
        ComponentReagent() = default;
        ComponentReagent(
                double coefficient, ChemicalComponent* component)
            : coefficient(coefficient), component(component) {
            }
    };
 
    struct ChemicalSpeciesReagent {
        double coefficient;
        ChemicalSpecies* species;
 
        ChemicalSpeciesReagent() = default;
 
        ChemicalSpeciesReagent(
                double coefficient, ChemicalSpecies* species)
            : coefficient(coefficient), species(species) {
            }
    };
 
    class InvalidSpecies : public std::exception {
        private:
            const ChemicalSystem* chemical_system;
            const std::string& name;
            Phase phase;
            
        protected:
            InvalidSpecies(
                    const ChemicalSystem* chemical_system,
                    const std::string& name,
                    Phase phase) :
                chemical_system(chemical_system),
                name(name), phase(phase) {
                }
 
        public:
            const ChemicalSystem& getChemicalSystem() const {
                return *chemical_system;
            }
 
            const std::string& getName() const {
                return name;
            }
 
            Phase getPhase() const {
                return phase;
            }
    };
 
    class InvalidMineralSpeciesWithUndefinedSitesCount : public InvalidSpecies {
        private:
            std::string message;
 
        public:
            InvalidMineralSpeciesWithUndefinedSitesCount(
                    const ChemicalSystem* chemical_system,
                    const std::string& name
                    );
 
            const char* what() const noexcept override {
                return message.c_str();
            }
    };
 
    class ChemicalSystem {
        private:
            struct CompiledReaction {
                const Reaction* reaction;
                ChemicalSpeciesReagent produced_species;
                std::vector<ComponentReagent> components;
 
                CompiledReaction() = default;
                CompiledReaction(const Reaction* reaction)
                    : reaction(reaction) {
                    }
            };
 
            std::size_t component_index = 0;
            std::size_t species_index = 0;
            std::size_t reaction_index = 0;
 
            /* species indexes */
            std::vector<std::unique_ptr<ChemicalSpecies>> species;
            std::unordered_map<std::string, const ChemicalSpecies*> species_by_name;
 
            /* components indexes */
            std::vector<std::unique_ptr<ChemicalComponent>> components;
            std::unordered_map<std::string, const ChemicalComponent*> components_by_name;
 
            /* reaction indexes */
            std::vector<std::unique_ptr<Reaction>> reactions;
            std::unordered_map<std::string, const Reaction*> reactions_by_name;
 
            /* Internal structures initialized by the setUp() method. */
            std::vector<std::vector<double>> reaction_matrix;
            std::vector<CompiledReaction> compiled_reactions;
 
            void addSpecies(ChemicalSpecies* component, std::size_t index);
 
            void addComponent(
                    ChemicalComponent* component,
                    std::size_t species_index,
                    std::size_t component_index);
 
            const Reaction& addReaction(Reaction* reaction, std::size_t index);
 
            std::size_t max_iteration = 200;
 
            // Adsorption model parameters
            double solid_concentration = 0;
            double specific_surface_area = 0;
            double site_concentration = 0;
 
            void addComponent(
                    const std::string& name,
                    Phase phase,
                    double concentration,
                    std::size_t species_index,
                    std::size_t component_index
                    );
            void fixComponent(
                    const std::string& name,
                    Phase phase,
                    double concentration,
                    std::size_t species_index,
                    std::size_t component_index
                    );
 
            void addSpecies(
                    const std::string& name,
                    Phase phase,
                    double concentration,
                    std::size_t index
                    );
 
            const Reaction& addReaction(
                    const std::string& name,
                    double K,
                    const std::vector<Reagent>& reagents,
                    std::size_t index);
 
            void compile(const Reaction* reaction);
 
            void initReactionMatrix();
 
            template<typename T>
                void _massConservationLaw(
                        const std::vector<T>& activities,
                        std::vector<T>& result) const;
 
            template<typename T>
                T _distanceToEquilibrium(
                        const std::vector<T>& activities,
                        const Reaction& reaction) const;
 
        public:
            ChemicalSystem(const ChemicalSystem& other);
 
            /*
             * TODO: ChemicalSystem copy assignment operator.
             */
            const ChemicalSystem& operator=(const ChemicalSystem& other) = delete;
 
            /*
             * TODO: ChemicalSystem move constructor.
             */
            ChemicalSystem(ChemicalSystem&& other) = delete;
 
            /*
             * TODO: ChemicalSystem move assignment operator.
             */
            const ChemicalSystem& operator=(ChemicalSystem&& other) = delete;
 
            ChemicalSystem() = default;
 
            // TODO: possibility to initialize surface components with a non
            // null initial molar fraction.
            ChemicalSystem(
                    double solid_concentration,
                    double specific_surface_area,
                    double site_concentration
                    );
 
            void setUp();
 
            const Reaction& addReaction(
                    const std::string& name,
                    double log_K,
                    const std::vector<Reagent>& reagents
                    );
 
            void addComponent(
                    const std::string& name,
                    double total_concentration
                    );
 
            void addComponent(
                    const std::string& name,
                    Phase phase,
                    double total_concentration
                    );
            void fixComponent(
                    const std::string& name,
                    double total_concentration
                    );
 
            void fixComponent(
                    const std::string& name,
                    Phase phase,
                    double total_concentration
                    );
 
            void addSolvent(const std::string& name);
 
            void initPH(double pH, const std::string& h_component_name);
 
            void initPH(double pH) {
                initPH(pH, "H+");
            }
 
            void fixPH(double pH, const std::string& h_component_name);
 
            void fixPH(double pH) {
                fixPH(pH, "H+");
            };
 
            double getPH(const std::string& h_component_name) const;
 
            double getPH() const {
                return getPH("H+");
            }
 
            double sitesQuantity() const {
                return MineralSpecies::sites_quantity(
                        solid_concentration, specific_surface_area, site_concentration
                        );
            }
 
            void setTotalQuantity(const ChemicalComponent& component, double quantity);
 
            void setTotalConcentration(const ChemicalComponent& component, double concentration);
 
            const Reaction& getReaction(const std::string& name) const;
 
            const ChemicalSpecies& getSpecies(const std::string& name) const;
 
            const ChemicalSpecies& getSpecies(const std::size_t& index) const;
 
            const ChemicalComponent& getComponent(const std::string& name) const;
 
            const ChemicalComponent& getComponent(const std::size_t& index) const;
 
            bool isComponent(const std::string& species_name) const;
 
            bool isComponent(const ChemicalSpecies& species) const;
 
            const std::vector<std::unique_ptr<ChemicalComponent>>& getComponents() const {
                return components;
            }
 
            const std::vector<std::unique_ptr<ChemicalSpecies>>& getSpecies() const {
                return species;
            }
 
            const std::vector<std::unique_ptr<Reaction>>& getReactions() const {
                return reactions;
            }
 
            const std::vector<ComponentReagent>& getComponentReagents(
                    const Reaction& reaction) const;
 
            const ChemicalSpeciesReagent& getSpeciesReagent(
                    const Reaction& reaction) const;
 
            const std::vector<std::vector<double>>& getReactionMatrix() const {
                return reaction_matrix;
            }
 
            void proceed(const Reaction& reaction, double extent);
 
            double distanceToEquilibrium(
                    const std::vector<double>& activities,
                    const Reaction& reaction) const;
 
            std::complex<double> distanceToEquilibrium(
                    const std::vector<std::complex<double>>& activities,
                    const Reaction& reaction) const;
 
            double degree(const Reaction& reaction) const;
 
            double distanceToEquilibrium(const Reaction& reaction) const;
 
            void solveEquilibrium();
 
            void solveEquilibrium(const solver::Solver& solver);
 
            double reactionQuotient(const Reaction& reaction) const;
 
            double reactionQuotient(const std::string& name) const;
 
            void massConservationLaw(
                    const std::vector<double>& activities,
                    std::vector<double>& result) const;
 
            void massConservationLaw(
                    const std::vector<std::complex<double>>& activities,
                    std::vector<std::complex<double>>& result) const;
 
            void massConservationLaw(std::vector<double>& result) const;
 
            std::size_t getMaxIteration() const {
                return max_iteration;
            }
 
            void setMaxIteration(const std::size_t& max_iteration) {
                this->max_iteration = max_iteration;
            }
    };
 
    template<typename T>
        void ChemicalSystem::_massConservationLaw(
                const std::vector<T>& activities,
                std::vector<T>& result) const {
            for(auto& component : getComponents()) {
                if(component->isFixed()) {
                    result[component->getIndex()] = 0.0;
                } else {
                    result[component->getIndex()]
                        = component->getSpecies()->quantity(
                                activities[component->getSpecies()->getIndex()]
                                );
                }
            }
            for(auto& reaction : compiled_reactions) {
                for(const ComponentReagent& reagent
                        : reaction.components) {
                    if(!reagent.component->isFixed()) {
                        result[reagent.component->getIndex()] +=
                            reagent.coefficient / (-reaction.produced_species.coefficient)
                            * reaction.produced_species.species->quantity(
                                    activities[reaction.produced_species.species->getIndex()]
                                    );
                    }
                }
            }
            for(auto& component : getComponents()) {
                if(!component->isFixed()) {
                    result[component->getIndex()]
                        -= component->getTotalQuantity();
                }
            }
        }
 
    template<typename T>
        T ChemicalSystem::_distanceToEquilibrium(
                const std::vector<T>& activities,
                const Reaction& reaction) const {
            T reactives = {reaction.getK()};
            T products {1.0};
            auto& compiled_reaction = compiled_reactions[reaction.getIndex()];
            for(const auto& reagent : compiled_reaction.components) {
                T activity = activities[reagent.component->getSpecies()->getIndex()];
                if(activity != 0.0) {
                    if(reagent.coefficient > 0) {
                        reactives *= std::pow(activity, reagent.coefficient);
                    } else {
                        products *= std::pow(activity, -reagent.coefficient);
                    }
                } else {
                    if(reagent.coefficient > 0) {
                        reactives = 0.0;
                    } else {
                        products = 0.0;
                    }
                }
            }
            T activity = activities[compiled_reaction.produced_species.species->getIndex()];
            if(activity != 0.0)
                // This coefficient is necessarily negative
                products *= std::pow(
                        activity,
                        -compiled_reaction.produced_species.coefficient);
            else
                products = 0.0;
            return products - reactives;
        }
}
#endif /*CHEMMISOL_SYSTEM_H*/