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 "LinearFVAnisotropicDiffusion.h"
      11             : #include "Assembly.h"
      12             : #include "SubProblem.h"
      13             : #include "LinearFVAdvectionDiffusionBC.h"
      14             : 
      15             : registerMooseObject("MooseApp", LinearFVAnisotropicDiffusion);
      16             : 
      17             : InputParameters
      18        3121 : LinearFVAnisotropicDiffusion::validParams()
      19             : {
      20        3121 :   InputParameters params = LinearFVFluxKernel::validParams();
      21        6242 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23        9363 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25        6242 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       12484 :   params.addParam<bool>(
      28             :       "use_nonorthogonal_correction_on_boundary",
      29             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      30        9363 :   params.addRequiredParam<MooseFunctorName>("diffusion_tensor",
      31             :                                             "Functor describing a diagonal diffusion tensor.");
      32        3121 :   return params;
      33           0 : }
      34             : 
      35          30 : LinearFVAnisotropicDiffusion::LinearFVAnisotropicDiffusion(const InputParameters & params)
      36             :   : LinearFVFluxKernel(params),
      37          30 :     _diffusion_tensor(getFunctor<RealVectorValue>("diffusion_tensor")),
      38          60 :     _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
      39          30 :     _use_nonorthogonal_correction_on_boundary(
      40          60 :         isParamValid("use_nonorthogonal_correction_on_boundary")
      41          60 :             ? getParam<bool>("use_nonorthogonal_correction_on_boundary")
      42          30 :             : _use_nonorthogonal_correction),
      43          30 :     _flux_matrix_contribution(0.0),
      44          30 :     _flux_rhs_contribution(0.0)
      45             : {
      46          30 :   _var.computeCellGradients();
      47          30 : }
      48             : 
      49             : void
      50          30 : LinearFVAnisotropicDiffusion::initialSetup()
      51             : {
      52         150 :   for (const auto bc : _var.getBoundaryConditionMap())
      53         120 :     if (!dynamic_cast<const LinearFVAdvectionDiffusionBC *>(bc.second))
      54           0 :       mooseError(
      55           0 :           bc.second->type(), " is not a compatible boundary condition with ", this->type(), "!");
      56          30 : }
      57             : 
      58             : Real
      59       72075 : LinearFVAnisotropicDiffusion::computeElemMatrixContribution()
      60             : {
      61       72075 :   return computeFluxMatrixContribution();
      62             : }
      63             : 
      64             : Real
      65       72075 : LinearFVAnisotropicDiffusion::computeNeighborMatrixContribution()
      66             : {
      67       72075 :   return -computeFluxMatrixContribution();
      68             : }
      69             : 
      70             : Real
      71       72075 : LinearFVAnisotropicDiffusion::computeElemRightHandSideContribution()
      72             : {
      73       72075 :   return computeFluxRHSContribution();
      74             : }
      75             : 
      76             : Real
      77       72075 : LinearFVAnisotropicDiffusion::computeNeighborRightHandSideContribution()
      78             : {
      79       72075 :   return -computeFluxRHSContribution();
      80             : }
      81             : 
      82             : Real
      83      144150 : LinearFVAnisotropicDiffusion::computeFluxMatrixContribution()
      84             : {
      85             :   // If we don't have the value yet, we compute it
      86      144150 :   if (!_cached_matrix_contribution)
      87             :   {
      88       72075 :     const auto face_arg = makeCDFace(*_current_face_info);
      89             : 
      90             :     // If we requested nonorthogonal correction, we use the normal component of the
      91             :     // cell to face vector.
      92       72075 :     const auto d = _use_nonorthogonal_correction
      93       72075 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
      94        7626 :                        : _current_face_info->dCNMag();
      95             : 
      96       72075 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, determineState());
      97             : 
      98      288300 :     for (const auto i : make_range(Moose::dim))
      99      216225 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     100             : 
     101       72075 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     102             : 
     103             :     // Cache the matrix contribution
     104       72075 :     _flux_matrix_contribution = normal_scaled_diff_tensor / d * _current_face_area;
     105       72075 :     _cached_matrix_contribution = true;
     106             :   }
     107             : 
     108      144150 :   return _flux_matrix_contribution;
     109             : }
     110             : 
     111             : Real
     112      144150 : LinearFVAnisotropicDiffusion::computeFluxRHSContribution()
     113             : {
     114             :   // Cache the RHS contribution
     115      144150 :   if (!_cached_rhs_contribution)
     116             :   {
     117       72075 :     const auto face_arg = makeCDFace(*_current_face_info);
     118       72075 :     const auto state_arg = determineState();
     119             : 
     120             :     // Get the gradients from the adjacent cells
     121       72075 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo(), state_arg);
     122       72075 :     const auto grad_neighbor = _var.gradSln(*_current_face_info->neighborInfo(), state_arg);
     123             : 
     124             :     // Interpolate the two gradients to the face
     125             :     const auto interp_coeffs =
     126       72075 :         interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
     127             : 
     128             :     const auto interpolated_gradient =
     129       72075 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor);
     130             : 
     131       72075 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, state_arg);
     132             : 
     133      288300 :     for (const auto i : make_range(Moose::dim))
     134      216225 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     135             : 
     136       72075 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     137             : 
     138       72075 :     _flux_rhs_contribution =
     139       72075 :         (scaled_diff_tensor - normal_scaled_diff_tensor * _current_face_info->normal()) *
     140             :         interpolated_gradient;
     141             : 
     142       72075 :     if (_use_nonorthogonal_correction)
     143             :     {
     144             :       // Compute correction vector. Potential optimization: this only depends on the geometry
     145             :       // so we can cache it in FaceInfo at some point.
     146             :       const auto correction_vector =
     147       64449 :           _current_face_info->normal() -
     148       64449 :           1 / (_current_face_info->normal() * _current_face_info->eCN()) *
     149      193347 :               _current_face_info->eCN();
     150             : 
     151       64449 :       _flux_rhs_contribution +=
     152       64449 :           normal_scaled_diff_tensor * interpolated_gradient * correction_vector;
     153             :     }
     154       72075 :     _flux_rhs_contribution *= _current_face_area;
     155       72075 :     _cached_rhs_contribution = true;
     156             :   }
     157             : 
     158      144150 :   return _flux_rhs_contribution;
     159             : }
     160             : 
     161             : Real
     162        7254 : LinearFVAnisotropicDiffusion::computeBoundaryMatrixContribution(
     163             :     const LinearFVBoundaryCondition & bc)
     164             : {
     165        7254 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     166             :   mooseAssert(diff_bc, "This should be a valid BC!");
     167             : 
     168        7254 :   auto grad_contrib = diff_bc->computeBoundaryGradientMatrixContribution() * _current_face_area;
     169             :   // If the boundary condition does not include the diffusivity contribution then
     170             :   // add it here.
     171        7254 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     172             :   {
     173        7254 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     174             : 
     175        7254 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, determineState());
     176             : 
     177       29016 :     for (const auto i : make_range(Moose::dim))
     178       21762 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     179             : 
     180        7254 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     181             : 
     182        7254 :     grad_contrib *= normal_scaled_diff_tensor;
     183             :   }
     184             : 
     185        7254 :   return grad_contrib;
     186             : }
     187             : 
     188             : Real
     189        7254 : LinearFVAnisotropicDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     190             : {
     191        7254 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     192             :   mooseAssert(diff_bc, "This should be a valid BC!");
     193             : 
     194        7254 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     195        7254 :   const auto state_arg = determineState();
     196        7254 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution();
     197             : 
     198        7254 :   auto scaled_diff_tensor = _diffusion_tensor(face_arg, state_arg);
     199             : 
     200       29016 :   for (const auto i : make_range(Moose::dim))
     201       21762 :     scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     202             : 
     203        7254 :   auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     204        7254 :   const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     205        7254 :                              ? _current_face_info->elemInfo()
     206           0 :                              : _current_face_info->neighborInfo();
     207             :   mooseAssert(elem_info, "We should always have an element info for the current face");
     208             : 
     209        7254 :   auto boundary_grad = _var.gradSln(*elem_info, state_arg);
     210             : 
     211             :   // If the boundary condition does not include the diffusivity contribution then
     212             :   // add it here.
     213        7254 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     214        7254 :     grad_contrib *= normal_scaled_diff_tensor;
     215             : 
     216             :   // We allow internal boundaries as well, in that case we have to make sure the normals point in
     217             :   // the right direction
     218        7254 :   const Real boundary_normal_multiplier =
     219        7254 :       (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     220             : 
     221           0 :   grad_contrib += (scaled_diff_tensor - normal_scaled_diff_tensor * boundary_normal_multiplier *
     222        7254 :                                             _current_face_info->normal()) *
     223             :                   boundary_grad;
     224             : 
     225             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     226             :   // this in the boundary condition too. For now, however, we keep it like this.
     227        7254 :   if (diff_bc->useBoundaryGradientExtrapolation() && _use_nonorthogonal_correction_on_boundary)
     228             :   {
     229        6138 :     const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
     230             :     const auto correction_vector =
     231        6138 :         _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
     232             : 
     233        6138 :     grad_contrib +=
     234        6138 :         normal_scaled_diff_tensor * boundary_grad * boundary_normal_multiplier * correction_vector;
     235             :   }
     236             : 
     237        7254 :   return grad_contrib * _current_face_area;
     238             : }