Skip to content

Commit

Permalink
Origin bucket (#18)
Browse files Browse the repository at this point in the history 
* Database object can now be constructed using a file or a string containing a ThermoDataSet in valid JSON

* fixed division by zero problem in the propagation of derivatives

* fix calculation of derivatives in Cp f(T) method at reference T-P

* example properties reaction

* bug-fix entropy of reaction
@gdmiron
gdmiron authored Feb 17, 2020

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
GPG key ID: 4AEE18F83AFDEB23
Expired
Verified
Learn about vigilant mode
1 parent ae8feb6 commit e3c2458
Showing 7 changed files with 157 additions and 129 deletions.
8 changes: 4 additions & 4 deletions ThermoFun/Common/ThermoScalar.hpp
View file
Original file line number Diff line number Diff line change
@@ -403,7 +403,7 @@ inline auto sqrt(const ThermoScalarBase<V>& l) -> ThermoScalarBase<double>
const double tmp1 = std::sqrt(l.val);
const double tmp2 = 0.5 * tmp1/l.val;
if (l.val == 0)
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, 0, status(l)};
return {tmp1, /*tmp2 * l.ddt, tmp2 * l.ddp,*/ 0.0,0.0,0.0, status(l)};
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, 0.5*(l.err/l.val), status(l)};
}

@@ -414,7 +414,7 @@ inline auto pow(const ThermoScalarBase<V>& l, double power) -> ThermoScalarBase<
const double tmp1 = std::pow(l.val, power);
const double tmp2 = power * tmp1/l.val;
if (l.val == 0)
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, 0.0, status(l)};
return {tmp1, /*tmp2 * l.ddt, tmp2 * l.ddp,*/0.0,0.0, 0.0, status(l)};
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, std::fabs(power)*(l.err/l.val), status(l)};
}

@@ -426,7 +426,7 @@ inline auto pow(const ThermoScalarBase<VL>& l, const ThermoScalarBase<VR>& power
const double powl = std::pow(l.val, power.val);
const double tmp = power.val/l.val;
if (l.val == 0)
return {powl, powl * (logl * power.ddt + tmp * l.ddt), powl * (logl * power.ddp + tmp * l.ddp), 0.0, status(l,power)};
return {powl, powl * (logl * power.ddt + /*tmp * l.ddt*/0.0), powl * (logl * power.ddp + /*tmp * l.ddp*/0.0), 0.0, status(l,power)};
return {powl, powl * (logl * power.ddt + tmp * l.ddt), powl * (logl * power.ddp + tmp * l.ddp), powl*(l.err/l.val), status(l,power)};
}

@@ -445,7 +445,7 @@ inline auto log(const ThermoScalarBase<V>& l) -> ThermoScalarBase<double>
const double tmp1 = std::log(l.val);
const double tmp2 = 1.0/l.val;
if (l.val == 0)
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, 0.0, status(l)};
return {tmp1, /*tmp2 * l.ddt, tmp2 * l.ddp,*/ 0.0,0.0, 0.0, status(l)};
return {tmp1, tmp2 * l.ddt, tmp2 * l.ddp, 0.434*(l.err/l.val), status(l)};
}

20 changes: 15 additions & 5 deletions ThermoFun/Database.cpp
View file
Original file line number Diff line number Diff line change
@@ -259,16 +259,26 @@ struct Database::Impl
}
}

