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 "VariableTimeIntegrationAux.h"
      11             : 
      12             : registerMooseObject("MooseApp", VariableTimeIntegrationAux);
      13             : 
      14             : InputParameters
      15       14340 : VariableTimeIntegrationAux::validParams()
      16             : {
      17       14340 :   InputParameters params = AuxKernel::validParams();
      18       14340 :   params.addClassDescription("Integrates a field variable in time.");
      19       14340 :   params.addRequiredCoupledVar("variable_to_integrate", "The variable to be integrated");
      20       14340 :   params.addParam<Real>("coefficient", 1.0, "A simple coefficient multiplying the integral");
      21       43020 :   params.addParam<unsigned int>(
      22       28680 :       "order", 2, "The order of global truncation error: midpoint=1, trapezoidal=2, Simpson=3");
      23       14340 :   return params;
      24           0 : }
      25             : 
      26          39 : VariableTimeIntegrationAux::VariableTimeIntegrationAux(const InputParameters & parameters)
      27             :   : AuxKernel(parameters),
      28          39 :     _coef(getParam<Real>("coefficient")),
      29          39 :     _order(getParam<unsigned int>("order")),
      30          39 :     _u_old(_order != 3 ? uOld() : genericZeroValue<false>()),
      31         117 :     _u_older(_order == 3 ? uOlder() : genericZeroValue<false>())
      32             : {
      33          39 :   switch (_order)
      34             :   {
      35          13 :     case 1:
      36          13 :       _integration_coef.push_back(1.0);
      37          13 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      38          13 :       break;
      39          13 :     case 2:
      40          13 :       _integration_coef.push_back(0.5);
      41          13 :       _integration_coef.push_back(0.5);
      42          13 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      43          13 :       _coupled_vars.push_back(&coupledValueOld("variable_to_integrate"));
      44          13 :       break;
      45          13 :     case 3:
      46          13 :       _integration_coef.push_back(1.0 / 3.0);
      47          13 :       _integration_coef.push_back(4.0 / 3.0);
      48          13 :       _integration_coef.push_back(1.0 / 3.0);
      49          13 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      50          13 :       _coupled_vars.push_back(&coupledValueOld("variable_to_integrate"));
      51          13 :       _coupled_vars.push_back(&coupledValueOlder("variable_to_integrate"));
      52          13 :       break;
      53           0 :     default:
      54           0 :       mooseError("VariableTimeIntegrationAux: unknown time integration order specified");
      55             :   }
      56          39 : }
      57             : 
      58             : Real
      59      748506 : VariableTimeIntegrationAux::computeValue()
      60             : {
      61      748506 :   Real integral = getIntegralValue();
      62             : 
      63      748506 :   if (_order == 3)
      64      249502 :     return _u_older[_qp] + _coef * integral;
      65             : 
      66      499004 :   return _u_old[_qp] + _coef * integral;
      67             : }
      68             : 
      69             : Real
      70      748506 : VariableTimeIntegrationAux::getIntegralValue()
      71             : {
      72      748506 :   Real integral_value = 0.0;
      73             : 
      74     2245518 :   for (unsigned int i = 0; i < _order; ++i)
      75     1497012 :     integral_value += _integration_coef[i] * (*_coupled_vars[i])[_qp] * _dt;
      76             : 
      77             :   /**
      78             :    * Subsequent timesteps may be unequal, so the standard Simpson rule
      79             :    * cannot be used. Use a different set of coefficients here.
      80             :    * J. McNAMEE, "A PROGRAM TO INTEGRATE A FUNCTION TABULATED AT
      81             :    * UNEQUAL INTERVALS," Internation Journal for Numerical Methods in
      82             :    * Engineering, Vol. 17, 217-279. (1981).
      83             :    */
      84      748506 :   if (_order == 3 && _dt != _dt_old)
      85             :   {
      86      171094 :     Real x0 = 0.0;
      87      171094 :     Real x1 = _dt_old;
      88      171094 :     Real x2 = _dt + _dt_old;
      89      171094 :     Real y0 = (*_coupled_vars[2])[_qp];
      90      171094 :     Real y1 = (*_coupled_vars[1])[_qp];
      91      171094 :     Real y2 = (*_coupled_vars[0])[_qp];
      92      171094 :     Real term1 = (x2 - x0) * (y0 + (x2 - x0) * (y1 - y0) / (2.0 * (x1 - x0)));
      93      171094 :     Real term2 = (2.0 * x2 * x2 - x0 * x2 - x0 * x0 + 3.0 * x0 * x1 - 3.0 * x1 * x2) / 6.0;
      94      171094 :     Real term3 = (y2 - y1) / (x2 - x1) - (y1 - y0) / (x1 - x0);
      95      171094 :     integral_value = term1 + term2 * term3;
      96             :   }
      97             : 
      98      748506 :   return integral_value;
      99             : }