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 "FVPorousFlowHeatAdvection.h"
      11             : #include "PorousFlowDictator.h"
      12             : 
      13             : registerADMooseObject("PorousFlowApp", FVPorousFlowHeatAdvection);
      14             : 
      15             : InputParameters
      16          19 : FVPorousFlowHeatAdvection::validParams()
      17             : {
      18          19 :   InputParameters params = FVFluxKernel::validParams();
      19             :   RealVectorValue g(0, 0, -9.81);
      20          38 :   params.addParam<RealVectorValue>("gravity", g, "Gravity vector. Defaults to (0, 0, -9.81)");
      21          38 :   params.addRequiredParam<UserObjectName>("PorousFlowDictator",
      22             :                                           "The PorousFlowDictator UserObject");
      23          19 :   params.set<unsigned short>("ghost_layers") = 2;
      24          19 :   params.addClassDescription("Heat flux advected by the fluid");
      25          19 :   return params;
      26           0 : }
      27             : 
      28          10 : FVPorousFlowHeatAdvection::FVPorousFlowHeatAdvection(const InputParameters & params)
      29             :   : FVFluxKernel(params),
      30          10 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      31          10 :     _num_phases(_dictator.numPhases()),
      32          20 :     _density(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      33          10 :     _density_neighbor(
      34          10 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      35          20 :     _viscosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      36          10 :     _viscosity_neighbor(
      37          10 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      38          20 :     _enthalpy(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_qp")),
      39          10 :     _enthalpy_neighbor(
      40          10 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_qp")),
      41          20 :     _relperm(getADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      42          10 :     _relperm_neighbor(
      43          10 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      44          20 :     _permeability(getADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      45          10 :     _permeability_neighbor(
      46          10 :         getNeighborADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      47          20 :     _pressure(getADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      48          10 :     _pressure_neighbor(
      49          10 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      50          20 :     _grad_p(getADMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      51          30 :     _gravity(getParam<RealVectorValue>("gravity"))
      52             : {
      53          10 : }
      54             : 
      55             : ADReal
      56       59313 : FVPorousFlowHeatAdvection::computeQpResidual()
      57             : {
      58       59313 :   ADReal flux = 0.0;
      59             :   ADRealGradient pressure_grad;
      60             :   ADRealTensorValue mobility;
      61             : 
      62      118626 :   for (const auto p : make_range(_num_phases))
      63             :   {
      64             :     // If we are on a boundary face, use the gradient computed in _grad_p
      65       59313 :     if (onBoundary(*_face_info))
      66             :     {
      67        2326 :       const auto & gradp = -_grad_p[_qp][p];
      68        2326 :       pressure_grad = gradp + _density[_qp][p] * _gravity;
      69             : 
      70        4652 :       mobility = _enthalpy[_qp][p] * _relperm[_qp][p] * _permeability[_qp] * _density[_qp][p] /
      71        4652 :                  _viscosity[_qp][p];
      72             :     }
      73             :     else
      74             :     {
      75             :       // If we are on an internal face, calculate the gradient explicitly
      76       56987 :       const auto & p_elem = _pressure[_qp][p];
      77       56987 :       const auto & p_neighbor = _pressure_neighbor[_qp][p];
      78             : 
      79      113974 :       const auto gradp = (p_elem - p_neighbor) * _face_info->eCN() / _face_info->dCNMag();
      80             : 
      81      113974 :       const auto mobility_element = _enthalpy[_qp][p] * _relperm[_qp][p] * _permeability[_qp] *
      82      113974 :                                     _density[_qp][p] / _viscosity[_qp][p];
      83             : 
      84           0 :       const auto mobility_neighbor = _enthalpy_neighbor[_qp][p] * _relperm_neighbor[_qp][p] *
      85      113974 :                                      _permeability_neighbor[_qp] * _density_neighbor[_qp][p] /
      86      113974 :                                      _viscosity_neighbor[_qp][p];
      87             : 
      88       56987 :       pressure_grad = gradp + _density[_qp][p] * _gravity;
      89             : 
      90       56987 :       interpolate(Moose::FV::InterpMethod::Upwind,
      91             :                   mobility,
      92             :                   mobility_element,
      93             :                   mobility_neighbor,
      94             :                   pressure_grad,
      95       56987 :                   *_face_info,
      96             :                   true);
      97             :     }
      98             : 
      99       59313 :     flux += mobility * pressure_grad * _normal;
     100             :   }
     101             : 
     102       59313 :   return flux;
     103             : }