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