ChemicalPropsPhase.hpp

// Reaktoro is a unified framework for modeling chemically reactive phases.
//
// Copyright © 2014-2024 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/ArrayStream.hpp>
#include <Reaktoro/Common/Constants.hpp>
#include <Reaktoro/Common/TypeOp.hpp>
#include <Reaktoro/Core/Phase.hpp>
#include <Reaktoro/Core/StateOfMatter.hpp>
 
namespace Reaktoro {
 
template<template<typename> typename TypeOp>
struct ChemicalPropsPhaseBaseData
{
    TypeOp<real> T;
 
    TypeOp<real> P;
 
    TypeOp<ArrayXr> n;
 
    TypeOp<real> nsum;
 
    TypeOp<real> msum;
 
    TypeOp<ArrayXr> x;
 
    TypeOp<ArrayXr> G0;
 
    TypeOp<ArrayXr> H0;
 
    TypeOp<ArrayXr> V0;
 
    TypeOp<ArrayXr> VT0;
 
    TypeOp<ArrayXr> VP0;
 
    TypeOp<ArrayXr> Cp0;
 
    TypeOp<real> Vx;
 
    TypeOp<real> VxT;
 
    TypeOp<real> VxP;
 
    TypeOp<ArrayXr> Vxi;
 
    TypeOp<real> Gx;
 
    TypeOp<real> Hx;
 
    TypeOp<real> Cpx;
 
    TypeOp<ArrayXr> ln_g;
 
    TypeOp<ArrayXr> ln_a;
 
    TypeOp<ArrayXr> u;
 
    TypeOp<StateOfMatter> som;
 
    template<template<typename> typename OtherTypeOp>
    auto operator=(const ChemicalPropsPhaseBaseData<OtherTypeOp>& other)
    {
        T    = other.T;
        P    = other.P;
        n    = other.n;
        nsum = other.nsum;
        msum = other.msum;
        x    = other.x;
        G0   = other.G0;
        H0   = other.H0;
        V0   = other.V0;
        VT0  = other.VT0;
        VP0  = other.VP0;
        Cp0  = other.Cp0;
        Vx   = other.Vx;
        VxT  = other.VxT;
        VxP  = other.VxP;
        Vxi  = other.Vxi;
        Gx   = other.Gx;
        Hx   = other.Hx;
        Cpx  = other.Cpx;
        ln_g = other.ln_g;
        ln_a = other.ln_a;
        u    = other.u;
        som  = other.som;
        return *this;
    }
 
    template<template<typename> typename OtherTypeOp>
    operator ChemicalPropsPhaseBaseData<OtherTypeOp>()
    {
        return { T, P, n, nsum, msum, x, G0, H0, V0, VT0, VP0, Cp0, Vx, VxT, VxP, Vxi, Gx, Hx, Cpx, ln_g, ln_a, u, som };
    }
 
    template<template<typename> typename OtherTypeOp>
    operator ChemicalPropsPhaseBaseData<OtherTypeOp>() const
    {
        return { T, P, n, nsum, msum, x, G0, H0, V0, VT0, VP0, Cp0, Vx, VxT, VxP, Vxi, Gx, Hx, Cpx, ln_g, ln_a, u, som };
    }
 
    auto operator=(const ArrayStream<real>& array)
    {
        array.to(T, P, n, nsum, msum, x, G0, H0, V0, VT0, VP0, Cp0, Vx, VxT, VxP, Vxi, Gx, Hx, Cpx, ln_g, ln_a, u, som);
        return *this;
    }
 
    operator ArrayStream<real>() const
    {
        return { T, P, n, nsum, msum, x, G0, H0, V0, VT0, VP0, Cp0, Vx, VxT, VxP, Vxi, Gx, Hx, Cpx, ln_g, ln_a, u, som };
    }
};
 
using ChemicalPropsPhaseData = ChemicalPropsPhaseBaseData<TypeOpIdentity>;
 
using ChemicalPropsPhaseDataRef = ChemicalPropsPhaseBaseData<TypeOpRef>;
 
using ChemicalPropsPhaseDataConstRef = ChemicalPropsPhaseBaseData<TypeOpConstRef>;
 
using ChemicalPropsPhaseFn = Fn<void(ChemicalPropsPhaseDataRef, const real&, const real&, ArrayXrConstRef)>;
 
template<template<typename> typename TypeOp>
class ChemicalPropsPhaseBase
{
public:
    explicit ChemicalPropsPhaseBase(const Phase& phase)
    : mphase(phase)
    {
        const auto N = phase.species().size();
 
        mdata.n    = ArrayXr::Zero(N);
        mdata.x    = ArrayXr::Zero(N);
        mdata.G0   = ArrayXr::Zero(N);
        mdata.H0   = ArrayXr::Zero(N);
        mdata.V0   = ArrayXr::Zero(N);
        mdata.VT0  = ArrayXr::Zero(N);
        mdata.VP0  = ArrayXr::Zero(N);
        mdata.Cp0  = ArrayXr::Zero(N);
        mdata.Vxi  = ArrayXr::Zero(N);
        mdata.ln_g = ArrayXr::Zero(N);
        mdata.ln_a = ArrayXr::Zero(N);
        mdata.u    = ArrayXr::Zero(N);
    }
 
