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 "FVPorousFlowAdvectiveFluxBC.h"
      11             : #include "PorousFlowDictator.h"
      12             : #include "FVDirichletBC.h"
      13             : 
      14             : registerADMooseObject("PorousFlowApp", FVPorousFlowAdvectiveFluxBC);
      15             : 
      16             : InputParameters
      17          57 : FVPorousFlowAdvectiveFluxBC::validParams()
      18             : {
      19          57 :   InputParameters params = FVFluxBC::validParams();
      20             :   RealVectorValue g(0, 0, -9.81);
      21         114 :   params.addParam<RealVectorValue>("gravity", g, "Gravity vector. Defaults to (0, 0, -9.81)");
      22         114 :   params.addRequiredParam<UserObjectName>("PorousFlowDictator",
      23             :                                           "The PorousFlowDictator UserObject");
      24         114 :   params.addParam<unsigned int>("phase", 0, "The fluid phase for this BC");
      25         114 :   params.addParam<unsigned int>("fluid_component", 0, "The fluid component for this BC");
      26          57 :   params.addClassDescription("Advective Darcy flux boundary condition");
      27         114 :   params.addRequiredParam<Real>("porepressure_value", "The porepressure value on the boundary");
      28          57 :   return params;
      29           0 : }
      30             : 
      31          30 : FVPorousFlowAdvectiveFluxBC::FVPorousFlowAdvectiveFluxBC(const InputParameters & params)
      32             :   : FVFluxBC(params),
      33          30 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      34          30 :     _num_phases(_dictator.numPhases()),
      35          60 :     _phase(getParam<unsigned int>("phase")),
      36          60 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      37          60 :     _density(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      38          30 :     _density_neighbor(
      39          30 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      40          60 :     _viscosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      41          30 :     _viscosity_neighbor(
      42          30 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      43          60 :     _relperm(getADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      44          30 :     _relperm_neighbor(
      45          30 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      46          30 :     _mass_fractions(
      47          30 :         getADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      48          30 :     _mass_fractions_neighbor(
      49          30 :         getNeighborADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      50          60 :     _permeability(getADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      51          30 :     _permeability_neighbor(
      52          30 :         getNeighborADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      53          60 :     _pressure(getADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      54          30 :     _pressure_neighbor(
      55          30 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      56          60 :     _gravity(getParam<RealVectorValue>("gravity")),
      57          90 :     _pp_value(getParam<Real>("porepressure_value"))
      58             : {
      59          30 :   if (_phase >= _num_phases)
      60           0 :     paramError(
      61             :         "phase",
      62             :         "The Dictator proclaims that the maximum fluid phase index in this simulation is ",
      63             :         _num_phases - 1,
      64             :         " whereas you have used ",
      65             :         _phase,
      66             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
      67             : 
      68          30 :   if (_fluid_component >= _dictator.numComponents())
      69           0 :     paramError(
      70             :         "fluid_component",
      71             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
      72           0 :         _dictator.numComponents() - 1,
      73             :         " whereas you have used ",
      74           0 :         _fluid_component,
      75             :         ". Remember that indexing starts at 0.");
      76             : 
      77             :   // Add a FVDirichletBC to boundary for cell gradient computation
      78          30 :   if (_tid == 0)
      79             :   {
      80          27 :     auto diri_params = FVDirichletBC::validParams();
      81          27 :     diri_params.applySpecificParameters(_pars, {"variable", "boundary"});
      82          27 :     diri_params.addPrivateParam(MooseBase::app_param, &_app);
      83          27 :     diri_params.set<Real>("value") = _pp_value;
      84          54 :     _fv_problem.addFVBC("FVDirichletBC", name() + "_diri", diri_params);
      85          27 :     _fv_problem.fvBCsIntegrityCheck(false);
      86          27 :   }
      87          30 : }
      88             : 
      89             : ADReal
      90        1590 : FVPorousFlowAdvectiveFluxBC::computeQpResidual()
      91             : {
      92        1590 :   const bool out_of_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM);
      93             : 
      94        1590 :   const auto p_interior = out_of_elem ? _pressure[_qp][_phase] : _pressure_neighbor[_qp][_phase];
      95        1590 :   const auto delta_p = out_of_elem ? p_interior - _pp_value : _pp_value - p_interior;
      96        1590 :   const auto gradp = delta_p * _face_info->eCN() / _face_info->dCNMag();
      97             : 
      98             :   const auto mobility =
      99        3180 :       out_of_elem ? _mass_fractions[_qp][_phase][_fluid_component] * _relperm[_qp][_phase] *
     100        4770 :                         _permeability[_qp] * _density[_qp][_phase] / _viscosity[_qp][_phase]
     101        1590 :                   : _mass_fractions_neighbor[_qp][_phase][_fluid_component] *
     102        1590 :                         _relperm_neighbor[_qp][_phase] * _permeability_neighbor[_qp] *
     103        1590 :                         _density_neighbor[_qp][_phase] / _viscosity_neighbor[_qp][_phase];
     104             : 
     105        1590 :   const auto pressure_grad = gradp + _density[_qp][_phase] * _gravity;
     106             : 
     107        1590 :   return mobility * pressure_grad * _normal;
     108             : }