Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "EquilibriumConstantInterpolator.h"
      11             : #include "MooseError.h"
      12             : #include "MooseUtils.h"
      13             : #include "libmesh/utility.h"
      14             : #include <memory> // std::make_unique
      15             : 
      16       63739 : EquilibriumConstantInterpolator::EquilibriumConstantInterpolator(
      17             :     const std::vector<Real> & temperature,
      18             :     const std::vector<Real> & logk,
      19             :     const std::string type,
      20       63739 :     const Real no_value)
      21       63739 :   : LeastSquaresFitBase(), _linear_interp(nullptr)
      22             : {
      23             :   // The number of temperature points and logk points must be equal
      24       63739 :   if (temperature.size() != logk.size())
      25           1 :     mooseError("The temperature and logk data sets must be equal in length in "
      26             :                "EquilibriumConstantInterpolator");
      27             : 
      28             :   // A value of 500.0 in logk signifies 'no value', so should not be used in the fit
      29             :   std::vector<Real> useful_temperature, useful_logk;
      30      573598 :   for (unsigned int i = 0; i < temperature.size(); ++i)
      31      509860 :     if (!MooseUtils::absoluteFuzzyEqual(logk[i], no_value))
      32             :     {
      33      505856 :       useful_temperature.push_back(temperature[i]);
      34      505856 :       useful_logk.push_back(logk[i]);
      35             :     }
      36             : 
      37       63738 :   LeastSquaresFitBase::setVariables(useful_temperature, useful_logk);
      38             : 
      39             :   // Set the type of fit
      40       63738 :   if (type == "fourth-order")
      41       24880 :     _fit_type = FitTypeEnum::FOURTHORDER;
      42       38858 :   else if (type == "maier-kelly")
      43           4 :     _fit_type = FitTypeEnum::MAIERKELLY;
      44       38854 :   else if (type == "piecewise-linear")
      45             :   {
      46       38853 :     _fit_type = FitTypeEnum::PIECEWISELINEAR;
      47             :     try
      48             :     {
      49       77705 :       _linear_interp = std::make_unique<LinearInterpolation>(useful_temperature, useful_logk);
      50             :     }
      51           1 :     catch (std::domain_error & e)
      52             :     {
      53           1 :       mooseError("EquilibriumConstantInterpolation: ", e.what());
      54           1 :     }
      55             :   }
      56             :   else
      57           1 :     mooseError("Type ", type, " is not supported in EquilibriumConstantInterpolator");
      58             : 
      59             :   // Set the number of coefficients for the fit
      60       63736 :   if (useful_temperature.size() >= 5)
      61       62957 :     _num_coeff = 5;
      62         779 :   else if (useful_temperature.size() >= 1 && _fit_type == FitTypeEnum::PIECEWISELINEAR)
      63          63 :     _num_coeff = 2;
      64             :   else
      65             :   {
      66         716 :     _num_coeff = 2;
      67         716 :     _fit_type = FitTypeEnum::LINEAR;
      68             :   }
      69             : 
      70             :   // If the type is Maier-Kelly and the first temperature point is zero, suggest that
      71             :   // the user should use a fourth-order polynomial
      72       62959 :   if (MooseUtils::absoluteFuzzyEqual(useful_temperature[0], 0.0) && type == "maier-kelly")
      73           1 :     mooseError("A Maier-Kelly fit cannot be used when the temperature points include 0. Use a "
      74             :                "fourth-order fit instead");
      75       63749 : }
      76             : 
      77             : void
      78       63730 : EquilibriumConstantInterpolator::fillMatrix()
      79             : {
      80       63730 :   const unsigned int num_rows = _x.size();
      81       63730 :   const unsigned int num_cols = _num_coeff;
      82       63730 :   _matrix.resize(num_rows * num_cols);
      83             : 
      84             :   // Type of function depends on fit type
      85       63730 :   switch (_fit_type)
      86             :   {
      87             :     case FitTypeEnum::FOURTHORDER:
      88             :     {
      89      217465 :       for (unsigned int row = 0; row < num_rows; ++row)
      90             :       {
      91      193300 :         _matrix[row] = 1.0;
      92      193300 :         _matrix[num_rows + row] = _x[row];
      93      193300 :         _matrix[(2 * num_rows) + row] = Utility::pow<2>(_x[row]);
      94      193300 :         _matrix[(3 * num_rows) + row] = Utility::pow<3>(_x[row]);
      95      193300 :         _matrix[(4 * num_rows) + row] = Utility::pow<4>(_x[row]);
      96             :       }
      97             :       break;
      98             :     }
      99             : 
     100             :     case FitTypeEnum::MAIERKELLY:
     101             :     {
     102           9 :       for (unsigned int row = 0; row < num_rows; ++row)
     103             :       {
     104           8 :         _matrix[row] = std::log(_x[row]);
     105           8 :         _matrix[num_rows + row] = 1.0;
     106           8 :         _matrix[(2 * num_rows) + row] = _x[row];
     107           8 :         _matrix[(3 * num_rows) + row] = 1.0 / _x[row];
     108           8 :         _matrix[(4 * num_rows) + row] = 1.0 / Utility::pow<2>(_x[row]);
     109             :       }
     110             :       break;
     111             :     }
     112             : 
     113             :     case FitTypeEnum::PIECEWISELINEAR:
     114             :       break; // _matrix is not used in piecewiselinear
     115             : 
     116             :     default: // FitTypeEnum::LINEAR
     117             :     {
     118        2857 :       for (unsigned int row = 0; row < num_rows; ++row)
     119             :       {
     120        2142 :         _matrix[row] = 1.0;
     121        2142 :         _matrix[num_rows + row] = _x[row];
     122             :       }
     123             :       break;
     124             :     }
     125             :   }
     126       63730 : }
     127             : 
     128             : Real
     129      523598 : EquilibriumConstantInterpolator::sample(Real T)
     130             : {
     131      523598 :   switch (_fit_type)
     132             :   {
     133      316173 :     case FitTypeEnum::FOURTHORDER:
     134      316173 :       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
     135      316173 :              _coeffs[3] * Utility::pow<3>(T) + _coeffs[4] * Utility::pow<4>(T);
     136             : 
     137           6 :     case FitTypeEnum::MAIERKELLY:
     138           6 :       return _coeffs[0] * std::log(T) + _coeffs[1] + _coeffs[2] * T + _coeffs[3] / T +
     139           6 :              _coeffs[4] / Utility::pow<2>(T);
     140             : 
     141      203116 :     case FitTypeEnum::PIECEWISELINEAR:
     142      203116 :       return _linear_interp->sample(T);
     143             : 
     144        4303 :     default: // FitTypeEnum::LINEAR
     145        4303 :       return _coeffs[0] + _coeffs[1] * T;
     146             :   }
     147             : }
     148             : 
     149             : Real
     150          36 : EquilibriumConstantInterpolator::sampleDerivative(Real T)
     151             : {
     152          36 :   switch (_fit_type)
     153             :   {
     154          21 :     case FitTypeEnum::FOURTHORDER:
     155          21 :       return _coeffs[1] + 2.0 * _coeffs[2] * T + 3.0 * _coeffs[3] * Utility::pow<2>(T) +
     156          21 :              4.0 * _coeffs[4] * Utility::pow<3>(T);
     157             : 
     158           1 :     case FitTypeEnum::MAIERKELLY:
     159           1 :       return _coeffs[0] / T + _coeffs[2] - _coeffs[3] / Utility::pow<2>(T) -
     160           1 :              2.0 * _coeffs[4] / Utility::pow<3>(T);
     161             : 
     162          10 :     case FitTypeEnum::PIECEWISELINEAR:
     163          10 :       return _linear_interp->sampleDerivative(T);
     164             : 
     165           4 :     default: // FitTypeEnum::LINEAR:
     166           4 :       return _coeffs[1];
     167             :   }
     168             : }
     169             : 
     170             : ADReal
     171           7 : EquilibriumConstantInterpolator::sample(const ADReal & T)
     172             : {
     173             :   using std::log;
     174           7 :   switch (_fit_type)
     175             :   {
     176           3 :     case FitTypeEnum::FOURTHORDER:
     177           3 :       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
     178           3 :              _coeffs[3] * Utility::pow<3>(T) + _coeffs[4] * Utility::pow<4>(T);
     179             : 
     180           1 :     case FitTypeEnum::MAIERKELLY:
     181           1 :       return _coeffs[0] * log(T) + _coeffs[1] + _coeffs[2] * T + _coeffs[3] / T +
     182           1 :              _coeffs[4] / Utility::pow<2>(T);
     183             : 
     184           2 :     case FitTypeEnum::LINEAR:
     185           2 :       return _coeffs[0] + _coeffs[1] * T;
     186             : 
     187           1 :     default:
     188           1 :       mooseError("Dual cannot be used for specified fit type in EquilibriumConstantInterpolator");
     189             :   }
     190             : }