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       14346 : VariableTimeIntegrationAux::validParams()
      16             : {
      17       14346 :   InputParameters params = AuxKernel::validParams();
      18       14346 :   params.addClassDescription("Integrates a field variable in time.");
      19       14346 :   params.addRequiredCoupledVar("variable_to_integrate", "The variable to be integrated");
      20       14346 :   params.addParam<Real>("coefficient", 1.0, "A simple coefficient multiplying the integral");
      21       43038 :   params.addParam<unsigned int>(
      22       28692 :       "order", 2, "The order of global truncation error: midpoint=1, trapezoidal=2, Simpson=3");
      23       14346 :   return params;
      24           0 : }
      25             : 
      26          42 : VariableTimeIntegrationAux::VariableTimeIntegrationAux(const InputParameters & parameters)
      27             :   : AuxKernel(parameters),
      28          42 :     _coef(getParam<Real>("coefficient")),
      29          42 :     _order(getParam<unsigned int>("order")),
      30          42 :     _u_old(_order != 3 ? uOld() : genericZeroValue<false>()),
      31         126 :     _u_older(_order == 3 ? uOlder() : genericZeroValue<false>())
      32             : {
      33          42 :   switch (_order)
      34             :   {
      35          14 :     case 1:
      36          14 :       _integration_coef.push_back(1.0);
      37          14 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      38          14 :       break;
      39          14 :     case 2:
      40          14 :       _integration_coef.push_back(0.5);
      41          14 :       _integration_coef.push_back(0.5);
      42          14 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      43          14 :       _coupled_vars.push_back(&coupledValueOld("variable_to_integrate"));
      44          14 :       break;
      45          14 :     case 3:
      46          14 :       _integration_coef.push_back(1.0 / 3.0);
      47          14 :       _integration_coef.push_back(4.0 / 3.0);
      48          14 :       _integration_coef.push_back(1.0 / 3.0);
      49          14 :       _coupled_vars.push_back(&coupledValue("variable_to_integrate"));
      50          14 :       _coupled_vars.push_back(&coupledValueOld("variable_to_integrate"));
      51          14 :       _coupled_vars.push_back(&coupledValueOlder("variable_to_integrate"));
      52          14 :       break;
      53           0 :     default:
      54           0 :       mooseError("VariableTimeIntegrationAux: unknown time integration order specified");
      55             :   }
      56          42 : }
      57             : 
      58             : Real
      59      849420 : VariableTimeIntegrationAux::computeValue()
      60             : {
      61      849420 :   Real integral = getIntegralValue();
      62             : 
      63      849420 :   if (_order == 3)
      64      283140 :     return _u_older[_qp] + _coef * integral;
      65             : 
      66      566280 :   return _u_old[_qp] + _coef * integral;
      67             : }
      68             : 
      69             : Real
      70      849420 : VariableTimeIntegrationAux::getIntegralValue()
      71             : {
      72      849420 :   Real integral_value = 0.0;
      73             : 
      74     2548260 :   for (unsigned int i = 0; i < _order; ++i)
      75     1698840 :     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      849420 :   if (_order == 3 && _dt != _dt_old)
      85             :   {
      86      195294 :     Real x0 = 0.0;
      87      195294 :     Real x1 = _dt_old;
      88      195294 :     Real x2 = _dt + _dt_old;
      89      195294 :     Real y0 = (*_coupled_vars[2])[_qp];
      90      195294 :     Real y1 = (*_coupled_vars[1])[_qp];
      91      195294 :     Real y2 = (*_coupled_vars[0])[_qp];
      92      195294 :     Real term1 = (x2 - x0) * (y0 + (x2 - x0) * (y1 - y0) / (2.0 * (x1 - x0)));
      93      195294 :     Real term2 = (2.0 * x2 * x2 - x0 * x2 - x0 * x0 + 3.0 * x0 * x1 - 3.0 * x1 * x2) / 6.0;
      94      195294 :     Real term3 = (y2 - y1) / (x2 - x1) - (y1 - y0) / (x1 - x0);
      95      195294 :     integral_value = term1 + term2 * term3;
      96             :   }
      97             : 
      98      849420 :   return integral_value;
      99             : }