EquilibriumSpecs.hpp

// Reaktoro is a unified framework for modeling chemically reactive systems.
//
// Copyright © 2014-2022 Allan Leal
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library. If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
// Reaktoro includes
#include <Reaktoro/Common/Matrix.hpp>
#include <Reaktoro/Common/Types.hpp>
#include <Reaktoro/Core/ChemicalFormula.hpp>
#include <Reaktoro/Core/ChemicalSystem.hpp>
 
namespace Reaktoro {
 
// Forward declarations
class ChemicalProps;
class ChemicalState;
 
struct ControlVariableQ
{
    using ChemicalPotentialFn = Fn<real(ChemicalProps const& props, VectorXrConstRef const& p, VectorXrConstRef const& w)>;
 
    String name;
 
    ChemicalFormula substance;
 
    String id;
 
    ChemicalPotentialFn fn;
};
 
struct ControlVariableP
{
    using ChemicalPotentialFn = Fn<real(ChemicalProps const& props, real const& pk)>;
 
    String name;
 
    ChemicalFormula substance;
 
    Index ispecies = Index(-1);
 
    ChemicalPotentialFn fn;
};
 
struct EquationConstraintFn
{
    using Func1 = Fn<real(ChemicalProps const& props, VectorXrConstRef const& p, VectorXrConstRef const& w)>;
 
    using Func2 = Fn<real(ChemicalProps const& props, VectorXrConstRef const& w)>;
 
    using Func3 = Fn<real(ChemicalProps const& props)>;
 
    EquationConstraintFn() = default;
 
    EquationConstraintFn(Func1 const& f) : _fn(f) {}
 
    EquationConstraintFn(Func2 const& f) : _fn([=](ChemicalProps const& props, VectorXrConstRef const& p, VectorXrConstRef const& w) -> real { return f(props, w); }) {}
 
    EquationConstraintFn(Func3 const& f) : _fn([=](ChemicalProps const& props, VectorXrConstRef const& p, VectorXrConstRef const& w) -> real { return f(props); }) {}
 
    auto operator()(ChemicalProps const& props, VectorXrConstRef const& p, VectorXrConstRef const& w) const -> real { return _fn(props, p, w); }
 
    auto operator=(Func1 const& f) { return *this = EquationConstraintFn(f); }
 
    auto operator=(Func2 const& f) { return *this = EquationConstraintFn(f); }
 
    auto operator=(Func3 const& f) { return *this = EquationConstraintFn(f); }
 
    auto initialized() const { return static_cast<bool>(_fn); }
 
    operator bool() const { return initialized(); }
 
    Func1 _fn;
};
 
struct EquationConstraint
{
    String id;
 
    EquationConstraintFn fn;
};
 
struct [[deprecated("ConstraintEquation has been renamed to EquationConstraint. Please make this change in your code.")]] ConstraintEquation : EquationConstraint
{};
 
struct EquationConstraints
{
    using Func = Fn<VectorXr(ChemicalProps const& props, VectorXrConstRef const& w, VectorXrConstRef const& p)>;
 
    Strings ids;
 
    Func fn;
};
 
struct ReactivityConstraint
{
    String id;
 
    VectorXd Kn;
 
    VectorXd Kp;
};
 
struct ReactivityConstraints
{
    Strings ids;
 
    MatrixXd Kn;
 
    MatrixXd Kp;
};
 
class EquilibriumSpecs
{
public:
    explicit EquilibriumSpecs(ChemicalSystem const& system);
 
    //=================================================================================================
    //
    // STATIC METHODS TO CREATE PREDEFINED CHEMICAL EQUILIBRIUM SPECIFICATIONS
    //
    //=================================================================================================
 
