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         123 : FVPorousFlowAdvectiveFluxBC::validParams()
      18             : {
      19         123 :   InputParameters params = FVFluxBC::validParams();
      20             :   RealVectorValue g(0, 0, -9.81);
      21         246 :   params.addParam<RealVectorValue>("gravity", g, "Gravity vector. Defaults to (0, 0, -9.81)");
      22         246 :   params.addRequiredParam<UserObjectName>("PorousFlowDictator",
      23             :                                           "The PorousFlowDictator UserObject");
      24         246 :   params.addParam<unsigned int>("phase", 0, "The fluid phase for this BC");
      25         246 :   params.addParam<unsigned int>("fluid_component", 0, "The fluid component for this BC");
      26         123 :   params.addClassDescription("Advective Darcy flux boundary condition");
      27         246 :   params.addRequiredParam<Real>("porepressure_value", "The porepressure value on the boundary");
      28         123 :   return params;
      29           0 : }
      30             : 
      31          66 : FVPorousFlowAdvectiveFluxBC::FVPorousFlowAdvectiveFluxBC(const InputParameters & params)
      32             :   : FVFluxBC(params),
      33          66 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      34          66 :     _num_phases(_dictator.numPhases()),
      35         132 :     _phase(getParam<unsigned int>("phase")),
      36         132 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      37         132 :     _density(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      38          66 :     _density_neighbor(
      39          66 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      40         132 :     _viscosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      41          66 :     _viscosity_neighbor(
      42          66 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      43         132 :     _relperm(getADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      44          66 :     _relperm_neighbor(
      45          66 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      46          66 :     _mass_fractions(
      47          66 :         getADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      48          66 :     _mass_fractions_neighbor(
      49          66 :         getNeighborADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      50         132 :     _permeability(getADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      51          66 :     _permeability_neighbor(
      52          66 :         getNeighborADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      53         132 :     _pressure(getADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      54          66 :     _pressure_neighbor(
      55          66 :         getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      56         132 :     _gravity(getParam<RealVectorValue>("gravity")),
      57         198 :     _pp_value(getParam<Real>("porepressure_value"))
      58             : {
      59          66 :   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          66 :   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          66 :   if (_tid == 0)
      79             :   {
      80          57 :     auto diri_params = FVDirichletBC::validParams();
      81          57 :     diri_params.applySpecificParameters(_pars, {"variable", "boundary"});
      82          57 :     diri_params.addPrivateParam(MooseBase::app_param, &_app);
      83          57 :     diri_params.set<Real>("value") = _pp_value;
      84         114 :     _fv_problem.addFVBC("FVDirichletBC", name() + "_diri", diri_params);
      85          57 :     _fv_problem.fvBCsIntegrityCheck(false);
      86          57 :   }
      87          66 : }
      88             : 
      89             : ADReal
      90        2310 : FVPorousFlowAdvectiveFluxBC::computeQpResidual()
      91             : {
      92        2310 :   const bool out_of_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM);
      93             : 
      94        2310 :   const auto p_interior = out_of_elem ? _pressure[_qp][_phase] : _pressure_neighbor[_qp][_phase];
      95        2310 :   const auto delta_p = out_of_elem ? p_interior - _pp_value : _pp_value - p_interior;
      96        2310 :   const auto gradp = delta_p * _face_info->eCN() / _face_info->dCNMag();
      97             : 
      98             :   const auto mobility =
      99        4620 :       out_of_elem ? _mass_fractions[_qp][_phase][_fluid_component] * _relperm[_qp][_phase] *
     100        6930 :                         _permeability[_qp] * _density[_qp][_phase] / _viscosity[_qp][_phase]
     101        2310 :                   : _mass_fractions_neighbor[_qp][_phase][_fluid_component] *
     102        2310 :                         _relperm_neighbor[_qp][_phase] * _permeability_neighbor[_qp] *
     103        2310 :                         _density_neighbor[_qp][_phase] / _viscosity_neighbor[_qp][_phase];
     104             : 
     105        2310 :   const auto pressure_grad = gradp + _density[_qp][_phase] * _gravity;
     106             : 
     107        2310 :   return mobility * pressure_grad * _normal;
     108             : }