/// Parses the JSON file and puts the data into the internal data structure
/// Parses the JSON file (or string) and puts the data into the internal data structure
/// @param filename name of the file (in the working directory)
auto fromFile(std::string filename) -> void
{
try {
std::ifstream ifs(filename);
if (!ifs.good())
funError("File reading error", std::string("Database file "+ filename +" not found!"), __LINE__, __FILE__);

json j = json::parse(ifs);
json j;

try
{
// trying to see if the string is a valid json
j = json::parse(filename);
}
catch (json::parse_error &e)
{
std::ifstream ifs(filename);
if (!ifs.good())
funError("File reading error", std::string("Database file " + filename + " not found!"), __LINE__, __FILE__);
j = json::parse(ifs);
}

if (j.contains("elements"))
addRecords(j["elements"], "element");
11 changes: 7 additions & 4 deletions ThermoFun/Database.h
View file
Original file line number Diff line number Diff line change
@@ -29,15 +29,18 @@ class Database
/// Construct default database instance
Database();

/// Construct a database instance by parsing a "json" file
/// containg the exported substances and reactions
/**
* @brief Construct a new Database object
*
* @param filename name/path of the file or a string containing a ThermoDataSet in JSON format
*/
explicit Database(std::string filename);

/**
* @brief Database constructs a database instace from a vector of records in json format
* Records with the same symbol will be overwritten!
* @param jsonRecords vector of records in JSON string format
* @param _label, oprional, (element, substance, reactions),
* @param _label, optional, (element, substance, reactions),
* used when the vector of records are of one type and do not contain themselves the key "_label"
*/
Database(std::vector<std::string> jsonRecords, std::string _label);
@@ -51,7 +54,7 @@ class Database
/**
* @brief appendData append records to the database from a file.
* Records with the same symbol will be overwritten!
* @param filename path to the file with recors
* @param filename name/path of the file or a string containing a ThermoDataSet in JSON format
*/
auto appendData(std::string filename) -> void;

221 changes: 114 additions & 107 deletions ThermoFun/Substances/EmpiricalCpIntegration.cpp
View file
Original file line number Diff line number Diff line change
@@ -5,25 +5,26 @@
#include "Substance.h"
#include <iomanip>

namespace ThermoFun {
namespace ThermoFun
{

auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pressure /*Pbar*/, Substance substance) -> ThermoPropertiesSubstance
{
ThermoPropertiesSubstance thermo_properties_PT = substance.thermoProperties();
ThermoPropertiesSubstance thermo_properties_PrTr = substance.thermoReferenceProperties();
SubstanceClass::type substance_class = substance.substanceClass();
ThermoParametersSubstance thermo_parameters = substance.thermoParameters();
ThermoPropertiesSubstance thermo_properties_PT = substance.thermoProperties();
ThermoPropertiesSubstance thermo_properties_PrTr = substance.thermoReferenceProperties();
SubstanceClass::type substance_class = substance.substanceClass();
ThermoParametersSubstance thermo_parameters = substance.thermoParameters();

Reaktoro_::ThermoScalar V;
int k=-1;
int k = -1;
vector<double> ac;
auto TK_=TK;
for (unsigned i = 0; i <16; i++)
auto TK_ = TK;
for (unsigned i = 0; i < 16; i++)
{
ac.push_back(0.0);
}

auto TrK = substance.referenceT()/* + C_to_K*/;
auto TrK = substance.referenceT() /* + C_to_K*/;

auto S = thermo_properties_PrTr.entropy;
auto G = thermo_properties_PrTr.gibbs_energy;
@@ -35,16 +36,8 @@ auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pres
return thermo_properties_PrTr;
}

// aW.twp->devG = dc[q].Gs[1];
// ??????????????????????????????????????'
// if( !dc[q].Cp )
// {
// aW.twp->Cp = (double)dc[q].Cps[0];
// goto NEXT;
// }

// get Cp interval -> this has to go!!!!
for (size_t i=0; i<thermo_parameters.temperature_intervals.size(); i++)
for (size_t i = 0; i < thermo_parameters.temperature_intervals.size(); i++)
{
if (thermo_parameters.temperature_intervals[i].size() > 0)
{
@@ -60,120 +53,134 @@ auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pres
}
}

if (k<0)
if (k < 0)
{
if (TK_ < thermo_parameters.temperature_intervals[0][0])
k = 0;
if (TK_ > thermo_parameters.temperature_intervals[thermo_parameters.temperature_intervals.size()-1][1])
k = thermo_parameters.temperature_intervals.size()-1;
std::cout << "The given temperature: "<< TK_ <<" is not inside the specified interval/s for the Cp calculation.";
std::cout << "The temperature is not inside the specified interval for the substance "<< substance.symbol() << "." << std::endl << __FILE__ << std::endl << __LINE__;
// Exception exception;
// exception.error << "The given temperature: "<< TK_ <<" is not inside the specified interval/s for the Cp calculation.";
// exception.reason << "The temperature is not inside the specified interval for the substance "<< substance.symbol() << ".";
// exception.file = __FILE__;
// exception.line = __LINE__;
// RaiseError(exception)
if (TK_ > thermo_parameters.temperature_intervals[thermo_parameters.temperature_intervals.size() - 1][1])
k = thermo_parameters.temperature_intervals.size() - 1;
std::cout << "The given temperature: " << TK_ << " is not inside the specified interval/s for the Cp calculation.";
std::cout << "The temperature is not inside the specified interval for the substance " << substance.symbol() << "." << std::endl
<< __FILE__ << std::endl
<< __LINE__;
}

k = 0;


k=0;

for (unsigned i=0; i<thermo_parameters.Cp_coeff[k].size(); i++)
for (unsigned i = 0; i < thermo_parameters.Cp_coeff[k].size(); i++)
{
if (i== 16) break;
if (i == 16)
break;
ac[i] = thermo_parameters.Cp_coeff[k][i];
}

// ac = thermo_parameters.Cp_coeff[k];

auto Cp = ( ac[0] + ac[1]*TK + ac[2]/(TK*TK) + ac[3]/(pow(TK,0.5)) + ac[4]*(TK*TK)
+ ac[5]*(TK*TK*TK) + ac[6]*(TK*TK*TK*TK) + ac[7]/(TK*TK*TK) + ac[8]/TK + ac[9]*(pow(TK,(1/2))) /*+ ac[10]*log(T)*/);

// Phase transitions
if (fabs(TK.val-TrK) > TEMPER_PREC)
auto Cp = (ac[0] +
ac[1] * TK +
ac[2] / (TK * TK) +
ac[3] / (pow(TK, 0.5)) +
ac[4] * (TK * TK) +
ac[5] * (TK * TK * TK) +
ac[6] * (TK * TK * TK * TK) +
ac[7] / (TK * TK * TK) +
ac[8] / TK +
ac[9] * (pow(TK, (1 / 2))) /*+ ac[10]*log(T)*/);

for (unsigned j = 0, ft = 0; j <= k; j++)
{
for (unsigned j=0, ft = 0; j<=k; j++)
{
if ( j == k )
TK = TK_.val/* + C_to_K*/; // current T is the end T for phase transition Cp calculations
else TK = thermo_parameters.temperature_intervals[j][1] /*+ C_to_K*/; // takes the upper bound from the j-th Tinterval

if( !j )
TrK = substance.referenceT() /*+ C_to_K*/; // if j=0 the first interval should contain the reference T (Tcr)
else TrK = thermo_parameters.temperature_intervals[j][0] /*+ C_to_K*/; // if j>0 then we are in a different Tinterval and the reference T becomes the lower bound of the interval

auto Tst2 = TrK * TrK;
auto Tst3 = Tst2 * TrK;
auto Tst4 = Tst3 * TrK;
auto Tst05 = std::sqrt( TrK );

// going trough the phase transitions parameters in FtP
// for (unsigned ft = 0; ft < thermo_parameters.phase_transition_prop.size(); ft++)

if (j && thermo_parameters.phase_transition_prop.size() == 0)
{
Exception exception;
exception.error << "No phase transition properties present in the record.";
exception.reason << "For substance "<< substance.symbol() << ".";
exception.line = __LINE__;
RaiseError(exception);
}
if (j == k)
TK = TK_.val /* + C_to_K*/; // current T is the end T for phase transition Cp calculations
else
TK = thermo_parameters.temperature_intervals[j][1] /*+ C_to_K*/; // takes the upper bound from the j-th Tinterval

if ( j && thermo_parameters.phase_transition_prop[ft][0] <= TrK/*-C_to_K*/ )
{ // Adding parameters of phase transition
if ( thermo_parameters.phase_transition_prop[ft].size() > 1 ) // dS
S += thermo_parameters.phase_transition_prop[ft][1];
if ( thermo_parameters.phase_transition_prop[ft].size() > 2 ) // dH
H += thermo_parameters.phase_transition_prop[ft][2];
if ( thermo_parameters.phase_transition_prop[ft].size() > 3 ) // dV
V += thermo_parameters.phase_transition_prop[ft][3];
// More to be added ?
ft++;
}
if (!j)
TrK = substance.referenceT() /*+ C_to_K*/; // if j=0 the first interval should contain the reference T (Tcr)
else
TrK = thermo_parameters.temperature_intervals[j][0] /*+ C_to_K*/; // if j>0 then we are in a different Tinterval and the reference T becomes the lower bound of the interval

G -= S * (TK - TrK);
auto Tst2 = TrK * TrK;
auto Tst3 = Tst2 * TrK;
auto Tst4 = Tst3 * TrK;
auto Tst05 = std::sqrt(TrK);

for (unsigned i=0; i<thermo_parameters.Cp_coeff[j].size(); i++)
{
if (i== 16) break;
ac[i] = thermo_parameters.Cp_coeff[j][i];
}
// going trough the phase transitions parameters in FtP
// for (unsigned ft = 0; ft < thermo_parameters.phase_transition_prop.size(); ft++)

if (j && thermo_parameters.phase_transition_prop.size() == 0)
{
Exception exception;
exception.error << "No phase transition properties present in the record.";
exception.reason << "For substance " << substance.symbol() << ".";
exception.line = __LINE__;
RaiseError(exception)
}

S += ( ac[0] * log( (TK/TrK) ) + ac[1] * (TK - TrK) + ac[2] * ( 1./Tst2 - 1./(TK*TK) ) / 2.
+ ac[3] * 2. * ( 1./Tst05 - 1./(pow(TK,0.5)) ) + ac[4] * ( (TK*TK) - Tst2 ) / 2.
+ ac[5] * ( (pow(TK,3)) - Tst3 ) / 3. + ac[6] * ( (pow(TK,4)) - Tst4 ) / 4.
+ ac[7] * ( 1./ Tst3 - 1./ (pow(TK,3)) ) / 3. + ac[8] * (1. / TrK - 1. / TK )
+ ac[9] * 2.* ( (pow(TK,0.5)) - Tst05 ) );
if (j && thermo_parameters.phase_transition_prop[ft][0] <= TrK /*-C_to_K*/)
{ // Adding parameters of phase transition
if (thermo_parameters.phase_transition_prop[ft].size() > 1) // dS
S += thermo_parameters.phase_transition_prop[ft][1];
if (thermo_parameters.phase_transition_prop[ft].size() > 2) // dH
H += thermo_parameters.phase_transition_prop[ft][2];
if (thermo_parameters.phase_transition_prop[ft].size() > 3) // dV
V += thermo_parameters.phase_transition_prop[ft][3];
// More to be added ?
ft++;
}

G -= ( ac[0] * ( TK * log( (TK/TrK) ) - (TK - TrK) ) + ac[1] * (pow((TK - TrK),2)) / 2. + ac[2] * (pow((TK - TrK),2)) / TK / Tst2 / 2.
+ ac[3] * 2. * ((pow(TK,0.5)) - Tst05)*((pow(TK,0.5)) - Tst05) / Tst05 + ac[4] * ( (pow(TK,3)) + 2.*Tst3 - 3.* TK * Tst2 ) / 6.
+ ac[5] * ( (pow(TK,4)) + 3.*Tst4 - 4.*TK * Tst3 ) / 12. + ac[6] * ( (pow(TK,4))*TK + 4.*Tst4*TrK - 5.*TK*Tst4 ) / 20.
+ ac[7] * ( Tst3 - 3.* (pow(TK,2)) * TrK + 2.*(pow(TK,3)) ) / 6./ (pow(TK,2)) / Tst3 + ac[8] * ( (TK/TrK) - 1. - log( (TK/TrK) ))
+ ac[9] * 2.* ( 2.* TK * (pow(TK,0.5)) - 3.* TK * Tst05 + TrK * Tst05 ) / 3. );
G -= S * (TK - TrK);

H += ( ac[0] * (TK - TrK) + ac[1] * ( (pow(TK,2)) - Tst2 ) / 2. + ac[2] * ( 1./TrK - 1./TK )
+ ac[3] * 2. * ( (pow(TK,0.5)) - Tst05 ) + ac[4] * ( (pow(TK,3)) - Tst3 ) / 3.
+ ac[5] * ( (pow(TK,4)) - Tst4 ) / 4. + ac[6] * ( (pow(TK,4)) * TK - Tst4 * TrK ) / 5
+ ac[7] * ( 1./ Tst2 - 1./ (pow(TK,2)) ) / 2. + ac[8] * log( (TK/TrK) )
+ ac[9] * 2.* ( TK * (pow(TK,0.5)) - TrK * Tst05 ) / 3. );
for (unsigned i = 0; i < thermo_parameters.Cp_coeff[j].size(); i++)
{
if (i == 16)
break;
ac[i] = thermo_parameters.Cp_coeff[j][i];
}

S += (ac[0] * log((TK / TrK)) +
ac[1] * (TK - TrK) +
ac[2] * (1. / Tst2 - 1. / (TK * TK)) / 2. +
ac[3] * 2. * (1. / Tst05 - 1. / (pow(TK, 0.5))) +
ac[4] * ((TK * TK) - Tst2) / 2. +
ac[5] * ((pow(TK, 3)) - Tst3) / 3. +
ac[6] * ((pow(TK, 4)) - Tst4) / 4. +
ac[7] * (1. / Tst3 - 1. / (pow(TK, 3))) / 3. +
ac[8] * (1. / TrK - 1. / TK) +
ac[9] * 2. * ((pow(TK, 0.5)) - Tst05));

G -= (ac[0] * (TK * log((TK / TrK)) - (TK - TrK)) +
ac[1] * (pow((TK - TrK), 2)) / 2. +
ac[2] * (pow((TK - TrK), 2)) / TK / Tst2 / 2. +
ac[3] * 2. * ((pow(TK, 0.5)) - Tst05) * ((pow(TK, 0.5)) - Tst05) / Tst05 +
ac[4] * ((pow(TK, 3)) + 2. * Tst3 - 3. * TK * Tst2) / 6. +
ac[5] * ((pow(TK, 4)) + 3. * Tst4 - 4. * TK * Tst3) / 12. +
ac[6] * ((pow(TK, 4)) * TK + 4. * Tst4 * TrK - 5. * TK * Tst4) / 20. +
ac[7] * (Tst3 - 3. * (pow(TK, 2)) * TrK + 2. * (pow(TK, 3))) / 6. / (pow(TK, 2)) / Tst3 +
ac[8] * ((TK / TrK) - 1. - log((TK / TrK))) +
ac[9] * 2. * (2. * TK * (pow(TK, 0.5)) - 3. * TK * Tst05 + TrK * Tst05) / 3.);

H += (ac[0] * (TK - TrK) +
ac[1] * ((pow(TK, 2)) - Tst2) / 2. +
ac[2] * (1. / TrK - 1. / TK) +
ac[3] * 2. * ((pow(TK, 0.5)) - Tst05) +
ac[4] * ((pow(TK, 3)) - Tst3) / 3. +
ac[5] * ((pow(TK, 4)) - Tst4) / 4. +
ac[6] * ((pow(TK, 4)) * TK - Tst4 * TrK) / 5 +
ac[7] * (1. / Tst2 - 1. / (pow(TK, 2))) / 2. +
ac[8] * log((TK / TrK)) +
ac[9] * 2. * (TK * (pow(TK, 0.5)) - TrK * Tst05) / 3.);
}

thermo_properties_PT.heat_capacity_cp = Cp;
thermo_properties_PT.gibbs_energy = G;
thermo_properties_PT.enthalpy = H;
thermo_properties_PT.entropy = S;
thermo_properties_PT.volume = V;
thermo_properties_PT.heat_capacity_cp = Cp;
thermo_properties_PT.gibbs_energy = G;
thermo_properties_PT.enthalpy = H;
thermo_properties_PT.entropy = S;
thermo_properties_PT.volume = V;

if (k<0)
if (k < 0)
{
setMessage(Reaktoro_::Status::calculated,"Empirical Cp integration: Outside temperature bounds", thermo_properties_PT );
setMessage(Reaktoro_::Status::calculated, "Empirical Cp integration: Outside temperature bounds", thermo_properties_PT);
}

return thermo_properties_PT;
}

}
} // namespace ThermoFun
2 changes: 1 addition & 1 deletion ThermoFun/ThermoEngine.cpp
View file
Original file line number Diff line number Diff line change
@@ -605,7 +605,7 @@ struct ThermoEngine::Impl
tpr.reaction_enthalpy += VP;
auto Vref = pref.workReaction.thermo_ref_prop().reaction_volume;
tpr.log_equilibrium_constant -= Vref *(P_-Pref)/(R_CONSTANT*T)/lg_to_ln;
tpr.reaction_entropy = tpr.reaction_gibbs_energy - T*tpr.reaction_entropy;
tpr.reaction_entropy = (tpr.reaction_enthalpy - tpr.reaction_gibbs_energy)/T;
tpr.reaction_internal_energy = tpr.reaction_enthalpy - P_*tpr.reaction_volume;
tpr.reaction_helmholtz_energy = tpr.reaction_internal_energy - T*tpr.reaction_entropy;

Loading

0 comments on commit e3c2458

Please sign in to comment.