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       79158 : EquilibriumConstantInterpolator::EquilibriumConstantInterpolator(
      17             :     const std::vector<Real> & temperature,
      18             :     const std::vector<Real> & logk,
      19             :     const std::string type,
      20       79158 :     const Real no_value)
      21       79158 :   : LeastSquaresFitBase(), _linear_interp(nullptr)
      22             : {
      23             :   // The number of temperature points and logk points must be equal
      24       79158 :   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      712369 :   for (unsigned int i = 0; i < temperature.size(); ++i)
      31      633212 :     if (!MooseUtils::absoluteFuzzyEqual(logk[i], no_value))
      32             :     {
      33      628038 :       useful_temperature.push_back(temperature[i]);
      34      628038 :       useful_logk.push_back(logk[i]);
      35             :     }
      36             : 
      37       79157 :   LeastSquaresFitBase::setVariables(useful_temperature, useful_logk);
      38             : 
      39             :   // Set the type of fit
      40       79157 :   if (type == "fourth-order")
      41       30326 :     _fit_type = FitTypeEnum::FOURTHORDER;
      42       48831 :   else if (type == "maier-kelly")
      43           4 :     _fit_type = FitTypeEnum::MAIERKELLY;
      44       48827 :   else if (type == "piecewise-linear")
      45             :   {
      46       48826 :     _fit_type = FitTypeEnum::PIECEWISELINEAR;
      47             :     try
      48             :     {
      49       97651 :       _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       79155 :   if (useful_temperature.size() >= 5)
      61       78152 :     _num_coeff = 5;
      62        1003 :   else if (useful_temperature.size() >= 1 && _fit_type == FitTypeEnum::PIECEWISELINEAR)
      63          91 :     _num_coeff = 2;
      64             :   else
      65             :   {
      66         912 :     _num_coeff = 2;
      67         912 :     _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       78154 :   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       79162 : }
      76             : 
      77             : void
      78       79149 : EquilibriumConstantInterpolator::fillMatrix()
      79             : {
      80       79149 :   const unsigned int num_rows = _x.size();
      81       79149 :   const unsigned int num_cols = _num_coeff;
      82       79149 :   _matrix.resize(num_rows * num_cols);
      83             : 
      84             :   // Type of function depends on fit type
      85       79149 :   switch (_fit_type)
      86             :   {
      87             :     case FitTypeEnum::FOURTHORDER:
      88             :     {
      89      264709 :       for (unsigned int row = 0; row < num_rows; ++row)
      90             :       {
      91      235294 :         _matrix[row] = 1.0;
      92      235294 :         _matrix[num_rows + row] = _x[row];
      93      235294 :         _matrix[(2 * num_rows) + row] = Utility::pow<2>(_x[row]);
      94      235294 :         _matrix[(3 * num_rows) + row] = Utility::pow<3>(_x[row]);
      95      235294 :         _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        3641 :       for (unsigned int row = 0; row < num_rows; ++row)
     119             :       {
     120        2730 :         _matrix[row] = 1.0;
     121        2730 :         _matrix[num_rows + row] = _x[row];
     122             :       }
     123             :       break;
     124             :     }
     125             :   }
     126       79149 : }
     127             : 
     128             : Real
     129      598938 : EquilibriumConstantInterpolator::sample(Real T)
     130             : {
     131      598938 :   switch (_fit_type)
     132             :   {
     133      366338 :     case FitTypeEnum::FOURTHORDER:
     134      366338 :       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
     135      366338 :              _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      227499 :     case FitTypeEnum::PIECEWISELINEAR:
     142      227499 :       return _linear_interp->sample(T);
     143             : 
     144        5095 :     default: // FitTypeEnum::LINEAR
     145        5095 :       return _coeffs[0] + _coeffs[1] * T;
     146             :   }
     147             : }
     148             : 
     149             : Real
     150          42 : EquilibriumConstantInterpolator::sampleDerivative(Real T)
     151             : {
     152          42 :   switch (_fit_type)
     153             :   {
     154          27 :     case FitTypeEnum::FOURTHORDER:
     155          27 :       return _coeffs[1] + 2.0 * _coeffs[2] * T + 3.0 * _coeffs[3] * Utility::pow<2>(T) +
     156          27 :              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           7 :   switch (_fit_type)
     174             :   {
     175           3 :     case FitTypeEnum::FOURTHORDER:
     176           3 :       return _coeffs[0] + _coeffs[1] * T + _coeffs[2] * Utility::pow<2>(T) +
     177           3 :              _coeffs[3] * Utility::pow<3>(T) + _coeffs[4] * Utility::pow<4>(T);
     178             : 
     179           1 :     case FitTypeEnum::MAIERKELLY:
     180           1 :       return _coeffs[0] * std::log(T) + _coeffs[1] + _coeffs[2] * T + _coeffs[3] / T +
     181           1 :              _coeffs[4] / Utility::pow<2>(T);
     182             : 
     183           2 :     case FitTypeEnum::LINEAR:
     184           2 :       return _coeffs[0] + _coeffs[1] * T;
     185             : 
     186           1 :     default:
     187           1 :       mooseError("Dual cannot be used for specified fit type in EquilibriumConstantInterpolator");
     188             :   }
     189             : }