    ChemicalPropsPhaseBase(const Phase& phase, const ChemicalPropsPhaseBaseData<TypeOp>& data)
    : mphase(phase), mdata(data)
    {}
 
    template<template<typename> typename OtherTypeOp>
    ChemicalPropsPhaseBase(ChemicalPropsPhaseBase<OtherTypeOp>& other)
    : mphase(other.mphase), mdata(other.mdata)
    {}
 
    template<template<typename> typename OtherTypeOp>
    ChemicalPropsPhaseBase(const ChemicalPropsPhaseBase<OtherTypeOp>& other)
    : mphase(other.mphase), mdata(other.mdata)
    {}
 
    auto update(const real& T, const real& P, ArrayXrConstRef n, Map<String, Any>& extra)
    {
        _update<false>(T, P, n, extra);
    }
 
    auto updateIdeal(const real& T, const real& P, ArrayXrConstRef n, Map<String, Any>& extra)
    {
        _update<true>(T, P, n, extra);
    }
 
    auto updateWithData(const ChemicalPropsPhaseBaseData<TypeOp>& data)
    {
        mdata = data;
    }
 
    auto phase() const -> const Phase&
    {
        return mphase;
    }
 
    auto data() const -> const ChemicalPropsPhaseBaseData<TypeOp>&
    {
        return mdata;
    }
 
    auto stateOfMatter() const -> StateOfMatter
    {
        return mdata.som;
    }
 
    auto temperature() const -> real
    {
        return mdata.T;
    }
 
    auto pressure() const -> real
    {
        return mdata.P;
    }
 
    auto speciesAmounts() const -> ArrayXrConstRef
    {
        return mdata.n;
    }
 
    auto speciesMoleFractions() const -> ArrayXrConstRef
    {
        return mdata.x;
    }
 
    auto speciesActivityCoefficientsLn() const -> ArrayXrConstRef
    {
        return mdata.ln_g;
    }
 
    auto speciesActivitiesLn() const -> ArrayXrConstRef
    {
        return mdata.ln_a;
    }
 
    auto speciesChemicalPotentials() const -> ArrayXrConstRef
    {
        return mdata.u;
    }
 
    auto speciesPartialMolarVolumes() const -> ArrayXrConstRef
    {
        return mdata.V0 + mdata.Vxi;
    }
 
    auto speciesStandardVolumes() const -> ArrayXrConstRef
    {
        return mdata.V0;
    }
 
    auto speciesStandardVolumesT() const -> ArrayXrConstRef
    {
        return mdata.VT0;
    }
 
    auto speciesStandardVolumesP() const -> ArrayXrConstRef
    {
        return mdata.VP0;
    }
 
    auto speciesStandardGibbsEnergies() const -> ArrayXrConstRef
    {
        return mdata.G0;
    }
 
    auto speciesStandardEnthalpies() const -> ArrayXrConstRef
    {
        return mdata.H0;
    }
 
    auto speciesStandardEntropies() const -> ArrayXr
    {
        return (mdata.H0 - mdata.G0)/mdata.T; // from G0 = H0 - T*S0
    }
 
    auto speciesStandardInternalEnergies() const -> ArrayXr
    {
        return mdata.H0 - mdata.P * mdata.V0; // from H0 = U0 + P*V0
    }
 
    auto speciesStandardHelmholtzEnergies() const -> ArrayXr
    {
        return mdata.G0 - mdata.P * mdata.V0; // from A0 = U0 - T*S0 = (H0 - P*V0) + (G0 - H0) = G0 - P*V0
    }
 
    auto speciesStandardHeatCapacitiesConstP() const -> ArrayXrConstRef
    {
        return mdata.Cp0;
    }
 
    auto speciesStandardHeatCapacitiesConstV() const -> ArrayXr
    {
        const auto& Cp0 = mdata.Cp0;
        const auto& T   = mdata.T;
        const auto& VT0 = mdata.VT0;
        const auto& VP0 = mdata.VP0;
        return (VP0 == 0.0).select(Cp0, Cp0 + T*VT0*VT0/VP0); // from Cv0 = Cp0 + T*VT0*VT0/VP0
    }
 
