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 "LinearFVDiffusion.h"
      11             : #include "Assembly.h"
      12             : #include "SubProblem.h"
      13             : #include "LinearFVAdvectionDiffusionBC.h"
      14             : 
      15             : registerMooseObject("MooseApp", LinearFVDiffusion);
      16             : 
      17             : InputParameters
      18        4151 : LinearFVDiffusion::validParams()
      19             : {
      20        4151 :   InputParameters params = LinearFVFluxKernel::validParams();
      21        8302 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23       12453 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25        8302 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       16604 :   params.addParam<MooseFunctorName>("diffusion_coeff", 1.0, "The diffusion coefficient.");
      28       12453 :   params.addParam<InterpolationMethodName>(
      29             :       "coeff_interp_method",
      30             :       "Optional finite volume interpolation method used to compute a face-centered diffusion "
      31             :       "coefficient. If omitted, the functor is evaluated directly on the face.");
      32        4151 :   return params;
      33           0 : }
      34             : 
      35         545 : LinearFVDiffusion::LinearFVDiffusion(const InputParameters & params)
      36             :   : LinearFVFluxKernel(params),
      37             :     FVInterpolationMethodInterface(this),
      38         545 :     _diffusion_coeff(getFunctor<Real>("diffusion_coeff")),
      39        1090 :     _coeff_interp_method(isParamValid("coeff_interp_method")
      40         545 :                              ? &getFVFaceInterpolationMethod(
      41             :                                    getParam<InterpolationMethodName>("coeff_interp_method"))
      42             :                              : nullptr),
      43        1090 :     _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
      44         545 :     _flux_matrix_contribution(0.0),
      45         545 :     _flux_rhs_contribution(0.0),
      46         545 :     _cached_face_diffusivity(false),
      47         545 :     _face_diffusivity(0.0)
      48             : {
      49         545 :   if (_use_nonorthogonal_correction)
      50         348 :     _var.computeCellGradients();
      51         545 : }
      52             : 
      53             : Real
      54      946323 : LinearFVDiffusion::faceDiffusivity() const
      55             : {
      56      946323 :   if (!_cached_face_diffusivity)
      57             :   {
      58             :     mooseAssert(_current_face_type == FaceInfo::VarFaceNeighbors::BOTH,
      59             :                 "faceDiffusivity() is only valid for two-sided internal faces.");
      60             : 
      61      521659 :     const auto state = determineState();
      62             : 
      63      521659 :     if (_coeff_interp_method)
      64           0 :       _face_diffusivity =
      65           0 :           _coeff_interp_method->interpolate(_diffusion_coeff, *_current_face_info, state);
      66             :     else
      67      521659 :       _face_diffusivity = _diffusion_coeff(makeCDFace(*_current_face_info), state);
      68             : 
      69      521659 :     _cached_face_diffusivity = true;
      70             :   }
      71             : 
      72      946323 :   return _face_diffusivity;
      73             : }
      74             : 
      75             : void
      76      578015 : LinearFVDiffusion::setupFaceData(const FaceInfo * face_info)
      77             : {
      78      578015 :   LinearFVFluxKernel::setupFaceData(face_info);
      79      578015 :   _cached_face_diffusivity = false;
      80      578015 : }
      81             : 
      82             : void
      83         565 : LinearFVDiffusion::initialSetup()
      84             : {
      85        2083 :   for (const auto bc : _var.getBoundaryConditionMap())
      86        1518 :     if (!dynamic_cast<const LinearFVAdvectionDiffusionBC *>(bc.second))
      87           0 :       mooseError(
      88           0 :           bc.second->type(), " is not a compatible boundary condition with ", this->type(), "!");
      89         565 : }
      90             : 
      91             : Real
      92      521659 : LinearFVDiffusion::computeElemMatrixContribution()
      93             : {
      94      521659 :   return computeFluxMatrixContribution();
      95             : }
      96             : 
      97             : Real
      98      521659 : LinearFVDiffusion::computeNeighborMatrixContribution()
      99             : {
     100      521659 :   return -computeFluxMatrixContribution();
     101             : }
     102             : 
     103             : Real
     104      521659 : LinearFVDiffusion::computeElemRightHandSideContribution()
     105             : {
     106      521659 :   return computeFluxRHSContribution();
     107             : }
     108             : 
     109             : Real
     110      521659 : LinearFVDiffusion::computeNeighborRightHandSideContribution()
     111             : {
     112      521659 :   return -computeFluxRHSContribution();
     113             : }
     114             : 
     115             : Real
     116     1043318 : LinearFVDiffusion::computeFluxMatrixContribution()
     117             : {
     118             :   // If we don't have the value yet, we compute it
     119     1043318 :   if (!_cached_matrix_contribution)
     120             :   {
     121             :     // If we requested nonorthogonal correction, we use the normal component of the
     122             :     // cell to face vector.
     123      521659 :     const auto d = _use_nonorthogonal_correction
     124      521659 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
     125       96995 :                        : _current_face_info->dCNMag();
     126             : 
     127             :     // Cache the matrix contribution
     128      521659 :     _flux_matrix_contribution = faceDiffusivity() / d * _current_face_area;
     129      521659 :     _cached_matrix_contribution = true;
     130             :   }
     131             : 
     132     1043318 :   return _flux_matrix_contribution;
     133             : }
     134             : 
     135             : Real
     136     1043318 : LinearFVDiffusion::computeFluxRHSContribution()
     137             : {
     138             :   // We only have contributions on the right hand side from internal faces
     139             :   // if the nonorthogonal correction is enabled.
     140     1043318 :   if (_use_nonorthogonal_correction && !_cached_rhs_contribution)
     141             :   {
     142      424664 :     const auto state = determineState();
     143             : 
     144             :     // Get the gradients from the adjacent cells
     145      424664 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo(), state);
     146      424664 :     const auto & grad_neighbor = _var.gradSln(*_current_face_info->neighborInfo(), state);
     147             : 
     148             :     // Interpolate the two gradients to the face
     149             :     const auto interp_coeffs =
     150      424664 :         interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
     151             : 
     152             :     // Compute correction vector. Potential optimization: this only depends on the geometry
     153             :     // so we can cache it in FaceInfo at some point.
     154             :     const auto correction_vector =
     155      424664 :         _current_face_info->normal() -
     156      849328 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     157             : 
     158             :     // Cache the matrix contribution
     159      424664 :     _flux_rhs_contribution =
     160      424664 :         faceDiffusivity() *
     161      849328 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor) *
     162      424664 :         correction_vector * _current_face_area;
     163      424664 :     _cached_rhs_contribution = true;
     164             :   }
     165             : 
     166     1043318 :   return _flux_rhs_contribution;
     167             : }
     168             : 
     169             : Real
     170       50516 : LinearFVDiffusion::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     171             : {
     172       50516 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     173             :   mooseAssert(diff_bc, "This should be a valid BC!");
     174             : 
     175       50516 :   auto grad_contrib = diff_bc->computeBoundaryGradientMatrixContribution() * _current_face_area;
     176             :   // If the boundary condition does not include the diffusivity contribution then
     177             :   // add it here.
     178       50516 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     179             :   {
     180       37247 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     181       37247 :     grad_contrib *= _diffusion_coeff(face_arg, determineState());
     182             :   }
     183             : 
     184       50516 :   return grad_contrib;
     185             : }
     186             : 
     187             : Real
     188       50516 : LinearFVDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     189             : {
     190       50516 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     191             :   mooseAssert(diff_bc, "This should be a valid BC!");
     192             : 
     193       50516 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     194       50516 :   const auto state = determineState();
     195       50516 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution() * _current_face_area;
     196             : 
     197             :   // If the boundary condition does not include the diffusivity contribution then
     198             :   // add it here.
     199       50516 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     200       37247 :     grad_contrib *= _diffusion_coeff(face_arg, state);
     201             : 
     202             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     203             :   // this in the boundary condition too. For now, however, we keep it like this.
     204             :   // This should only be used for BCs where the gradient of the value is computed and
     205             :   // not prescribed.
     206             : 
     207       50516 :   if (_use_nonorthogonal_correction && diff_bc->useBoundaryGradientExtrapolation())
     208             :   {
     209             :     // We support internal boundaries as well. In that case we have to decide on which side
     210             :     // of the boundary we are on.
     211       28645 :     const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     212       28645 :                                ? _current_face_info->elemInfo()
     213           0 :                                : _current_face_info->neighborInfo();
     214       28645 :     const Real boundary_normal_multiplier =
     215       28645 :         (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     216             : 
     217             :     // Unit vector to the boundary. Unfortunately, we have to recompute it because the value
     218             :     // stored in the face info is only correct for external boundaries
     219       28645 :     const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
     220             :     const auto correction_vector =
     221       28645 :         _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
     222             : 
     223       28645 :     grad_contrib += _diffusion_coeff(face_arg, state) * _var.gradSln(*elem_info, state) *
     224       28645 :                     boundary_normal_multiplier * correction_vector * _current_face_area;
     225             :   }
     226             : 
     227       50516 :   return grad_contrib;
     228             : }