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 "PorousFlowEnthalpySink.h"
      11             : #include "Function.h"
      12             : #include "PorousFlowDictator.h"
      13             : #include "SinglePhaseFluidProperties.h"
      14             : #include "MooseVariable.h"
      15             : #include "PorousFlowSinkBC.h"
      16             : 
      17             : #include "libmesh/quadrature.h"
      18             : 
      19             : #include <iostream>
      20             : 
      21             : registerMooseObject("PorousFlowApp", PorousFlowEnthalpySink);
      22             : 
      23             : InputParameters
      24         190 : PorousFlowEnthalpySink::validParams()
      25             : {
      26         190 :   InputParameters params = IntegratedBC::validParams();
      27         190 :   params += PorousFlowSinkBC::validParamsCommon();
      28         380 :   params.addRequiredParam<UserObjectName>(
      29             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      30         190 :   params.addClassDescription(
      31             :       "Applies a source equal to the product of the mass flux and the "
      32             :       "fluid enthalpy. The enthalpy is computed at temperature T_in and pressure equal to the "
      33             :       "porepressure in the porous medium, if fluid_phase is given, otherwise at the supplied "
      34             :       "porepressure. Hence this adds heat energy to the porous medium at rate corresponding to a "
      35             :       "fluid being injected at (porepressure, T_in) at rate (-flux_function).");
      36             : 
      37         190 :   return params;
      38           0 : }
      39             : 
      40         106 : PorousFlowEnthalpySink::PorousFlowEnthalpySink(const InputParameters & parameters)
      41             :   : IntegratedBC(parameters),
      42         106 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      43         141 :     _pressure(isCoupled("porepressure_var") ? &coupledValue("porepressure_var") : nullptr),
      44         354 :     _ph(isParamValid("fluid_phase") ? getParam<unsigned int>("fluid_phase")
      45             :                                     : libMesh::invalid_uint),
      46         106 :     _m_func(getFunction("flux_function")),
      47         212 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
      48         106 :     _dpp_dvar(
      49         106 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_nodal_dvar")),
      50         212 :     _T_in(getParam<Real>("T_in")),
      51         212 :     _fp(getUserObject<SinglePhaseFluidProperties>("fp"))
      52             : {
      53         176 :   if ((_pressure != nullptr) && (isParamValid("fluid_phase")))
      54           2 :     mooseError(name(), ": Cannot specify both pressure and pore pressure.");
      55             : 
      56         315 :   if ((_pressure == nullptr) && (!isParamValid("fluid_phase")))
      57           2 :     mooseError(name(), ": You have to specify either 'pressure' or 'fluid_phase'.");
      58             : 
      59         204 :   if (isParamValid("fluid_phase"))
      60          69 :     if (_ph >= _dictator.numPhases())
      61           2 :       mooseError(name(),
      62             :                  ": Specified 'fluid_phase' is larger than the number of phases available in the "
      63             :                  "simulation (",
      64           2 :                  _dictator.numPhases(),
      65             :                  ").");
      66         100 : }
      67             : 
      68             : Real
      69       48632 : PorousFlowEnthalpySink::computeQpResidual()
      70             : {
      71             :   Real h;
      72       48632 :   if (_pressure)
      73       10208 :     h = _fp.h_from_p_T((*_pressure)[_qp], _T_in);
      74             :   else
      75       38424 :     h = _fp.h_from_p_T(_pp[_i][_ph], _T_in);
      76             : 
      77       48632 :   return _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]) * h;
      78             : }
      79             : 
      80             : Real
      81      286100 : PorousFlowEnthalpySink::computeQpJacobian()
      82             : {
      83      286100 :   return jac(_var.number());
      84             : }
      85             : 
      86             : Real
      87      286100 : PorousFlowEnthalpySink::computeQpOffDiagJacobian(unsigned int jvar)
      88             : {
      89      286100 :   return jac(jvar);
      90             : }
      91             : 
      92             : Real
      93      572200 : PorousFlowEnthalpySink::jac(unsigned int jvar) const
      94             : {
      95      572200 :   if (_dictator.notPorousFlowVariable(jvar))
      96             :     return 0.0;
      97             : 
      98      572200 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
      99             : 
     100      572200 :   if (_i != _j)
     101             :     return 0.0;
     102             : 
     103             :   Real jac;
     104       71668 :   if (_pressure)
     105             :   {
     106             :     jac = 0.;
     107             :   }
     108             :   else
     109             :   {
     110             :     Real h, dh_dpp, dh_dT;
     111       56996 :     _fp.h_from_p_T(_pp[_i][_ph], _T_in, h, dh_dpp, dh_dT);
     112       56996 :     jac = dh_dpp * _dpp_dvar[_i][_ph][pvar];
     113             :   }
     114       71668 :   return _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]) * jac;
     115             : }