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       14603 : LinearFVAdvectionDiffusionFunctorRobinBC::validParams()
      16             : {
      17       14603 :   InputParameters params = LinearFVAdvectionDiffusionBC::validParams();
      18       29206 :   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       43809 :   params.addParam<MooseFunctorName>(
      25       29206 :       "alpha", 1.0, "Functor for the coefficient of the normal gradient term.");
      26       58412 :   params.addParam<MooseFunctorName>("beta", 1.0, "Functor for the coefficient of the scalar term.");
      27       29206 :   params.addParam<MooseFunctorName>(
      28       29206 :       "gamma", 1.0, "Functor for the constant term on the RHS of the Robin BC.");
      29       14603 :   return params;
      30           0 : }
      31             : 
      32         168 : LinearFVAdvectionDiffusionFunctorRobinBC::LinearFVAdvectionDiffusionFunctorRobinBC(
      33         168 :     const InputParameters & parameters)
      34             :   : LinearFVAdvectionDiffusionBC(parameters),
      35         168 :     _alpha(getFunctor<Real>("alpha")),
      36         336 :     _beta(getFunctor<Real>("beta")),
      37         504 :     _gamma(getFunctor<Real>("gamma"))
      38             : {
      39         168 :   _var.computeCellGradients();
      40             : 
      41         168 :   if (_alpha.isConstant())
      42             :   {
      43             :     // We check if we can parse the value to a number and if yes, we throw an error if it is 0
      44         336 :     std::istringstream ss(getParam<MooseFunctorName>("alpha"));
      45             :     Real real_value;
      46         168 :     if (ss >> real_value && ss.eof())
      47         168 :       if (MooseUtils::isZero(real_value))
      48           0 :         paramError("alpha",
      49             :                    "This value shall not be 0. Use a Dirichlet boundary condition instead!");
      50         168 :   }
      51         168 : }
      52             : 
      53             : Real
      54        9968 : LinearFVAdvectionDiffusionFunctorRobinBC::getAlpha(Moose::FaceArg face, Moose::StateArg state) const
      55             : {
      56        9968 :   return _alpha(face, state);
      57             : }
      58             : 
      59             : Real
      60        9968 : LinearFVAdvectionDiffusionFunctorRobinBC::getBeta(Moose::FaceArg face, Moose::StateArg state) const
      61             : {
      62        9968 :   return _beta(face, state);
      63             : }
      64             : 
      65             : Real
      66        7476 : LinearFVAdvectionDiffusionFunctorRobinBC::getGamma(Moose::FaceArg face, Moose::StateArg state) const
      67             : {
      68        7476 :   return _gamma(face, state);
      69             : }
      70             : 
      71             : Real
      72        4984 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValue() const
      73             : {
      74        4984 :   const auto face = singleSidedFaceArg(_current_face_info);
      75             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
      76             :               "This should not be assigned on an internal face!");
      77        4984 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
      78        4984 :                                ? _current_face_info->elemInfo()
      79           0 :                                : _current_face_info->neighborInfo();
      80        4984 :   const auto state = determineState();
      81             : 
      82        4984 :   const auto alpha = getAlpha(face, state);
      83        4984 :   const auto beta = getBeta(face, state);
      84        4984 :   const auto gamma = getGamma(face, state);
      85             : 
      86        4984 :   const auto phi = _var.getElemValue(*elem_info, state);
      87        4984 :   const auto grad_phi = _var.gradSln(*elem_info);
      88             : 
      89        4984 :   const auto & nhat = _current_face_info->normal();
      90             : 
      91        4984 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
      92        4984 :   const auto projection = d_cf * nhat;
      93        4984 :   const auto vc = d_cf - (projection * nhat);
      94        4984 :   return ((alpha * phi) + (alpha * grad_phi * vc) + (gamma * projection)) /
      95        4984 :          (alpha + (beta * projection));
      96             : }
      97             : 
      98             : Real
      99           0 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryNormalGradient() const
     100             : {
     101           0 :   const auto face = singleSidedFaceArg(_current_face_info);
     102           0 :   const auto state = determineState();
     103           0 :   const auto alpha = getAlpha(face, state);
     104             :   mooseAssert(!MooseUtils::isZero(alpha), "Alpha should not be 0!");
     105           0 :   const auto beta = _beta(face, state);
     106           0 :   const auto gamma = _gamma(face, state);
     107           0 :   return (gamma - beta * computeBoundaryValue()) / alpha;
     108             : }
     109             : 
     110             : // implicit terms for advection kernel
     111             : Real
     112         448 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValueMatrixContribution() const
     113             : {
     114         448 :   const auto face = singleSidedFaceArg(_current_face_info);
     115         448 :   const auto state = determineState();
     116         448 :   const auto alpha = getAlpha(face, state);
     117         448 :   const auto beta = getBeta(face, state);
     118         448 :   const auto & nhat = _current_face_info->normal();
     119             : 
     120         448 :   return alpha / (alpha + (beta * computeCellToFaceVector() * nhat));
     121             : }
     122             : 
     123             : // explicit terms for advection kernel
     124             : Real
     125         448 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryValueRHSContribution() const
     126             : {
     127         448 :   const auto face = singleSidedFaceArg(_current_face_info);
     128         448 :   const auto state = determineState();
     129             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
     130             :               "This should not be assigned on an internal face!");
     131         448 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     132         448 :                                ? _current_face_info->elemInfo()
     133           0 :                                : _current_face_info->neighborInfo();
     134             : 
     135         448 :   const auto alpha = getAlpha(face, state);
     136         448 :   const auto beta = getBeta(face, state);
     137         448 :   const auto gamma = getGamma(face, state);
     138             : 
     139         448 :   const auto & grad_phi = _var.gradSln(*elem_info);
     140             : 
     141         448 :   const auto & nhat = _current_face_info->normal();
     142             : 
     143         448 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
     144         448 :   const auto projection = d_cf * nhat;
     145         448 :   const auto vc = d_cf - (projection * nhat); // correction vector for non-orthogonal cells
     146             : 
     147         448 :   return (gamma * projection / (alpha + (beta * projection))) +
     148         448 :          (alpha * grad_phi * vc / (alpha + (beta * projection)));
     149             : }
     150             : 
     151             : // implicit terms for diffusion kernel
     152             : Real
     153        2044 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryGradientMatrixContribution() const
     154             : {
     155        2044 :   const auto face = singleSidedFaceArg(_current_face_info);
     156        2044 :   const auto state = determineState();
     157             : 
     158        2044 :   const auto alpha = getAlpha(face, state);
     159        2044 :   const auto beta = getBeta(face, state);
     160             : 
     161        2044 :   const auto & nhat = _current_face_info->normal();
     162             : 
     163        2044 :   return beta / (alpha + (beta * computeCellToFaceVector() * nhat));
     164             : }
     165             : 
     166             : // explicit terms for diffusion kernel
     167             : Real
     168        2044 : LinearFVAdvectionDiffusionFunctorRobinBC::computeBoundaryGradientRHSContribution() const
     169             : {
     170             :   mooseAssert(_current_face_type != FaceInfo::VarFaceNeighbors::BOTH,
     171             :               "This should not be assigned on an internal face!");
     172        2044 :   const auto & elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     173        2044 :                                ? _current_face_info->elemInfo()
     174           0 :                                : _current_face_info->neighborInfo();
     175        2044 :   const auto face = singleSidedFaceArg(_current_face_info);
     176        2044 :   const auto state = determineState();
     177        2044 :   const auto & grad_phi = _var.gradSln(*elem_info);
     178             : 
     179        2044 :   const auto alpha = getAlpha(face, state);
     180        2044 :   const auto beta = getBeta(face, state);
     181        2044 :   const auto gamma = getGamma(face, state);
     182             : 
     183        2044 :   const auto & nhat = _current_face_info->normal();
     184             : 
     185        2044 :   const auto d_cf = computeCellToFaceVector(); // vector from boundary cell centre to boundary face
     186        2044 :   const auto projection = d_cf * nhat;
     187        2044 :   const auto vc = d_cf - (projection * nhat); // correction vector for non-orthogonal cells
     188             : 
     189        2044 :   return (gamma / alpha) + (-beta * gamma * projection / alpha / (alpha + (beta * projection))) +
     190        2044 :          (-beta * grad_phi * vc / (alpha + (beta * projection)));
     191             : }