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         104 : PorousFlowEnthalpySink::validParams()
      25             : {
      26         104 :   InputParameters params = IntegratedBC::validParams();
      27         104 :   params += PorousFlowSinkBC::validParamsCommon();
      28         208 :   params.addRequiredParam<UserObjectName>(
      29             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      30         104 :   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         104 :   return params;
      38           0 : }
      39             : 
      40          58 : PorousFlowEnthalpySink::PorousFlowEnthalpySink(const InputParameters & parameters)
      41             :   : IntegratedBC(parameters),
      42          58 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      43          77 :     _pressure(isCoupled("porepressure_var") ? &coupledValue("porepressure_var") : nullptr),
      44         194 :     _ph(isParamValid("fluid_phase") ? getParam<unsigned int>("fluid_phase")
      45             :                                     : libMesh::invalid_uint),
      46          58 :     _m_func(getFunction("flux_function")),
      47         116 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
      48          58 :     _dpp_dvar(
      49          58 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_nodal_dvar")),
      50         116 :     _T_in(getParam<Real>("T_in")),
      51         116 :     _fp(getUserObject<SinglePhaseFluidProperties>("fp"))
      52             : {
      53          96 :   if ((_pressure != nullptr) && (isParamValid("fluid_phase")))
      54           2 :     mooseError(name(), ": Cannot specify both pressure and pore pressure.");
      55             : 
      56         171 :   if ((_pressure == nullptr) && (!isParamValid("fluid_phase")))
      57           2 :     mooseError(name(), ": You have to specify either 'pressure' or 'fluid_phase'.");
      58             : 
      59         108 :   if (isParamValid("fluid_phase"))
      60          37 :     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          52 : }
      67             : 
      68             : Real
      69       33072 : PorousFlowEnthalpySink::computeQpResidual()
      70             : {
      71             :   Real h;
      72       33072 :   if (_pressure)
      73        7040 :     h = _fp.h_from_p_T((*_pressure)[_qp], _T_in);
      74             :   else
      75       26032 :     h = _fp.h_from_p_T(_pp[_i][_ph], _T_in);
      76             : 
      77       33072 :   return _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]) * h;
      78             : }
      79             : 
      80             : Real
      81      193864 : PorousFlowEnthalpySink::computeQpJacobian()
      82             : {
      83      193864 :   return jac(_var.number());
      84             : }
      85             : 
      86             : Real
      87      193864 : PorousFlowEnthalpySink::computeQpOffDiagJacobian(unsigned int jvar)
      88             : {
      89      193864 :   return jac(jvar);
      90             : }
      91             : 
      92             : Real
      93      387728 : PorousFlowEnthalpySink::jac(unsigned int jvar) const
      94             : {
      95      387728 :   if (_dictator.notPorousFlowVariable(jvar))
      96             :     return 0.0;
      97             : 
      98      387728 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
      99             : 
     100      387728 :   if (_i != _j)
     101             :     return 0.0;
     102             : 
     103             :   Real jac;
     104       48584 :   if (_pressure)
     105             :   {
     106             :     jac = 0.;
     107             :   }
     108             :   else
     109             :   {
     110             :     Real h, dh_dpp, dh_dT;
     111       38536 :     _fp.h_from_p_T(_pp[_i][_ph], _T_in, h, dh_dpp, dh_dT);
     112       38536 :     jac = dh_dpp * _dpp_dvar[_i][_ph][pvar];
     113             :   }
     114       48584 :   return _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]) * jac;
     115             : }