research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Equil] Eliminate switching between dimensional / nondimensional in VCS

Browse source 
The solver always works in nondimensional units
This commit is contained in:
Ray Speth 2017-08-24 11:08:56 -04:00
parent b92bd24f92
commit 21c2215117
7 changed files with 8 additions and 131 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
include/cantera/equil/vcs_defs.h
View file

@ -77,17 +77,6 @@ namespace Cantera
//@}
/*!
* @name State of Dimensional Units for Gibbs free energies
* @{
*/
//! nondimensional
#define VCS_NONDIMENSIONAL_G 1
//! dimensioned
#define VCS_DIMENSIONAL_G 0
//@}
//! @name Species Categories used during the iteration
/*!
* These defines are valid values for spStatus()

37
include/cantera/equil/vcs_solve.h
View file

@ -445,9 +445,8 @@ public:
/*!
* The actual problem statement is assumed to be in the structure already.
* This is a wrapper around the solve_TP() function. In this wrapper, we
* nondimensionalize the system we calculate the standard state Gibbs free
* energies of the species, and we decide whether to we need to use the
* initial guess algorithm.
* calculate the standard state Gibbs free energies of the species
* and we decide whether to we need to use the initial guess algorithm.
*
* @param ipr = 1 -> Print results to standard output;
* 0 -> don't report on anything
@ -622,31 +621,6 @@ public:
*/
int vcs_report(int iconv);
//! Nondimensionalize the problem data
/*!
* Nondimensionalize the free energies using the divisor, R * T
*
* Essentially the internal data can either be in dimensional form or in
* nondimensional form. This routine switches the data from dimensional form
* into nondimensional form.
*
* @todo Add a scale factor based on the total mole numbers. The algorithm
* contains hard coded numbers based on the total mole number. If we
* ever were faced with a problem with significantly different total
* kmol numbers than one the algorithm would have problems.
*/
void vcs_nondim_TP();
//! Redimensionalize the problem data
/*!
* Redimensionalize the free energies using the multiplier R * T
*
* Essentially the internal data can either be in dimensional form or in
* nondimensional form. This routine switches the data from nondimensional
* form into dimensional form.
*/
void vcs_redim_TP();
//! Computes the current elemental abundances vector
/*!
* Computes the elemental abundances vector, m_elemAbundances[], and stores
@ -1419,13 +1393,6 @@ public:
//! Array of Phase Structures. Length = number of phases.
std::vector<std::unique_ptr<vcs_VolPhase>> m_VolPhaseList;
//! This specifies the current state of units for the Gibbs free energy
//! properties in the program.
/*!
* The default is to have this unitless
*/
char m_unitsState;
//! specifies the activity convention of the phase containing the species
/*!
* * 0 = molar based

12
src/equil/vcs_TP.cpp
View file

@ -13,15 +13,12 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
// Store the temperature and pressure in the private global variables
m_temperature = T_arg;
m_pressurePA = pres_arg;
m_Faraday_dim = m_Faraday_dim = Faraday / (m_temperature * GasConstant);
// Evaluate the standard state free energies
// at the current temperatures and pressures.
int iconv = vcs_evalSS_TP(ipr, ip1, m_temperature, pres_arg);
// Prepare the problem data: nondimensionalize the free energies using the
// divisor, R * T
vcs_nondim_TP();
// Prep the fe field
vcs_fePrep_TP();
@ -39,10 +36,6 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
// energies are now in dimensionless form.)
iconv = vcs_solve_TP(ipr, ip1, maxit);
// Redimensionalize the free energies using the reverse of vcs_nondim to add
// back units.
vcs_redim_TP();
// Return the convergence success flag.
return iconv;
}
@ -54,6 +47,9 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
vph->setState_TP(m_temperature, m_pressurePA);
vph->sendToVCS_GStar(&m_SSfeSpecies[0]);
}
for (size_t k = 0; k < m_nsp; k++) {
m_SSfeSpecies[k] /= GasConstant * m_temperature;
}
return VCS_SUCCESS;
}

54
src/equil/vcs_nondim.cpp
View file

@ -1,54 +0,0 @@
/**
* @file vcs_nondim.cpp
* Nondimensionalization routines within VCSnonideal
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at http://www.cantera.org/license.txt for license and copyright information.
#include "cantera/equil/vcs_solve.h"
#include "cantera/equil/vcs_VolPhase.h"
#include "cantera/base/stringUtils.h"
#include "cantera/base/ctexceptions.h"
namespace Cantera
{
void VCS_SOLVE::vcs_nondim_TP()
{
if (m_unitsState == VCS_DIMENSIONAL_G) {
m_unitsState = VCS_NONDIMENSIONAL_G;
double tf = 1.0 / (GasConstant * m_temperature);
for (size_t i = 0; i < m_nsp; ++i) {
// Modify the standard state and total chemical potential data,
// FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may
// divide it by the temperature.
m_SSfeSpecies[i] *= tf;
m_deltaGRxn_new[i] *= tf;
m_deltaGRxn_old[i] *= tf;
m_feSpecies_old[i] *= tf;
}
m_Faraday_dim = ElectronCharge * Avogadro / (m_temperature * GasConstant);
}
}
void VCS_SOLVE::vcs_redim_TP()
{
if (m_unitsState != VCS_DIMENSIONAL_G) {
m_unitsState = VCS_DIMENSIONAL_G;
double tf = m_temperature * GasConstant;
for (size_t i = 0; i < m_nsp; ++i) {
// Modify the standard state and total chemical potential data,
// FF(I), to make it have units, i.e. mu = RT * mu_star
m_SSfeSpecies[i] *= tf;
m_deltaGRxn_new[i] *= tf;
m_deltaGRxn_old[i] *= tf;
m_feSpecies_old[i] *= tf;
}
m_Faraday_dim *= tf;
}
}
}

18
src/equil/vcs_report.cpp
View file

@ -12,7 +12,6 @@ namespace Cantera
int VCS_SOLVE::vcs_report(int iconv)
{
bool inertYes = false;
char originalUnitsState = m_unitsState;
std::vector<size_t> sortindex(m_nsp, 0);
vector_fp xy(m_nsp, 0.0);
@ -32,13 +31,6 @@ int VCS_SOLVE::vcs_report(int iconv)
}
}
// Decide whether we have to nondimensionalize the equations. For the
// printouts from this routine, we will use nondimensional representations.
// This may be expanded in the future.
if (m_unitsState == VCS_DIMENSIONAL_G) {
vcs_nondim_TP();
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp);
@ -309,16 +301,6 @@ int VCS_SOLVE::vcs_report(int iconv)
writeline('-', 80);
writeline('-', 80);
// Set the Units state of the system back to where it was when we
// entered the program.
if (originalUnitsState != m_unitsState) {
if (originalUnitsState == VCS_DIMENSIONAL_G) {
vcs_redim_TP();
} else {
vcs_nondim_TP();
}
}
// Return a successful completion flag
return VCS_SUCCESS;
}

4
src/equil/vcs_solve.cpp
View file

@ -40,10 +40,9 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_tolmin(1.0E-6),
m_tolmaj2(1.0E-10),
m_tolmin2(1.0E-8),
m_unitsState(VCS_DIMENSIONAL_G),
m_useActCoeffJac(0),
m_totalVol(mphase->volume()),
m_Faraday_dim(ElectronCharge * Avogadro),
m_Faraday_dim(Faraday / (m_temperature * GasConstant)),
m_VCount(0),
m_debug_print_lvl(0),
m_timing_print_lvl(1)
@ -543,6 +542,7 @@ void VCS_SOLVE::vcs_prob_specifyFully()
// the call to the equilibrium solver.
m_temperature = m_mix->temperature();
m_pressurePA = m_mix->pressure();
m_Faraday_dim = Faraday / (m_temperature * GasConstant);
m_totalVol = m_mix->volume();
vector<size_t> invSpecies(m_nsp);

3
src/equil/vcs_solve_TP.cpp
View file

@ -2733,9 +2733,6 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
"Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc);
}
AssertThrowMsg(m_unitsState != VCS_DIMENSIONAL_G, "VCS_SOLVE::vcs_dfe",
"called with wrong units state");
if (m_debug_print_lvl >= 2) {
if (ll == 0) {
if (lbot != 0) {