Skip to content

Commit

Permalink
Fixes for formula parser and CpfT for phase transition (#24)
Browse files Browse the repository at this point in the history 
* fix for Cp=f(T) phase transition, problems with the derivatives

* added default valences and fixed formula parser bug

* tests for formula parser and CpfT with phase transition
  • Loading branch information
@gdmiron
gdmiron committed Jan 6, 2021
1 parent 5a781ae commit 626ca33
Show file tree
Hide file tree
Showing 10 changed files with 1,356 additions and 5 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
cmake_minimum_required(VERSION 3.9)

# Set the name of the project
project(ThermoFun VERSION 0.3.5 LANGUAGES CXX)
project(ThermoFun VERSION 0.3.6 LANGUAGES CXX)

# Set the cmake module path of the project
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules")
Expand Down
14 changes: 12 additions & 2 deletions ThermoFun/Common/formuladata.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,10 @@ void FormulaToken::unpack( list<ICTERM>& itt_ )
}
datamap.push_back( FormulaValues( key, itr->stoich, itr->val ));
elements.insert(key);
elements_map.insert(pair<ElementKey,double>(key,itr->stoich));
if (elements_map.find(key) != elements_map.end())
elements_map.at(key) += itr->stoich;
else
elements_map.insert(pair<ElementKey,double>(key,itr->stoich));
// coefficients.push_back(itr->stoich);
itr++;
}
Expand Down Expand Up @@ -407,9 +410,16 @@ void ChemicalFormula::addOneElement(Element e)
eldata.entropy = e.entropy();
eldata.heat_capacity = e.heatCapacity();
eldata.volume = e.volume();
if (e.valence()==777)
{
if (elements_valences.find(e.symbol()) != elements_valences.end())
e.setValence(elements_valences.at(e.symbol()));
else
e.setValence(0);
}
eldata.valence = e.valence();
eldata.number = e.number();
eldata.name = e.name();
eldata.name = e.symbol(); // was e.name();

