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 "FVPorousFlowAdvectiveFlux.h"
      11             : #include "PorousFlowDictator.h"
      12             : 
      13             : registerADMooseObject("PorousFlowApp", FVPorousFlowAdvectiveFlux);
      14             : 
      15             : InputParameters
      16         821 : FVPorousFlowAdvectiveFlux::validParams()
      17             : {
      18         821 :   InputParameters params = FVFluxKernel::validParams();
      19             :   RealVectorValue g(0, 0, -9.81);
      20        1642 :   params.addParam<RealVectorValue>("gravity", g, "Gravity vector. Defaults to (0, 0, -9.81)");
      21        1642 :   params.addRequiredParam<UserObjectName>("PorousFlowDictator",
      22             :                                           "The PorousFlowDictator UserObject");
      23        1642 :   params.addParam<unsigned int>("fluid_component", 0, "The fluid component for this kernel");
      24         821 :   params.set<unsigned short>("ghost_layers") = 2;
      25         821 :   params.addClassDescription("Advective Darcy flux");
      26         821 :   return params;
      27           0 : }
      28             : 
      29         441 : FVPorousFlowAdvectiveFlux::FVPorousFlowAdvectiveFlux(const InputParameters & params)
      30             :   : FVFluxKernel(params),
      31         441 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      32         441 :     _num_phases(_dictator.numPhases()),
      33         882 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      34         882 :     _density(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      35         441 :     _density_neighbor(
      36         441 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      37         882 :     _viscosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      38         441 :     _viscosity_neighbor(
      39         441 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      40         882 :     _relperm(getADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      41         441 :     _relperm_neighbor(
      42         441 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      43         441 :     _mass_fractions(
      44         441 :         getADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      45         441 :     _mass_fractions_neighbor(
      46         441 :         getNeighborADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      47         882 :     _permeability(getADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      48         441 :     _permeability_neighbor(
      49         441 :         getNeighborADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      50         882 :     _pressure(getADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      51         441 :     _pressure_neighbor(
      52         441 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      53         882 :     _grad_p(getADMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      54        1323 :     _gravity(getParam<RealVectorValue>("gravity"))
      55             : {
      56         441 :   if (_fluid_component >= _dictator.numComponents())
      57           0 :     paramError(
      58             :         "fluid_component",
      59             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
      60           0 :         _dictator.numComponents() - 1,
      61             :         " whereas you have used ",
      62           0 :         _fluid_component,
      63             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
      64         441 : }
      65             : 
      66             : ADReal
      67      308117 : FVPorousFlowAdvectiveFlux::computeQpResidual()
      68             : {
      69      308117 :   ADReal flux = 0.0;
      70             :   ADRealGradient pressure_grad;
      71             :   ADRealTensorValue mobility;
      72             : 
      73      803916 :   for (const auto p : make_range(_num_phases))
      74             :   {
      75             :     // If we are on a boundary face, use the gradient computed in _grad_p
      76      495799 :     if (onBoundary(*_face_info))
      77             :     {
      78       10563 :       const auto & gradp = -_grad_p[_qp][p];
      79       10563 :       pressure_grad = gradp + _density[_qp][p] * _gravity;
      80             : 
      81       21126 :       mobility = _mass_fractions[_qp][p][_fluid_component] * _relperm[_qp][p] * _permeability[_qp] *
      82       31689 :                  _density[_qp][p] / _viscosity[_qp][p];
      83             :     }
      84             :     else
      85             :     {
      86             :       // If we are on an internal face, calculate the gradient explicitly
      87      485236 :       const auto & p_elem = _pressure[_qp][p];
      88      485236 :       const auto & p_neighbor = _pressure_neighbor[_qp][p];
      89             : 
      90      970472 :       const auto gradp = (p_elem - p_neighbor) * _face_info->eCN() / _face_info->dCNMag();
      91             : 
      92      485236 :       const auto mobility_element = _mass_fractions[_qp][p][_fluid_component] * _relperm[_qp][p] *
      93      970472 :                                     _permeability[_qp] * _density[_qp][p] / _viscosity[_qp][p];
      94             : 
      95      970472 :       const auto mobility_neighbor = _mass_fractions_neighbor[_qp][p][_fluid_component] *
      96      970472 :                                      _relperm_neighbor[_qp][p] * _permeability_neighbor[_qp] *
      97      970472 :                                      _density_neighbor[_qp][p] / _viscosity_neighbor[_qp][p];
      98             : 
      99      485236 :       pressure_grad = gradp + _density[_qp][p] * _gravity;
     100             : 
     101      485236 :       interpolate(Moose::FV::InterpMethod::Upwind,
     102             :                   mobility,
     103             :                   mobility_element,
     104             :                   mobility_neighbor,
     105             :                   pressure_grad,
     106      485236 :                   *_face_info,
     107             :                   true);
     108             :     }
     109             : 
     110      495799 :     flux += mobility * pressure_grad * _normal;
     111             :   }
     112             : 
     113      308117 :   return flux;
     114             : }