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       15727 : LinearFVDiffusion::validParams()
      19             : {
      20       15727 :   InputParameters params = LinearFVFluxKernel::validParams();
      21       15727 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23       47181 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25       31454 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       15727 :   params.addParam<MooseFunctorName>("diffusion_coeff", 1.0, "The diffusion coefficient.");
      28       15727 :   return params;
      29           0 : }
      30             : 
      31         731 : LinearFVDiffusion::LinearFVDiffusion(const InputParameters & params)
      32             :   : LinearFVFluxKernel(params),
      33         731 :     _diffusion_coeff(getFunctor<Real>("diffusion_coeff")),
      34         731 :     _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
      35         731 :     _flux_matrix_contribution(0.0),
      36         731 :     _flux_rhs_contribution(0.0)
      37             : {
      38         731 :   if (_use_nonorthogonal_correction)
      39         369 :     _var.computeCellGradients();
      40         731 : }
      41             : 
      42             : void
      43           0 : LinearFVDiffusion::initialSetup()
      44             : {
      45           0 :   for (const auto bc : _var.getBoundaryConditionMap())
      46           0 :     if (!dynamic_cast<const LinearFVAdvectionDiffusionBC *>(bc.second))
      47           0 :       mooseError(
      48           0 :           bc.second->type(), " is not a compatible boundary condition with ", this->type(), "!");
      49           0 : }
      50             : 
      51             : Real
      52      637940 : LinearFVDiffusion::computeElemMatrixContribution()
      53             : {
      54      637940 :   return computeFluxMatrixContribution();
      55             : }
      56             : 
      57             : Real
      58      637940 : LinearFVDiffusion::computeNeighborMatrixContribution()
      59             : {
      60      637940 :   return -computeFluxMatrixContribution();
      61             : }
      62             : 
      63             : Real
      64      637940 : LinearFVDiffusion::computeElemRightHandSideContribution()
      65             : {
      66      637940 :   return computeFluxRHSContribution();
      67             : }
      68             : 
      69             : Real
      70      637940 : LinearFVDiffusion::computeNeighborRightHandSideContribution()
      71             : {
      72      637940 :   return -computeFluxRHSContribution();
      73             : }
      74             : 
      75             : Real
      76     1275880 : LinearFVDiffusion::computeFluxMatrixContribution()
      77             : {
      78             :   // If we don't have the value yet, we compute it
      79     1275880 :   if (!_cached_matrix_contribution)
      80             :   {
      81      637940 :     const auto face_arg = makeCDFace(*_current_face_info);
      82             : 
      83             :     // If we requested nonorthogonal correction, we use the normal component of the
      84             :     // cell to face vector.
      85      637940 :     const auto d = _use_nonorthogonal_correction
      86      637940 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
      87      224124 :                        : _current_face_info->dCNMag();
      88             : 
      89             :     // Cache the matrix contribution
      90      637940 :     _flux_matrix_contribution =
      91      637940 :         _diffusion_coeff(face_arg, determineState()) / d * _current_face_area;
      92      637940 :     _cached_matrix_contribution = true;
      93             :   }
      94             : 
      95     1275880 :   return _flux_matrix_contribution;
      96             : }
      97             : 
      98             : Real
      99     1275880 : LinearFVDiffusion::computeFluxRHSContribution()
     100             : {
     101             :   // We only have contributions on the right hand side from internal faces
     102             :   // if the nonorthogonal correction is enabled.
     103     1275880 :   if (_use_nonorthogonal_correction && !_cached_rhs_contribution)
     104             :   {
     105      413816 :     const auto face_arg = makeCDFace(*_current_face_info);
     106      413816 :     const auto state_arg = determineState();
     107             : 
     108             :     // Get the gradients from the adjacent cells
     109      413816 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo());
     110      413816 :     const auto & grad_neighbor = _var.gradSln(*_current_face_info->neighborInfo());
     111             : 
     112             :     // Interpolate the two gradients to the face
     113             :     const auto interp_coeffs =
     114      413816 :         interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
     115             : 
     116             :     // Compute correction vector. Potential optimization: this only depends on the geometry
     117             :     // so we can cache it in FaceInfo at some point.
     118             :     const auto correction_vector =
     119      413816 :         _current_face_info->normal() -
     120      827632 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     121             : 
     122             :     // Cache the matrix contribution
     123      413816 :     _flux_rhs_contribution =
     124      413816 :         _diffusion_coeff(face_arg, state_arg) *
     125      827632 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor) *
     126      413816 :         correction_vector * _current_face_area;
     127      413816 :     _cached_rhs_contribution = true;
     128             :   }
     129             : 
     130     1275880 :   return _flux_rhs_contribution;
     131             : }
     132             : 
     133             : Real
     134       60658 : LinearFVDiffusion::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     135             : {
     136       60658 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     137             :   mooseAssert(diff_bc, "This should be a valid BC!");
     138             : 
     139       60658 :   auto grad_contrib = diff_bc->computeBoundaryGradientMatrixContribution() * _current_face_area;
     140             :   // If the boundary condition does not include the diffusivity contribution then
     141             :   // add it here.
     142       60658 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     143             :   {
     144       47484 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     145       47484 :     grad_contrib *= _diffusion_coeff(face_arg, determineState());
     146             :   }
     147             : 
     148       60658 :   return grad_contrib;
     149             : }
     150             : 
     151             : Real
     152       60658 : LinearFVDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     153             : {
     154       60658 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     155             :   mooseAssert(diff_bc, "This should be a valid BC!");
     156             : 
     157       60658 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     158       60658 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution() * _current_face_area;
     159             : 
     160             :   // If the boundary condition does not include the diffusivity contribution then
     161             :   // add it here.
     162       60658 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     163       47484 :     grad_contrib *= _diffusion_coeff(face_arg, determineState());
     164             : 
     165             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     166             :   // this in the boundary condition too. For now, however, we keep it like this.
     167             :   // This should only be used for BCs where the gradient of the value is computed and
     168             :   // not prescribed.
     169             : 
     170       60658 :   if (_use_nonorthogonal_correction && diff_bc->useBoundaryGradientExtrapolation())
     171             :   {
     172             :     const auto correction_vector =
     173       28666 :         _current_face_info->normal() -
     174       57332 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     175             : 
     176             :     // We might be on a face which is an internal boundary so we want to make sure we
     177             :     // get the gradient from the right side.
     178       28666 :     const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     179       28666 :                                ? _current_face_info->elemInfo()
     180           0 :                                : _current_face_info->neighborInfo();
     181             : 
     182       28666 :     grad_contrib += _diffusion_coeff(face_arg, determineState()) * _var.gradSln(*elem_info) *
     183       28666 :                     correction_vector * _current_face_area;
     184             :   }
     185             : 
     186       60658 :   return grad_contrib;
     187             : }