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 "LinearFVAdvectionDiffusionFunctorRobinBC.h"
      11             : 
      12             : registerMooseObject("MooseApp", LinearFVAdvectionDiffusionFunctorRobinBC);
      13             : 
      14             : InputParameters
      15       14601 : LinearFVAdvectionDiffusionFunctorRobinBC::validParams()
      16             : {
      17       14601 :   InputParameters params = LinearFVAdvectionDiffusionBC::validParams();
      18       14601 :   params.addClassDescription(
      19             :       "Adds a Robin BC of the form \\alpha * \\nabla \\phi*n + \\beta * \\phi = \\gamma, "
      20             :       "which can be used for the assembly of linear "
      21             :       "finite volume system and whose face values are determined using "
      22             :       "three functors. This kernel is "
      23             :       "only designed to work with advection-diffusion problems.");
      24       14601 :   params.addRequiredParam<MooseFunctorName>(
      25             :       "alpha", "The functor which is the coefficient of the normal gradient term.");
      26       14601 :   params.addRequiredParam<MooseFunctorName>(
      27             :       "beta", "The functor which is the coefficient of the scalar term.");
      28       14601 :   params.addRequiredParam<MooseFunctorName>(
      29             :       "gamma", "The functor which is the constant term on the RHS of the Robin BC expression.");
      30       14601 :   return params;
      31           0 : }
      32             : 
      33         168 : LinearFVAdvectionDiffusionFunctorRobinBC::LinearFVAdvectionDiffusionFunctorRobinBC(
      34         168 :     const InputParameters & parameters)
      35             :   : LinearFVAdvectionDiffusionBC(parameters),
      36         168 :     _alpha(getFunctor<Real>("alpha")),
      37         168 :     _beta(getFunctor<Real>("beta")),
      38         336 :     _gamma(getFunctor<Real>("gamma"))
      39             : {
      40         168 :   _var.computeCellGradients();
      41             : 
      42         168 :   if (_alpha.isConstant())
      43             :   {
      44             :     // We check if we can parse the value to a number and if yes, we throw an error if it is 0
      45         168 :     std::istringstream ss(getParam<MooseFunctorName>("alpha"));
      46             :     Real real_value;
      47         168 :     if (ss >> real_value && ss.eof())
      48         168 :       if (MooseUtils::isZero(real_value))
      49           0 :         paramError("alpha",
      50             :                    "This value shall not be 0. Use a Dirichlet boundary condition instead!");
      51         168 :   }
      52         168 : }
      53             : 
      54             : Real
      55        4984 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValue() const
      56             : {
      57        4984 :   const auto face = singleSidedFaceArg(_current_face_info);
      58             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
      59             :               "This should not be assigned on an internal face!");
      60        4984 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
      61        4984 :                                ? _current_face_info->elemInfo()
      62           0 :                                : _current_face_info->neighborInfo();
      63        4984 :   const auto state = determineState();
      64             : 
      65        4984 :   const auto alpha = _alpha(face, state);
      66        4984 :   const auto beta = _beta(face, state);
      67        4984 :   const auto gamma = _gamma(face, state);
      68             : 
      69        4984 :   const auto phi = _var.getElemValue(*elem_info, state);
      70        4984 :   const auto grad_phi = _var.gradSln(*elem_info);
      71             : 
      72        4984 :   const auto & nhat = _current_face_info->normal();
      73             : 
      74        4984 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
      75        4984 :   const auto projection = d_cf * nhat;
      76        4984 :   const auto vc = d_cf - (projection * nhat);
      77        4984 :   return ((alpha * phi) + (alpha * grad_phi * vc) + (gamma * projection)) /
      78        4984 :          (alpha + (beta * projection));
      79             : }
      80             : 
      81             : Real
      82           0 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryNormalGradient() const
      83             : {
      84           0 :   const auto face = singleSidedFaceArg(_current_face_info);
      85           0 :   const auto state = determineState();
      86           0 :   const auto alpha = _alpha(face, state);
      87             :   mooseAssert(!MooseUtils::isZero(alpha), "Alpha should not be 0!");
      88           0 :   const auto beta = _beta(face, state);
      89           0 :   const auto gamma = _gamma(face, state);
      90           0 :   return (gamma - beta * computeBoundaryValue()) / alpha;
      91             : }
      92             : 
      93             : // implicit terms for advection kernel
      94             : Real
      95         448 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValueMatrixContribution() const
      96             : {
      97         448 :   const auto face = singleSidedFaceArg(_current_face_info);
      98         448 :   const auto state = determineState();
      99         448 :   const auto alpha = _alpha(face, state);
     100         448 :   const auto beta = _beta(face, state);
     101         448 :   const auto & nhat = _current_face_info->normal();
     102             : 
     103         448 :   return alpha / (alpha + (beta * computeCellToFaceVector() * nhat));
     104             : }
     105             : 
     106             : // explicit terms for advection kernel
     107             : Real
     108         448 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValueRHSContribution() const
     109             : {
     110         448 :   const auto face = singleSidedFaceArg(_current_face_info);
     111         448 :   const auto state = determineState();
     112             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
     113             :               "This should not be assigned on an internal face!");
     114         448 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     115         448 :                                ? _current_face_info->elemInfo()
     116           0 :                                : _current_face_info->neighborInfo();
     117         448 :   const auto alpha = _alpha(face, state);
     118         448 :   const auto beta = _beta(face, state);
     119         448 :   const auto gamma = _gamma(face, state);
     120         448 :   const auto & grad_phi = _var.gradSln(*elem_info);
     121             : 
     122         448 :   const auto & nhat = _current_face_info->normal();
     123             : 
     124         448 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
     125         448 :   const auto projection = d_cf * nhat;
     126         448 :   const auto vc = d_cf - (projection * nhat); // correction vector for non-orthogonal cells
     127             : 
     128         448 :   return (gamma * projection / (alpha + (beta * projection))) +
     129         448 :          (alpha * grad_phi * vc / (alpha + (beta * projection)));
     130             : }
     131             : 
     132             : // implicit terms for diffusion kernel
     133             : Real
     134        2044 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryGradientMatrixContribution() const
     135             : {
     136        2044 :   const auto face = singleSidedFaceArg(_current_face_info);
     137        2044 :   const auto state = determineState();
     138             : 
     139        2044 :   const auto alpha = _alpha(face, state);
     140        2044 :   const auto beta = _beta(face, state);
     141        2044 :   const auto & nhat = _current_face_info->normal();
     142             : 
     143        2044 :   return beta / (alpha + (beta * computeCellToFaceVector() * nhat));
     144             : }
     145             : 
     146             : // explicit terms for diffusion kernel
     147             : Real
     148        2044 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryGradientRHSContribution() const
     149             : {
     150             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
     151             :               "This should not be assigned on an internal face!");
     152        2044 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     153        2044 :                                ? _current_face_info->elemInfo()
     154           0 :                                : _current_face_info->neighborInfo();
     155        2044 :   const auto face = singleSidedFaceArg(_current_face_info);
     156        2044 :   const auto state = determineState();
     157        2044 :   const auto & grad_phi = _var.gradSln(*elem_info);
     158             : 
     159        2044 :   const auto alpha = _alpha(face, state);
     160        2044 :   const auto beta = _beta(face, state);
     161        2044 :   const auto gamma = _gamma(face, state);
     162             : 
     163        2044 :   const auto & nhat = _current_face_info->normal();
     164             : 
     165        2044 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
     166        2044 :   const auto projection = d_cf * nhat;
     167        2044 :   const auto vc = d_cf - (projection * nhat); // correction vector for non-orthogonal cells
     168             : 
     169        2044 :   return (gamma / alpha) + (-beta * gamma * projection / alpha / (alpha + (beta * projection))) +
     170        2044 :          (-beta * grad_phi * vc / (alpha + (beta * projection)));
     171             : }