    auto molarMass() const -> real
    {
        return mdata.msum / mdata.nsum;
    }
 
    auto molarVolume() const -> real
    {
        return (mdata.x * mdata.V0).sum() + mdata.Vx;
    }
 
    auto molarVolumeT() const -> real
    {
        return (mdata.x * mdata.VT0).sum() + mdata.VxT;
    }
 
    auto molarVolumeP() const -> real
    {
        return (mdata.x * mdata.VP0).sum() + mdata.VxP;
    }
 
    auto molarGibbsEnergy() const -> real
    {
        return (mdata.x * mdata.G0).sum() + mdata.Gx;
    }
 
    auto molarEnthalpy() const -> real
    {
        return (mdata.x * mdata.H0).sum() + mdata.Hx;
    }
 
    auto molarEntropy() const -> real
    {
        const auto S0 = (mdata.H0 - mdata.G0)/mdata.T; // from G0 = H0 - T*S0
        const auto Sx = (mdata.Hx - mdata.Gx)/mdata.T; // from Gx = Hx - T*Sx
        return (mdata.x * S0).sum() + Sx;
    }
 
    auto molarInternalEnergy() const -> real
    {
        const auto U0 = mdata.H0 - mdata.P * mdata.V0; // from H0 = U0 + P*V0
        const auto Ux = mdata.Hx - mdata.P * mdata.Vx; // from Hx = U0 + P*Vx
        return (mdata.x * U0).sum() + Ux;
    }
 
    auto molarHelmholtzEnergy() const -> real
    {
        const auto A0 = mdata.G0 - mdata.P * mdata.V0; // from A0 = U0 - T*S0 = (H0 - P*V0) + (G0 - H0) = G0 - P*V0
        const auto Ax = mdata.Gx - mdata.P * mdata.Vx; // from Ax = U0 - T*Sx = (Hx - P*Vx) + (Gx - Hx) = Gx - P*Vx
        return (mdata.x * A0).sum() + Ax;
    }
 
    auto molarHeatCapacityConstP() const -> real
    {
        return (mdata.x * mdata.Cp0).sum() + mdata.Cpx;
    }
 
    auto molarHeatCapacityConstV() const -> real
    {
        const auto Cp = molarHeatCapacityConstP();
        const auto VT = molarVolumeT();
        const auto VP = molarVolumeP();
        const auto T  = temperature();
        return VP == 0.0 ? Cp : Cp + T*VT*VT/VP;
    }
 
    auto specificVolume() const -> real
    {
        return molarVolume() / molarMass();
    }
 
    auto specificVolumeT() const -> real
    {
        return molarVolumeT() / molarMass();
    }
 
    auto specificVolumeP() const -> real
    {
        return molarVolumeP() / molarMass();
    }
 
    auto specificGibbsEnergy() const -> real
    {
        return molarGibbsEnergy() / molarMass();
    }
 
    auto specificEnthalpy() const -> real
    {
        return molarEnthalpy() / molarMass();
    }
 
    auto specificEntropy() const -> real
    {
        return molarEntropy() / molarMass();
    }
 
    auto specificInternalEnergy() const -> real
    {
        return molarInternalEnergy() / molarMass();
    }
 
    auto specificHelmholtzEnergy() const -> real
    {
        return molarHelmholtzEnergy() / molarMass();
    }
 
    auto specificHeatCapacityConstP() const -> real
    {
        return molarHeatCapacityConstP() / molarMass();
    }
 
    auto specificHeatCapacityConstV() const -> real
    {
        return molarHeatCapacityConstV() / molarMass();
    }
 
    auto density() const -> real
    {
        return molarMass() / molarVolume();
    }
 
    auto amount() const -> real
    {
        return mdata.nsum;
    }
 
    auto mass() const -> real
    {
        return mdata.msum;
    }
 
    auto volume() const -> real
    {
        return molarVolume() * amount();
    }
 
    auto volumeT() const -> real
    {
        return molarVolumeT() * amount();
    }
 
    auto volumeP() const -> real
    {
        return molarVolumeP() * amount();
    }
 
    auto gibbsEnergy() const -> real
    {
        return molarGibbsEnergy() * amount();
    }
 
    auto enthalpy() const -> real
    {
        return molarEnthalpy() * amount();
    }
 
    auto entropy() const -> real
    {
        return molarEntropy() * amount();
    }
 
    auto internalEnergy() const -> real
    {
        return molarInternalEnergy() * amount();
    }
 
    auto helmholtzEnergy() const -> real
    {
        return molarHelmholtzEnergy() * amount();
    }
 
    auto heatCapacityConstP() const -> real
    {
        return molarHeatCapacityConstP() * amount();
    }
 