dbElements[elkey] = eldata;
}
Expand Down
105 changes: 105 additions & 0 deletions ThermoFun/Common/formuladata.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,111 @@ namespace ThermoFun {

class Element;

const std::map<std::string, int> elements_valences = {
{"Ac", 3},
{"Ag", 1},
{"Al", 3},
{"Ar", 0},
{"Am", 3},
{"As", 5},
{"Au", 1},
{"B", 3},
{"Ba", 2},
{"Be", 2},
{"Bi", 3},
{"Br", -1},
{"C", 4},
{"Ca", 2},
{"Cd", 2},
{"Ce", 3},
{"Cf", 3},
{"Cit", -3},
{"Cl", -1},
{"Co", 2},
{"Cr", 3},
{"Cm", 3},
{"Cs", 1},
{"Cu", 2},
{"Dy", 3},
{"Edta", -4},
{"Er", 3},
{"Eu", 3},
{"F", -1},
{"Fr", 1},
{"Fe", 2},
{"Ga", 3},
{"Gd", 3},
{"Ge", 4},
{"H", 1},
{"He", 0},
{"Hf", 4},
{"Hg", 2},
{"Ho", 3},
{"I", -1},
{"In", 3},
{"Isa", -4},
{"Ir", 4},
{"K", 1},
{"Kr", 0},
{"La", 3},
{"Li", 1},
{"Lu", 3},
{"Mg", 2},
{"Mn", 2},
{"Mo", 6},
{"N", 5},
{"N_atm", 0},
{"Na", 1},
{"Nb", 5},
{"Nd", 3},
{"Ne", 0},
{"Ni", 2},
{"Np", 6},
{"O", -2},
{"Os", 4},
{"Ox", -2},
{"P", 5},
{"Pa", 5},
{"Pb", 2},
{"Pd", 2},
{"Po", 4},
{"Pu", 6},
{"Pr", 3},
{"Pm", 3},
{"Pt", 2},
{"Ra", 2},
{"Rb", 1},
{"Re", 4},
{"Rh", 2},
{"Rn", 0},
{"Ru", 2},
{"S", 6},
{"Sb", 3},
{"Sc", 3},
{"Se", 4},
{"Si", 4},
{"Sm", 3},
{"Sn", 2},
{"Sr", 2},
{"Ta", 5},
{"Tb", 3},
{"Tc", 7},
{"Te", 6},
{"Th", 4},
{"Ti", 4},
{"Tl", 1},
{"Tm", 3},
{"U", 6},
{"V", 5},
{"W", 6},
{"Xe", 0},
{"Y", 3},
{"Yb", 3},
{"Zn", 2},
{"Zr", 4},
{"Zz", 0}
};


/// Key fields of Element vertex
class ElementKey
Expand Down
3 changes: 3 additions & 0 deletions ThermoFun/Database.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -485,6 +485,9 @@ auto Database::parseSubstanceFormula(std::string formula_) const -> std::map<Ele

formula.setFormula(formula_);

// FormulaProperites props;
// formula.calcFormulaProperites(props);

for (auto element : formula.getElements_map())
{
Element e = elementKeyToElement(element.first);
Expand Down
2 changes: 1 addition & 1 deletion ThermoFun/Element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ struct Element::Impl
double volume;

/// The valence of the element
int valence;
int valence = 777;

/// The molar mass of the chemical element (in units of g/mol)
double molarMass;
Expand Down
79 changes: 78 additions & 1 deletion ThermoFun/Substances/EmpiricalCpIntegration.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,10 @@ auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pres

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

auto Sr = thermo_properties_PrTr.entropy;
auto Gr = thermo_properties_PrTr.gibbs_energy;
auto Hr = thermo_properties_PrTr.enthalpy;

auto S = thermo_properties_PrTr.entropy;
auto G = thermo_properties_PrTr.gibbs_energy;
auto H = thermo_properties_PrTr.enthalpy;
Expand Down Expand Up @@ -65,7 +69,7 @@ auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pres
<< __LINE__;
}

k = 0;
//k = 0; fix

for (unsigned i = 0; i < thermo_parameters.Cp_coeff[k].size(); i++)
{
Expand Down Expand Up @@ -180,6 +184,79 @@ auto thermoPropertiesEmpCpIntegration(Reaktoro_::Temperature TK, Reaktoro_::Pres
setMessage(Reaktoro_::Status::calculated, "Empirical Cp integration: Outside temperature bounds", thermo_properties_PT);
}


/// reaktoro implementation
/*
// Collect the temperature points used for the integrals along the pressure line P = Pr
std::vector<Reaktoro_::Temperature> Ti;
const auto& Tr = substance.referenceT();
std::vector<double> dHt;
std::vector<double> dVt;
Ti.push_back(substance.referenceT());
for(int i = 0; i < thermo_parameters.phase_transition_prop.size(); ++i)
{
if(TK_ > thermo_parameters.temperature_intervals[i][1])
{Ti.push_back(thermo_parameters.temperature_intervals[i][1]);
dHt.push_back(thermo_parameters.phase_transition_prop[i][2]);
dVt.push_back(thermo_parameters.phase_transition_prop[i][3]);}
}
Ti.push_back(TK_);
Reaktoro_::ThermoScalar xCp;
for(unsigned i = 0; i+1 < Ti.size(); ++i)
if(Ti[i] <= TK_ && TK_ <= Ti[i+1])
xCp = thermo_parameters.Cp_coeff[i][0] + thermo_parameters.Cp_coeff[i][1]*TK_ + thermo_parameters.Cp_coeff[i][2]/(TK_*TK_);
// Calculate the integrals of the heat capacity function of the mineral from Tr to T at constant pressure Pr
Reaktoro_::ThermoScalar CpdT;
Reaktoro_::ThermoScalar CpdlnT;
for(unsigned i = 0; i+1 < Ti.size(); ++i)
{
const auto T0 = Ti[i];
const auto T1 = Ti[i+1];
for (unsigned j = 0; j < thermo_parameters.Cp_coeff[i].size(); j++)
{
if (j == 16)
break;
ac[j] = thermo_parameters.Cp_coeff[i][j];
}
CpdT += ac[0]*(T1 - T0) + 0.5*ac[1]*(T1*T1 - T0*T0) - ac[2]*(1.0/T1 - 1.0/T0);
CpdlnT += ac[0]*log(T1/T0) + ac[1]*(T1 - T0) - 0.5*ac[2]*(1.0/(T1*T1) - 1.0/(T0*T0));
}
// Calculate the volume and other auxiliary quantities for the thermodynamic properties of the mineral
Reaktoro_::ThermoScalar xV(0.0);
Reaktoro_::ThermoScalar GdH;
Reaktoro_::ThermoScalar HdH;
Reaktoro_::ThermoScalar SdH;
for(unsigned i = 1; i+1 < Ti.size(); ++i)
{
GdH += dHt[i-1]*(TK_ - Ti[i])/Ti[i];
HdH += dHt[i-1];
SdH += dHt[i-1]/Ti[i];
xV += dVt[i-1];
}
// Calculate the standard molal thermodynamic properties of the mineral
auto xG = Gr - Sr * (TK_ - Tr) + CpdT - TK_ * CpdlnT - GdH; // + VdP
auto xH = Hr + CpdT + HdH; // + VdP
auto xS = Sr + CpdlnT + SdH;
// auto xU = xH - Pb*V;
// auto xA = xU - TK_*S;
*/

return thermo_properties_PT;
}

Expand Down
Loading

0 comments on commit 626ca33

Please sign in to comment.