    static auto TP(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    static auto HP(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    static auto TV(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    static auto UV(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    static auto SP(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    static auto SV(ChemicalSystem const& system) -> EquilibriumSpecs;
 
    //=================================================================================================
    //
    // METHODS TO SPECIFY THERMODYNAMIC CONSTRAINTS
    //
    //=================================================================================================
 
    auto temperature() -> void;
 
    auto pressure() -> void;
 
    auto volume() -> void;
 
    auto internalEnergy() -> void;
 
    auto enthalpy() -> void;
 
    auto gibbsEnergy() -> void;
 
    auto helmholtzEnergy() -> void;
 
    auto entropy() -> void;
 
    auto charge() -> void;
 
    auto elementAmount(StringOrIndex const& element) -> void;
 
    auto elementAmountInPhase(StringOrIndex const& element, StringOrIndex const& phase) -> void;
 
    auto elementMass(StringOrIndex const& element) -> void;
 
    auto elementMassInPhase(StringOrIndex const& element, StringOrIndex const& phase) -> void;
 
    auto phaseAmount(StringOrIndex const& phase) -> void;
 
    auto phaseMass(StringOrIndex const& phase) -> void;
 
    auto phaseVolume(StringOrIndex const& phase) -> void;
 
    //=================================================================================================
    //
    // METHODS TO SPECIFY UNKNOWN INPUT CONDITIONS
    //
    //=================================================================================================
 
    auto unknownTemperature() -> void;
 
    auto unknownPressure() -> void;
 
    //=================================================================================================
    //
    // METHODS TO SPECIFY CHEMICAL POTENTIAL CONSTRAINTS
    //
    //=================================================================================================
 
    auto chemicalPotential(String substance) -> void;
 
    auto lnActivity(Species const& species) -> void;
 
    auto lnActivity(String species) -> void;
 
    auto lgActivity(String species) -> void;
 
    auto activity(String species) -> void;
 
    auto fugacity(String gas) -> void;
 
    auto pH() -> void;
 
    auto pMg() -> void;
 
    auto pE() -> void;
 
    auto Eh() -> void;
 
    //=================================================================================================
    //
    // METHODS TO SPECIFY HOW THE CHEMICAL SYSTEM IS OPEN
    //
    //=================================================================================================
 
    auto openTo(ChemicalFormula const& substance) -> void;
 
    // auto openTo(ChemicalState const& state) -> void;
 
    // auto openTo(Material const& material) -> void;
 
    //=================================================================================================
    //
    // METHODS TO SPECIFY ADDITIONAL UNKNOWNS
    //
    //=================================================================================================
 
    auto addUnknownTitrantAmount(ChemicalFormula const& substance) -> void;
 
    auto addUnknownChemicalPotential(String const& species) -> void;
 
    auto addUnknownStandardChemicalPotential(String const& species) -> void;
 
    auto addUnknownActivity(String const& species) -> void;
 
    auto addUnknownActivityCoefficient(String const& species) -> void;
 
    //=================================================================================================
    //
    // METHODS TO GET THE NUMBER OF INTRODUCED CONSTRAINTS, INPUT VARIABLES, AND CONTROL VARIABLES
    //
    //=================================================================================================
 
    auto numInputs() const -> Index;
 
    auto numInputParams() const -> Index;
 
    auto numControlVariables() const -> Index;
 
    auto numControlVariablesP() const -> Index;
 
    auto numControlVariablesQ() const -> Index;
 
    auto numTitrants() const -> Index;
 
    auto numTitrantsExplicit() const -> Index;
 
    auto numTitrantsImplicit() const -> Index;
 
    auto numEquationConstraints() const -> Index;
 
    auto numReactivityConstraints() const -> Index;
 
    auto numConstraints() const -> Index;
 
    auto numConservativeComponents() const -> Index;
 
    //=================================================================================================
    //
    // METHODS TO GET THE NAMES OF INTRODUCED CONSTRAINTS, INPUT VARIABLES, AND CONTROL VARIABLES
    //
    //=================================================================================================
 
    auto namesInputs() const -> Strings;
 
    auto namesInputParams() const -> Strings;
 
    auto namesControlVariables() const -> Strings;
 
    auto namesControlVariablesP() const -> Strings;
 
    auto namesControlVariablesQ() const -> Strings;
 
    auto namesTitrants() const -> Strings;
 
    auto namesTitrantsExplicit() const -> Strings;
 
    auto namesTitrantsImplicit() const -> Strings;
 
    auto namesConstraints() const -> Strings;
 
    auto namesConservativeComponents() const -> Strings;
 
    //=================================================================================================
    //
    // METHODS TO ADD CONTROL VARIABLES, CONSTRAINTS, AND INPUT VARIABLES
    //
    //=================================================================================================
 
    auto addControlVariableQ(ControlVariableQ const& qvar) -> void;
 
    auto addControlVariableP(ControlVariableP const& pvar) -> void;
 
    auto addConstraint(EquationConstraint const& constraint) -> void;
 
    auto addConstraints(EquationConstraints const& constraints) -> void;
 
    auto addReactivityConstraint(ReactivityConstraint const& constraint) -> void;
 
    auto addReactivityConstraints(ReactivityConstraints const& constraints) -> void;
 
    auto addInput(String const& var) -> Index;
 
    auto addInput(Param const& param) -> Index;
 
    //=================================================================================================
    //
    // MISCELLANEOUS METHODS
    //
    //=================================================================================================
 
    auto system() const -> ChemicalSystem const&;
 
    auto inputs() const -> Strings const&;
 
    auto params() const -> Vec<Param> const&;
 
    auto indicesInputParams() const -> Vec<Index> const&;
 
    auto isTemperatureUnknown() const -> bool;
 
    auto isPressureUnknown() const -> bool;
 
    auto indexTemperatureAmongInputVariables() const -> Index;
 
    auto indexTemperatureAmongControlVariablesP() const -> Index;
 
    auto indexPressureAmongInputVariables() const -> Index;
 
    auto indexPressureAmongControlVariablesP() const -> Index;
 
    auto indexInputVariable(String const& name) const -> Index;
 
    auto indexControlVariableP(String const& name) const -> Index;
 
    auto indexControlVariableQ(String const& name) const -> Index;
 
    auto controlVariablesQ() const -> Vec<ControlVariableQ> const&;
 
    auto controlVariablesP() const -> Vec<ControlVariableP> const&;
 
    auto titrants() const -> Vec<ChemicalFormula>;
 
    auto titrantsExplicit() const -> Vec<ChemicalFormula>;
 
    auto titrantsImplicit() const -> Vec<ChemicalFormula>;
 
    auto equationConstraintsSingle() const -> Vec<EquationConstraint> const&;
 
    auto equationConstraintsSystem() const -> Vec<EquationConstraints> const&;
 
    auto reactivityConstraintsSingle() const -> Vec<ReactivityConstraint> const&;
 
    auto reactivityConstraintsSystem() const -> Vec<ReactivityConstraints> const&;
 
    //=================================================================================================
    //
    // METHODS THAT ASSEMBLE CONSTRAINTS AND MATRICES FOR CURRENT STATE OF EQUILIBRIUM SPECIFICATIONS
    //
    //=================================================================================================
 
    auto assembleEquationConstraints() const -> EquationConstraints;
 
    auto assembleReactivityConstraints() const -> ReactivityConstraints;
 
    auto assembleReactivityConstraintsMatrixKn() const -> MatrixXd;
 
    auto assembleReactivityConstraintsMatrixKp() const -> MatrixXd;
 
    auto assembleConservationMatrix() const -> MatrixXd;
 
    auto assembleConservationMatrixN() const -> MatrixXd;
 
    auto assembleConservationMatrixQ() const -> MatrixXd;
 
    auto assembleConservationMatrixP() const -> MatrixXd;
 
private:
    ChemicalSystem m_system;
 
    Strings m_inputs;
 
    Vec<Param> m_params;
 
    Vec<Index> m_params_idxs;
 
    Vec<ControlVariableQ> qvars;
 
    Vec<ControlVariableP> pvars;
 
    Vec<EquationConstraint> econstraints_single;
 
    Vec<EquationConstraints> econstraints_system;
 
    Strings econstraints_ids;
 
    Vec<ReactivityConstraint> rconstraints_single;
 
    Vec<ReactivityConstraints> rconstraints_system;
 
    Strings rconstraints_ids;
 
    // ----- AUXILIARY DATA MEMBERS ----- //
 
    Vec<ChemicalFormula> titrants_explicit;
 
    Vec<ChemicalFormula> titrants_implicit;
 
    Strings species_with_unknown_chemical_potentials;
 
private:
    auto throwErrorIfTitrantHasBeenRegistered(ChemicalFormula const& substance) const -> void;
};
 
} // namespace Reaktoro