    auto heatCapacityConstV() const -> real
    {
        return molarHeatCapacityConstV() * amount();
    }
 
    auto soundSpeed() const -> real
    {
        const auto rho = density();
        const auto Cp = molarHeatCapacityConstP();
        const auto Cv = molarHeatCapacityConstV();
        const auto VP = molarVolumeP();
        const auto MM = molarMass();
        return sqrt(-MM * Cp/Cv / VP) / rho;
    }
 
    auto operator=(const ArrayStream<real>& array)
    {
        mdata = array;
        return *this;
    }
 
    operator ArrayStream<real>() const
    {
        return mdata;
    }
 
    // Ensure other ChemicalPropsPhaseBase types are friend among themselves.
    template<template<typename> typename OtherTypeOp>
    friend class ChemicalPropsPhaseBase;
 
private:
    Phase mphase;
 
    ChemicalPropsPhaseBaseData<TypeOp> mdata;
 
private:
 
    template<bool use_ideal_activity_model>
    auto _update(const real& T, const real& P, ArrayXrConstRef n, Map<String, Any>& extra)
    {
        mdata.T = T;
        mdata.P = P;
        mdata.n = n;
 
        const auto R = universalGasConstant;
 
        auto& nsum = mdata.nsum;
        auto& msum = mdata.msum;
        auto& x    = mdata.x;
        auto& G0   = mdata.G0;
        auto& H0   = mdata.H0;
        auto& V0   = mdata.V0;
        auto& VT0  = mdata.VT0;
        auto& VP0  = mdata.VP0;
        auto& Cp0  = mdata.Cp0;
        auto& Vx   = mdata.Vx;
        auto& VxT  = mdata.VxT;
        auto& VxP  = mdata.VxP;
        auto& Gx   = mdata.Gx;
        auto& Hx   = mdata.Hx;
        auto& Cpx  = mdata.Cpx;
        auto& Vxi = mdata.Vxi;
        auto& ln_g = mdata.ln_g;
        auto& ln_a = mdata.ln_a;
        auto& u    = mdata.u;
        auto& som  = mdata.som;
 
        const auto& species = phase().species();
        const auto N = species.size();
 
        assert(    n.size() == N );
        assert(   G0.size() == N );
        assert(   H0.size() == N );
        assert(   V0.size() == N );
        assert(  VT0.size() == N );
        assert(  VP0.size() == N );
        assert(  Cp0.size() == N );
        assert( ln_g.size() == N );
        assert( ln_a.size() == N );
        assert(    u.size() == N );
        assert(   Vxi.size() == N );
 
        // Compute the standard thermodynamic properties of the species in the phase.
        StandardThermoProps aux;
        for(auto i = 0; i < N; ++i)
        {
            aux = species[i].standardThermoProps(T, P);
            G0[i]  = aux.G0;
            H0[i]  = aux.H0;
            V0[i]  = aux.V0;
            VT0[i] = aux.VT0;
            VP0[i] = aux.VP0;
            Cp0[i] = aux.Cp0;
        }
 
        // Compute the amount of the phase
        nsum = n.sum();
 
        // Compute the mass of the phase
        msum = (n * mphase.speciesMolarMasses()).sum();
 
        // Compute the mole fractions of the species
        if(nsum == 0.0)
            x = (N == 1) ? 1.0 : 0.0;
        else x = n / nsum;
 
        // Ensure there are no zero mole fractions
        error(x.minCoeff() == 0.0, "Could not compute the chemical properties of phase ",
            phase().name(), " because it has one or more species with zero amounts.");
 
        // Compute the activity properties of the phase
        ActivityPropsRef aprops{ Vx, VxT, VxP, Vxi, Gx, Hx, Cpx, ln_g, ln_a, som, extra };
        ActivityModelArgs args{ T, P, x };
        const ActivityModel& activity_model = use_ideal_activity_model ?  // IMPORTANT: Use `const ActivityModel&` here instead of `ActivityModel`, otherwise a new model is constructed without cache, and so memoization will not take effect.
            phase().idealActivityModel() : phase().activityModel();
 
        if(nsum == 0.0) aprops = 0.0;
        else activity_model(aprops, args);
 
        // Compute the chemical potentials of the species
        u = G0 + R*T*ln_a;
    }
};
 
using ChemicalPropsPhase = ChemicalPropsPhaseBase<TypeOpIdentity>;
 
using ChemicalPropsPhaseRef = ChemicalPropsPhaseBase<TypeOpRef>;
 
using ChemicalPropsPhaseConstRef = ChemicalPropsPhaseBase<TypeOpConstRef>;
 
} // namespace Reaktoro