Fit the equilibrium constant values read from the thermodynamic databse at specified temperature values with one of three types of fits: More...

#include <EquilibriumConstantInterpolator.h>

Inheritance diagram for EquilibriumConstantInterpolator:

[legend]

Public Member Functions
	EquilibriumConstantInterpolator (const std::vector< Real > &temperature, const std::vector< Real > &logk, const std::string type, const Real no_value=500.0)

virtual Real	sample (Real T) override

ADReal	sample (const ADReal &T)

Real	sampleDerivative (Real T)
	Sample derivative of function at temperature T. More...

virtual void	generate ()

unsigned int	getSampleSize ()

const std::vector< Real > &	getCoefficients ()

void	setVariables (const std::vector< Real > &x, const std::vector< Real > &y)

Protected Types
enum	FitTypeEnum { FitTypeEnum::FOURTHORDER, FitTypeEnum::MAIERKELLY, FitTypeEnum::LINEAR, FitTypeEnum::PIECEWISELINEAR }
	Functional form of fit. More...

Protected Member Functions
virtual void	fillMatrix () override

void	doLeastSquares ()

Protected Attributes
enum EquilibriumConstantInterpolator::FitTypeEnum	_fit_type

std::unique_ptr< LinearInterpolation >	_linear_interp
	helper object to perform the linear interpolation of the function data More...

std::vector< Real >	_x

std::vector< Real >	_y

std::vector< Real >	_matrix

std::vector< Real >	_coeffs

unsigned int	_num_coeff

Detailed Description

Fit the equilibrium constant values read from the thermodynamic databse at specified temperature values with one of three types of fits:

(1) a fourth-order polynomial fit (for GWB databses)

log(K)= a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4

(2) a Maier-Kelly type fit (for EQ3/6 databases)

log(K)= a_0 ln(T) + a_1 + a_2 T + a_3 / T + a_4 / T^2

(3) a piecewise linear function (for testing purposes) defined by the (temperature, logk) values, with no extrapolation (if T<T_min then the first data value is used, while if T>T_max then the last data value is used)

where T is the temperature in C.

Note: at least five data points must be provided to generate a fit of type (1) or (2), while at least one data point must be provided to generate a piecewise linear function. If fewer data points exist (where 500.0000 represents no value), then a linear fit is used.

The type of fit is specified during construction using the type parameter. The value of "no value" used in the database is provided using the optional no_value parameter. The default value of 500 is used if this parameter is not specified.

Definition at line 43 of file EquilibriumConstantInterpolator.h.

Member Enumeration Documentation

◆ FitTypeEnum

enum EquilibriumConstantInterpolator::FitTypeEnum

strongprotected

Functional form of fit.

Enumerator
FOURTHORDER
MAIERKELLY
LINEAR
PIECEWISELINEAR

Definition at line 65 of file EquilibriumConstantInterpolator.h.

   {
     FOURTHORDER,
     MAIERKELLY,
     LINEAR,
     PIECEWISELINEAR
   } _fit_type;

Constructor & Destructor Documentation

◆ EquilibriumConstantInterpolator()

EquilibriumConstantInterpolator::EquilibriumConstantInterpolator	(	const std::vector< Real > &	temperature,
		const std::vector< Real > &	logk,
		const std::string	type,
		const Real	no_value = `500.0`
	)

Definition at line 16 of file EquilibriumConstantInterpolator.C.

   : LeastSquaresFitBase(), _linear_interp(nullptr)
 {
   // The number of temperature points and logk points must be equal
   if (temperature.size() != logk.size())
     mooseError("The temperature and logk data sets must be equal in length in "
                "EquilibriumConstantInterpolator");
 
   // A value of 500.0 in logk signifies 'no value', so should not be used in the fit
   std::vector<Real> useful_temperature, useful_logk;
   for (unsigned int i = 0; i < temperature.size(); ++i)
     if (!MooseUtils::absoluteFuzzyEqual(logk[i], no_value))
     {
       useful_temperature.push_back(temperature[i]);
       useful_logk.push_back(logk[i]);
     }
 
   LeastSquaresFitBase::setVariables(useful_temperature, useful_logk);
 
   // Set the type of fit
   if (type == "fourth-order")
     _fit_type = FitTypeEnum::FOURTHORDER;
   else if (type == "maier-kelly")
     _fit_type = FitTypeEnum::MAIERKELLY;
   else if (type == "piecewise-linear")
   {
     _fit_type = FitTypeEnum::PIECEWISELINEAR;
     try
     {
       _linear_interp = std::make_unique<LinearInterpolation>(useful_temperature, useful_logk);
     }
     catch (std::domain_error & e)
     {
       mooseError("EquilibriumConstantInterpolation: ", e.what());
     }
   }
   else
     mooseError("Type ", type, " is not supported in EquilibriumConstantInterpolator");
 
   // Set the number of coefficients for the fit
   if (useful_temperature.size() >= 5)
     _num_coeff = 5;
   else if (useful_temperature.size() >= 1 && _fit_type == FitTypeEnum::PIECEWISELINEAR)
     _num_coeff = 2;
   else
   {
     _num_coeff = 2;
     _fit_type = FitTypeEnum::LINEAR;
   }
 
   // If the type is Maier-Kelly and the first temperature point is zero, suggest that
   // the user should use a fourth-order polynomial
   if (MooseUtils::absoluteFuzzyEqual(useful_temperature[0], 0.0) && type == "maier-kelly")
     mooseError("A Maier-Kelly fit cannot be used when the temperature points include 0. Use a "
                "fourth-order fit instead");
 }

