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       14405 : LinearFVAnisotropicDiffusion::validParams()
      19             : {
      20       14405 :   InputParameters params = LinearFVFluxKernel::validParams();
      21       14405 :   params.addClassDescription("Represents the matrix and right hand side contributions of a "
      22             :                              "diffusion term in a partial differential equation.");
      23       43215 :   params.addParam<bool>(
      24             :       "use_nonorthogonal_correction",
      25       28810 :       true,
      26             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      27       14405 :   params.addParam<bool>(
      28             :       "use_nonorthogonal_correction_on_boundary",
      29             :       "If the nonorthogonal correction should be used when computing the normal gradient.");
      30       14405 :   params.addRequiredParam<MooseFunctorName>("diffusion_tensor",
      31             :                                             "Functor describing a diagonal diffusion tensor.");
      32       14405 :   return params;
      33           0 : }
      34             : 
      35          70 : LinearFVAnisotropicDiffusion::LinearFVAnisotropicDiffusion(const InputParameters & params)
      36             :   : LinearFVFluxKernel(params),
      37          70 :     _diffusion_tensor(getFunctor<RealVectorValue>("diffusion_tensor")),
      38          70 :     _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
      39          70 :     _use_nonorthogonal_correction_on_boundary(
      40         140 :         isParamValid("use_nonorthogonal_correction_on_boundary")
      41         140 :             ? getParam<bool>("use_nonorthogonal_correction_on_boundary")
      42          70 :             : _use_nonorthogonal_correction),
      43          70 :     _flux_matrix_contribution(0.0),
      44          70 :     _flux_rhs_contribution(0.0)
      45             : {
      46          70 :   _var.computeCellGradients();
      47          70 : }
      48             : 
      49             : void
      50           0 : LinearFVAnisotropicDiffusion::initialSetup()
      51             : {
      52           0 :   for (const auto bc : _var.getBoundaryConditionMap())
      53           0 :     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           0 : }
      57             : 
      58             : Real
      59      168175 : LinearFVAnisotropicDiffusion::computeElemMatrixContribution()
      60             : {
      61      168175 :   return computeFluxMatrixContribution();
      62             : }
      63             : 
      64             : Real
      65      168175 : LinearFVAnisotropicDiffusion::computeNeighborMatrixContribution()
      66             : {
      67      168175 :   return -computeFluxMatrixContribution();
      68             : }
      69             : 
      70             : Real
      71      168175 : LinearFVAnisotropicDiffusion::computeElemRightHandSideContribution()
      72             : {
      73      168175 :   return computeFluxRHSContribution();
      74             : }
      75             : 
      76             : Real
      77      168175 : LinearFVAnisotropicDiffusion::computeNeighborRightHandSideContribution()
      78             : {
      79      168175 :   return -computeFluxRHSContribution();
      80             : }
      81             : 
      82             : Real
      83      336350 : LinearFVAnisotropicDiffusion::computeFluxMatrixContribution()
      84             : {
      85             :   // If we don't have the value yet, we compute it
      86      336350 :   if (!_cached_matrix_contribution)
      87             :   {
      88      168175 :     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      168175 :     const auto d = _use_nonorthogonal_correction
      93      168175 :                        ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
      94       17794 :                        : _current_face_info->dCNMag();
      95             : 
      96      168175 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, determineState());
      97             : 
      98      672700 :     for (const auto i : make_range(Moose::dim))
      99      504525 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     100             : 
     101      168175 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     102             : 
     103             :     // Cache the matrix contribution
     104      168175 :     _flux_matrix_contribution = normal_scaled_diff_tensor / d * _current_face_area;
     105      168175 :     _cached_matrix_contribution = true;
     106             :   }
     107             : 
     108      336350 :   return _flux_matrix_contribution;
     109             : }
     110             : 
     111             : Real
     112      336350 : LinearFVAnisotropicDiffusion::computeFluxRHSContribution()
     113             : {
     114             :   // Cache the RHS contribution
     115      336350 :   if (!_cached_rhs_contribution)
     116             :   {
     117      168175 :     const auto face_arg = makeCDFace(*_current_face_info);
     118      168175 :     const auto state_arg = determineState();
     119             : 
     120             :     // Get the gradients from the adjacent cells
     121      168175 :     const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo());
     122      168175 :     const auto grad_neighbor = _var.gradSln(*_current_face_info->neighborInfo());
     123             : 
     124             :     // Interpolate the two gradients to the face
     125             :     const auto interp_coeffs =
     126      168175 :         interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
     127             : 
     128             :     const auto interpolated_gradient =
     129      168175 :         (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor);
     130             : 
     131      168175 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, state_arg);
     132             : 
     133      672700 :     for (const auto i : make_range(Moose::dim))
     134      504525 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     135             : 
     136      168175 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     137             : 
     138      168175 :     _flux_rhs_contribution =
     139      168175 :         (scaled_diff_tensor - normal_scaled_diff_tensor * _current_face_info->normal()) *
     140             :         interpolated_gradient;
     141             : 
     142      168175 :     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      150381 :           _current_face_info->normal() -
     148      150381 :           1 / (_current_face_info->normal() * _current_face_info->eCN()) *
     149      451143 :               _current_face_info->eCN();
     150             : 
     151      150381 :       _flux_rhs_contribution +=
     152      150381 :           normal_scaled_diff_tensor * interpolated_gradient * correction_vector;
     153             :     }
     154      168175 :     _flux_rhs_contribution *= _current_face_area;
     155      168175 :     _cached_rhs_contribution = true;
     156             :   }
     157             : 
     158      336350 :   return _flux_rhs_contribution;
     159             : }
     160             : 
     161             : Real
     162       16926 : LinearFVAnisotropicDiffusion::computeBoundaryMatrixContribution(
     163             :     const LinearFVBoundaryCondition & bc)
     164             : {
     165       16926 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     166             :   mooseAssert(diff_bc, "This should be a valid BC!");
     167             : 
     168       16926 :   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       16926 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     172             :   {
     173       16926 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
     174             : 
     175       16926 :     auto scaled_diff_tensor = _diffusion_tensor(face_arg, determineState());
     176             : 
     177       67704 :     for (const auto i : make_range(Moose::dim))
     178       50778 :       scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     179             : 
     180       16926 :     auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     181             : 
     182       16926 :     grad_contrib *= normal_scaled_diff_tensor;
     183             :   }
     184             : 
     185       16926 :   return grad_contrib;
     186             : }
     187             : 
     188             : Real
     189       16926 : LinearFVAnisotropicDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     190             : {
     191       16926 :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     192             :   mooseAssert(diff_bc, "This should be a valid BC!");
     193             : 
     194       16926 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     195       16926 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution();
     196             : 
     197       16926 :   auto scaled_diff_tensor = _diffusion_tensor(face_arg, determineState());
     198             : 
     199       67704 :   for (const auto i : make_range(Moose::dim))
     200       50778 :     scaled_diff_tensor(i) = _current_face_info->normal()(i) * scaled_diff_tensor(i);
     201             : 
     202       16926 :   auto normal_scaled_diff_tensor = scaled_diff_tensor * _current_face_info->normal();
     203       16926 :   auto boundary_grad = _var.gradSln(*_current_face_info->elemInfo());
     204             : 
     205             :   // If the boundary condition does not include the diffusivity contribution then
     206             :   // add it here.
     207       16926 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     208             :   {
     209       16926 :     grad_contrib *= normal_scaled_diff_tensor;
     210             :   }
     211             : 
     212       16926 :   grad_contrib += (scaled_diff_tensor - normal_scaled_diff_tensor * _current_face_info->normal()) *
     213             :                   boundary_grad;
     214             : 
     215             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     216             :   // this in the boundary condition too. For now, however, we keep it like this.
     217       16926 :   if (_use_nonorthogonal_correction_on_boundary)
     218             :   {
     219             :     const auto correction_vector =
     220       14322 :         _current_face_info->normal() -
     221       28644 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     222             : 
     223       14322 :     grad_contrib += normal_scaled_diff_tensor * boundary_grad * correction_vector;
     224             :   }
     225             : 
     226       16926 :   return grad_contrib * _current_face_area;
     227             : }