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       15945 : LinearFVDiffusion::validParams()
      19             : {
      20       15945 :   InputParameters params = LinearFVFluxKernel::validParams();
      21       15945 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23       47835 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25       31890 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       15945 :   params.addParam<MooseFunctorName>("diffusion_coeff", 1.0, "The diffusion coefficient.");
      28       15945 :   return params;
      29           0 : }
      30             : 
      31         840 : LinearFVDiffusion::LinearFVDiffusion(const InputParameters & params)
      32             :   : LinearFVFluxKernel(params),
      33         840 :     _diffusion_coeff(getFunctor<Real>("diffusion_coeff")),
      34         840 :     _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
      35         840 :     _flux_matrix_contribution(0.0),
      36         840 :     _flux_rhs_contribution(0.0)
      37             : {
      38         840 :   if (_use_nonorthogonal_correction)
      39         449 :     _var.computeCellGradients();
      40         840 : }
      41             : 
      42             : void
      43         866 : LinearFVDiffusion::initialSetup()
      44             : {
      45        3303 :   for (const auto bc : _var.getBoundaryConditionMap())
      46        2437 :     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         866 : }
      50             : 
      51             : Real
      52      685792 : LinearFVDiffusion::computeElemMatrixContribution()
      53             : {
      54      685792 :   return computeFluxMatrixContribution();
      55             : }
      56             : 
      57             : Real
      58      685792 : LinearFVDiffusion::computeNeighborMatrixContribution()
      59             : {
      60      685792 :   return -computeFluxMatrixContribution();
      61             : }
      62             : 
      63             : Real
      64      685792 : LinearFVDiffusion::computeElemRightHandSideContribution()
      65             : {
      66      685792 :   return computeFluxRHSContribution();
      67             : }
      68             : 
      69             : Real
      70      685792 : LinearFVDiffusion::computeNeighborRightHandSideContribution()
      71             : {
      72      685792 :   return -computeFluxRHSContribution();
      73             : }
      74             : 
      75             : Real
      76     1371584 : LinearFVDiffusion::computeFluxMatrixContribution()
      77             : {
      78             :   // If we don't have the value yet, we compute it
      79     1371584 :   if (!_cached_matrix_contribution)
      80             :   {
      81      685792 :     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      685792 :     const auto d = _use_nonorthogonal_correction
      86      685792 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
      87      224346 :                        : _current_face_info->dCNMag();
      88             : 
      89             :     // Cache the matrix contribution
      90      685792 :     _flux_matrix_contribution =
      91      685792 :         _diffusion_coeff(face_arg, determineState()) / d * _current_face_area;
      92      685792 :     _cached_matrix_contribution = true;
      93             :   }
      94             : 
      95     1371584 :   return _flux_matrix_contribution;
      96             : }
      97             : 
      98             : Real
      99     1371584 : LinearFVDiffusion::computeFluxRHSContribution()
     100             : {
     101             :   // We only have contributions on the right hand side from internal faces
     102             :   // if the nonorthogonal correction is enabled.
     103     1371584 :   if (_use_nonorthogonal_correction && !_cached_rhs_contribution)
     104             :   {
     105      461446 :     const auto face_arg = makeCDFace(*_current_face_info);
     106      461446 :     const auto state_arg = determineState();
     107             : 
     108             :     // Get the gradients from the adjacent cells
     109      461446 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo());
     110      461446 :     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      461446 :         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      461446 :         _current_face_info->normal() -
     120      922892 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     121             : 
     122             :     // Cache the matrix contribution
     123      461446 :     _flux_rhs_contribution =
     124      461446 :         _diffusion_coeff(face_arg, state_arg) *
     125      922892 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor) *
     126      461446 :         correction_vector * _current_face_area;
     127      461446 :     _cached_rhs_contribution = true;
     128             :   }
     129             : 
     130     1371584 :   return _flux_rhs_contribution;
     131             : }
     132             : 
     133             : Real
     134       69492 : LinearFVDiffusion::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     135             : {
     136       69492 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     137             :   mooseAssert(diff_bc, "This should be a valid BC!");
     138             : 
     139       69492 :   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       69492 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     143             :   {
     144       56318 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     145       56318 :     grad_contrib *= _diffusion_coeff(face_arg, determineState());
     146             :   }
     147             : 
     148       69492 :   return grad_contrib;
     149             : }
     150             : 
     151             : Real
     152       69492 : LinearFVDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     153             : {
     154       69492 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     155             :   mooseAssert(diff_bc, "This should be a valid BC!");
     156             : 
     157       69492 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     158       69492 :   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       69492 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     163       56318 :     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       69492 :   if (_use_nonorthogonal_correction && diff_bc->useBoundaryGradientExtrapolation())
     171             :   {
     172             :     // We support internal boundaries as well. In that case we have to decide on which side
     173             :     // of the boundary we are on.
     174       35426 :     const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     175       35426 :                                ? _current_face_info->elemInfo()
     176           0 :                                : _current_face_info->neighborInfo();
     177       35426 :     const Real boundary_normal_multiplier =
     178       35426 :         (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     179             : 
     180             :     // Unit vector to the boundary. Unfortunately, we have to recompute it because the value
     181             :     // stored in the face info is only correct for external boundaries
     182       35426 :     const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
     183             :     const auto correction_vector =
     184       35426 :         _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
     185             : 
     186       35426 :     grad_contrib += _diffusion_coeff(face_arg, determineState()) * _var.gradSln(*elem_info) *
     187       35426 :                     boundary_normal_multiplier * correction_vector * _current_face_area;
     188             :   }
     189             : 
     190       69492 :   return grad_contrib;
     191             : }