ThermoPropsPhase.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/AutoDiff.hpp>
#include <Reaktoro/Common/TypeOp.hpp>
#include <Reaktoro/Core/Phase.hpp>
 
namespace Reaktoro {
 
template<template<typename> typename TypeOp>
struct ThermoPropsPhaseBaseData
{
    TypeOp<real> T;
 
    TypeOp<real> P;
 
    TypeOp<ArrayXr> G0;
 
    TypeOp<ArrayXr> H0;
 
    TypeOp<ArrayXr> V0;
 
    TypeOp<ArrayXr> VT0;
 
    TypeOp<ArrayXr> VP0;
 
    TypeOp<ArrayXr> Cp0;
 
    template<template<typename> typename OtherTypeOp>
    auto operator=(const ThermoPropsPhaseBaseData<OtherTypeOp>& other)
    {
        T   = other.T;
        P   = other.P;
        G0  = other.G0;
        H0  = other.H0;
        V0  = other.V0;
        VT0 = other.VT0;
        VP0 = other.VP0;
        Cp0 = other.Cp0;
        return *this;
    }
 
    template<template<typename> typename OtherTypeOp>
    operator ThermoPropsPhaseBaseData<OtherTypeOp>()
    {
        return { T, P, G0, H0, V0, VT0, VP0, Cp0 };
    }
 
    template<template<typename> typename OtherTypeOp>
    operator ThermoPropsPhaseBaseData<OtherTypeOp>() const
    {
        return { T, P, G0, H0, V0, VT0, VP0, Cp0 };
    }
 
    auto operator=(const ArrayStream<real>& array)
    {
        array.to(T, P, G0, H0, V0, VT0, VP0, Cp0);
        return *this;
    }
 
    operator ArrayStream<real>() const
    {
        return {T, P, G0, H0, V0, VT0, VP0, Cp0};
    }
};
 
using ThermoPropsPhaseData = ThermoPropsPhaseBaseData<TypeOpIdentity>;
 
using ThermoPropsPhaseDataRef = ThermoPropsPhaseBaseData<TypeOpRef>;
 
using ThermoPropsPhaseDataConstRef = ThermoPropsPhaseBaseData<TypeOpConstRef>;
 
using ThermoPropsPhaseFn = Fn<void(ThermoPropsPhaseDataRef, const real&, const real&, ArrayXrConstRef)>;
 
template<template<typename> typename TypeOp>
class ThermoPropsPhaseBase
{
public:
    explicit ThermoPropsPhaseBase(const Phase& phase)
    : mphase(phase)
    {
        const auto numspecies = phase.species().size();
 
        mdata.G0   = ArrayXr::Zero(numspecies);
        mdata.H0   = ArrayXr::Zero(numspecies);
        mdata.V0   = ArrayXr::Zero(numspecies);
        mdata.VT0  = ArrayXr::Zero(numspecies);
        mdata.VP0  = ArrayXr::Zero(numspecies);
        mdata.Cp0  = ArrayXr::Zero(numspecies);
    }
 
    ThermoPropsPhaseBase(const Phase& phase, const ThermoPropsPhaseBaseData<TypeOp>& data)
    : mphase(phase), mdata(data)
    {}
 
    template<template<typename> typename OtherTypeOp>
    ThermoPropsPhaseBase(ThermoPropsPhaseBase<OtherTypeOp>& other)
    : mphase(other.mphase), mdata(other.mdata)
    {}
 
    template<template<typename> typename OtherTypeOp>
    ThermoPropsPhaseBase(const ThermoPropsPhaseBase<OtherTypeOp>& other)
    : mphase(other.mphase), mdata(other.mdata)
    {}
 
    auto update(const real& T, const real& P)
    {
        // Check if this update call can be skipped if T, P conditions remain the same
        if(T == mdata.T && P == mdata.P)
            return;
 
        mdata.T = T;
        mdata.P = P;
 
        auto& G0   = mdata.G0;
        auto& H0   = mdata.H0;
        auto& V0   = mdata.V0;
        auto& VT0  = mdata.VT0;
        auto& VP0  = mdata.VP0;
        auto& Cp0  = mdata.Cp0;
 
        const auto& species = mphase.species();
        const auto size = species.size();
 
        assert(G0.size()   == size);
        assert(H0.size()   == size);
        assert(V0.size()   == size);
        assert(Cp0.size()  == size);
 
        // Compute the standard thermodynamic properties of the species in the phase.
        StandardThermoProps aux;
        for(auto i = 0; i < size; ++i)
        {
            const auto& standard_thermo_model = mphase.species(i).standardThermoModel();
            aux = standard_thermo_model ? standard_thermo_model(T, P) : StandardThermoProps{};
            G0[i]  = aux.G0;
            H0[i]  = aux.H0;
            V0[i]  = aux.V0;
            VT0[i] = aux.VT0;
            VP0[i] = aux.VP0;
            Cp0[i] = aux.Cp0;
        }
    }
 
    auto update(const real& T, const real& P, Wrt<real&> wrtvar)
    {
        autodiff::seed(wrtvar);
        update(T, P);
        autodiff::unseed(wrtvar);
    }
 
    auto phase() const -> const Phase&
    {
        return mphase;
    }
 
    auto data() const -> const ThermoPropsPhaseBaseData<TypeOp>&
    {
        return mdata;
    }
 
    auto temperature() const -> real
    {
        return mdata.T;
    }
 
    auto pressure() const -> real
    {
        return mdata.P;
    }
 
    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 -> ArrayXrConstRef
    {
        return mdata.Cp0 + mdata.T * mdata.VT0 * mdata.VT0 / mdata.VP0; // from Cv0 = Cp0 + T*VT0*VT0/VP0
    }
 
    auto operator=(const ArrayStream<real>& array)
    {
        mdata = array;
        return *this;
    }
 
    operator ArrayStream<real>() const
    {
        return mdata;
    }
 
    // Ensure other ThermoPropsPhaseBase types are friend among themselves.
    template<template<typename> typename OtherTypeOp>
    friend class ThermoPropsPhaseBase;
 
private:
    Phase mphase;
 
    ThermoPropsPhaseBaseData<TypeOp> mdata;
};
 
using ThermoPropsPhase = ThermoPropsPhaseBase<TypeOpIdentity>;
 
using ThermoPropsPhaseRef = ThermoPropsPhaseBase<TypeOpRef>;
 
using ThermoPropsPhaseConstRef = ThermoPropsPhaseBase<TypeOpConstRef>;
 
} // namespace Reaktoro