Member Function Documentation

◆ fillMatrix()

void EquilibriumConstantInterpolator::fillMatrix ( )

overrideprotectedvirtual

Implements LeastSquaresFitBase.

Definition at line 78 of file EquilibriumConstantInterpolator.C.

 {
   const unsigned int num_rows = _x.size();
   const unsigned int num_cols = _num_coeff;
   _matrix.resize(num_rows * num_cols);
 
   // Type of function depends on fit type
   switch (_fit_type)
   {
     case FitTypeEnum::FOURTHORDER:
     {
       for (unsigned int row = 0; row < num_rows; ++row)
       {
         _matrix[row] = 1.0;
         _matrix[num_rows + row] = _x[row];
         _matrix[(2 * num_rows) + row] = Utility::pow<2>(_x[row]);
         _matrix[(3 * num_rows) + row] = Utility::pow<3>(_x[row]);
         _matrix[(4 * num_rows) + row] = Utility::pow<4>(_x[row]);
       }
       break;
     }
 
     case FitTypeEnum::MAIERKELLY:
     {
       for (unsigned int row = 0; row < num_rows; ++row)
       {
         _matrix[row] = std::log(_x[row]);
         _matrix[num_rows + row] = 1.0;
         _matrix[(2 * num_rows) + row] = _x[row];
         _matrix[(3 * num_rows) + row] = 1.0 / _x[row];
         _matrix[(4 * num_rows) + row] = 1.0 / Utility::pow<2>(_x[row]);
       }
       break;
     }
 
     case FitTypeEnum::PIECEWISELINEAR:
       break; // _matrix is not used in piecewiselinear
 
     default: // FitTypeEnum::LINEAR
     {
       for (unsigned int row = 0; row < num_rows; ++row)
       {
         _matrix[row] = 1.0;
         _matrix[num_rows + row] = _x[row];
       }
       break;
     }
   }
 }

◆ sample() [1/2]

Real EquilibriumConstantInterpolator::sample ( Real T )

overridevirtual

Implements LeastSquaresFitBase.

Definition at line 129 of file EquilibriumConstantInterpolator.C.

Referenced by GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters(), and TEST().

 {
   switch (_fit_type)
   {
     case FitTypeEnum::FOURTHORDER:
       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
              _coeffs[3] * Utility::pow<3>(T) + _coeffs[4] * Utility::pow<4>(T);
 
     case FitTypeEnum::MAIERKELLY:
       return _coeffs[0] * std::log(T) + _coeffs[1] + _coeffs[2] * T + _coeffs[3] / T +
              _coeffs[4] / Utility::pow<2>(T);
 
     case FitTypeEnum::PIECEWISELINEAR:
       return _linear_interp->sample(T);
 
     default: // FitTypeEnum::LINEAR
       return _coeffs[0] + _coeffs[1] * T;
   }
 }

◆ sample() [2/2]

ADReal EquilibriumConstantInterpolator::sample ( const ADReal & T )

Definition at line 171 of file EquilibriumConstantInterpolator.C.

 {
   switch (_fit_type)
   {
     case FitTypeEnum::FOURTHORDER:
       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
              _coeffs[3] * Utility::pow<3>(T) + _coeffs[4] * Utility::pow<4>(T);
 
     case FitTypeEnum::MAIERKELLY:
       return _coeffs[0] * std::log(T) + _coeffs[1] + _coeffs[2] * T + _coeffs[3] / T +
              _coeffs[4] / Utility::pow<2>(T);
 
     case FitTypeEnum::LINEAR:
       return _coeffs[0] + _coeffs[1] * T;
 
     default:
       mooseError("Dual cannot be used for specified fit type in EquilibriumConstantInterpolator");
   }
 }

◆ sampleDerivative()

Real EquilibriumConstantInterpolator::sampleDerivative ( Real T )

Sample derivative of function at temperature T.

Parameters

T	temperature

Returns: derivative of fit wrt T

Definition at line 150 of file EquilibriumConstantInterpolator.C.

Referenced by TEST().

 {
   switch (_fit_type)
   {
     case FitTypeEnum::FOURTHORDER:
       return _coeffs[1] + 2.0 * _coeffs[2] * T + 3.0 * _coeffs[3] * Utility::pow<2>(T) +
              4.0 * _coeffs[4] * Utility::pow<3>(T);
 
     case FitTypeEnum::MAIERKELLY:
       return _coeffs[0] / T + _coeffs[2] - _coeffs[3] / Utility::pow<2>(T) -
              2.0 * _coeffs[4] / Utility::pow<3>(T);
 
     case FitTypeEnum::PIECEWISELINEAR:
       return _linear_interp->sampleDerivative(T);
 
     default: // FitTypeEnum::LINEAR:
       return _coeffs[1];
   }
 }