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        4179 : LinearFVDiffusion::validParams()
      19             : {
      20        4179 :   InputParameters params = LinearFVFluxKernel::validParams();
      21        8358 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23       12537 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25        8358 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       16716 :   params.addParam<MooseFunctorName>("diffusion_coeff", 1.0, "The diffusion coefficient.");
      28       12537 :   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        4179 :   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      788873 : LinearFVDiffusion::faceDiffusivity() const
      55             : {
      56      788873 :   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      418885 :     const auto state = determineState();
      62             : 
      63      418885 :     if (_coeff_interp_method)
      64           0 :       _face_diffusivity =
      65           0 :           _coeff_interp_method->interpolate(_diffusion_coeff, *_current_face_info, state);
      66             :     else
      67      418885 :       _face_diffusivity = _diffusion_coeff(makeCDFace(*_current_face_info), state);
      68             : 
      69      418885 :     _cached_face_diffusivity = true;
      70             :   }
      71             : 
      72      788873 :   return _face_diffusivity;
      73             : }
      74             : 
      75             : void
      76      465498 : LinearFVDiffusion::setupFaceData(const FaceInfo * face_info)
      77             : {
      78      465498 :   LinearFVFluxKernel::setupFaceData(face_info);
      79      465498 :   _cached_face_diffusivity = false;
      80      465498 : }
      81             : 
      82             : void
      83         543 : LinearFVDiffusion::initialSetup()
      84             : {
      85        2040 :   for (const auto bc : _var.getBoundaryConditionMap())
      86        1497 :     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         543 : }
      90             : 
      91             : Real
      92      418885 : LinearFVDiffusion::computeElemMatrixContribution()
      93             : {
      94      418885 :   return computeFluxMatrixContribution();
      95             : }
      96             : 
      97             : Real
      98      418885 : LinearFVDiffusion::computeNeighborMatrixContribution()
      99             : {
     100      418885 :   return -computeFluxMatrixContribution();
     101             : }
     102             : 
     103             : Real
     104      418885 : LinearFVDiffusion::computeElemRightHandSideContribution()
     105             : {
     106      418885 :   return computeFluxRHSContribution();
     107             : }
     108             : 
     109             : Real
     110      418885 : LinearFVDiffusion::computeNeighborRightHandSideContribution()
     111             : {
     112      418885 :   return -computeFluxRHSContribution();
     113             : }
     114             : 
     115             : Real
     116      837770 : LinearFVDiffusion::computeFluxMatrixContribution()
     117             : {
     118             :   // If we don't have the value yet, we compute it
     119      837770 :   if (!_cached_matrix_contribution)
     120             :   {
     121             :     // If we requested nonorthogonal correction, we use the normal component of the
     122             :     // cell to face vector.
     123      418885 :     const auto d = _use_nonorthogonal_correction
     124      418885 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
     125       48897 :                        : _current_face_info->dCNMag();
     126             : 
     127             :     // Cache the matrix contribution
     128      418885 :     _flux_matrix_contribution = faceDiffusivity() / d * _current_face_area;
     129      418885 :     _cached_matrix_contribution = true;
     130             :   }
     131             : 
     132      837770 :   return _flux_matrix_contribution;
     133             : }
     134             : 
     135             : Real
     136      837770 : LinearFVDiffusion::computeFluxRHSContribution()
     137             : {
     138             :   // We only have contributions on the right hand side from internal faces
     139             :   // if the nonorthogonal correction is enabled.
     140      837770 :   if (_use_nonorthogonal_correction && !_cached_rhs_contribution)
     141             :   {
     142      369988 :     const auto state = determineState();
     143             : 
     144             :     // Get the gradients from the adjacent cells
     145      369988 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo(), state);
     146      369988 :     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      369988 :         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      369988 :         _current_face_info->normal() -
     156      739976 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     157             : 
     158             :     // Cache the matrix contribution
     159      369988 :     _flux_rhs_contribution =
     160      369988 :         faceDiffusivity() *
     161      739976 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor) *
     162      369988 :         correction_vector * _current_face_area;
     163      369988 :     _cached_rhs_contribution = true;
     164             :   }
     165             : 
     166      837770 :   return _flux_rhs_contribution;
     167             : }
     168             : 
     169             : Real
     170       40773 : LinearFVDiffusion::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     171             : {
     172       40773 :   const auto * const diff_bc = libMesh::cast_ptr<const LinearFVAdvectionDiffusionBC *>(&bc);
     173             : 
     174       40773 :   auto grad_contrib = diff_bc->computeBoundaryGradientMatrixContribution() * _current_face_area;
     175             :   // If the boundary condition does not include the diffusivity contribution then
     176             :   // add it here.
     177       40773 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     178             :   {
     179       29577 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     180       29577 :     grad_contrib *= _diffusion_coeff(face_arg, determineState());
     181             :   }
     182             : 
     183       40773 :   return grad_contrib;
     184             : }
     185             : 
     186             : Real
     187       40773 : LinearFVDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     188             : {
     189       40773 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     190             :   mooseAssert(diff_bc, "This should be a valid BC!");
     191             : 
     192       40773 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     193       40773 :   const auto state = determineState();
     194       40773 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution() * _current_face_area;
     195             : 
     196             :   // If the boundary condition does not include the diffusivity contribution then
     197             :   // add it here.
     198       40773 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     199       29577 :     grad_contrib *= _diffusion_coeff(face_arg, state);
     200             : 
     201             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     202             :   // this in the boundary condition too. For now, however, we keep it like this.
     203             :   // This should only be used for BCs where the gradient of the value is computed and
     204             :   // not prescribed.
     205             : 
     206       40773 :   if (_use_nonorthogonal_correction && diff_bc->useBoundaryGradientExtrapolation())
     207             :   {
     208             :     // We support internal boundaries as well. In that case we have to decide on which side
     209             :     // of the boundary we are on.
     210       24972 :     const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     211       24972 :                                ? _current_face_info->elemInfo()
     212           0 :                                : _current_face_info->neighborInfo();
     213       24972 :     const Real boundary_normal_multiplier =
     214       24972 :         (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     215             : 
     216             :     // Unit vector to the boundary. Unfortunately, we have to recompute it because the value
     217             :     // stored in the face info is only correct for external boundaries
     218       24972 :     const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
     219             :     const auto correction_vector =
     220       24972 :         _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
     221             : 
     222       24972 :     grad_contrib += _diffusion_coeff(face_arg, state) * _var.gradSln(*elem_info, state) *
     223       24972 :                     boundary_normal_multiplier * correction_vector * _current_face_area;
     224             :   }
     225             : 
     226       40773 :   return grad_contrib;
     227             : }