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 "PINSFVRhieChowInterpolator.h"
      11             : #include "Reconstructions.h"
      12             : #include "NS.h"
      13             : #include "Assembly.h"
      14             : #include "BernoulliPressureVariable.h"
      15             : 
      16             : registerMooseObject("NavierStokesApp", PINSFVRhieChowInterpolator);
      17             : 
      18             : InputParameters
      19        5890 : PINSFVRhieChowInterpolator::validParams()
      20             : {
      21        5890 :   auto params = INSFVRhieChowInterpolator::validParams();
      22        5890 :   params.addClassDescription("Performs interpolations and reconstructions of porosity and computes "
      23             :                              "the Rhie-Chow face velocities.");
      24        5890 :   ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", /*quiet=*/true);
      25        5890 :   exec_enum.setAdditionalValue(EXEC_INITIAL);
      26        5890 :   params.addRequiredParam<MooseFunctorName>(NS::porosity, "The porosity");
      27       17670 :   params.addParam<unsigned short>(
      28             :       "smoothing_layers",
      29       11780 :       0,
      30             :       "The number of interpolation-reconstruction operations to perform on the porosity");
      31       11780 :   params.addRelationshipManager("ElementSideNeighborLayers",
      32             :                                 Moose::RelationshipManagerType::GEOMETRIC,
      33        6311 :                                 [](const InputParameters & obj_params, InputParameters & rm_params)
      34             :                                 {
      35        6311 :                                   rm_params.set<unsigned short>("layers") =
      36       12622 :                                       obj_params.get<unsigned short>("smoothing_layers");
      37        6311 :                                   rm_params.set<bool>("use_displaced_mesh") =
      38        6311 :                                       obj_params.get<bool>("use_displaced_mesh");
      39        6311 :                                 });
      40        5890 :   return params;
      41           0 : }
      42             : 
      43        2946 : PINSFVRhieChowInterpolator::PINSFVRhieChowInterpolator(const InputParameters & params)
      44             :   : INSFVRhieChowInterpolator(params),
      45        5892 :     _eps(getFunctor<ADReal>(NS::porosity)),
      46        2946 :     _smoothed_eps(_moose_mesh, NS::smoothed_porosity, /*extrapolated_boundary*/ true),
      47        2946 :     _epss(libMesh::n_threads(), nullptr),
      48        2946 :     _smoothed_epss(libMesh::n_threads(), nullptr),
      49        8838 :     _smoothing_layers(getParam<unsigned short>("smoothing_layers"))
      50             : {
      51        2946 :   if (_smoothing_layers && _eps.wrapsType<MooseVariableBase>())
      52           2 :     paramError(
      53             :         NS::porosity,
      54             :         "If we are reconstructing porosity, then the input porosity to this user object cannot "
      55             :         "be a Moose variable. There are issues with reconstructing Moose variables: 1) initial "
      56             :         "conditions are run after user objects initial setup 2) reconstructing from a variable "
      57             :         "requires ghosting the solution vectors 3) it's difficult to restrict the face "
      58             :         "informations we evaluate interpolations and reconstructions on such that we never query "
      59             :         "an algebraically remote element due to things like two-term extrapolated boundary faces "
      60             :         "which trigger gradient evaluations which trigger neighbor element evaluation");
      61             : 
      62        2944 :   const auto porosity_name = deduceFunctorName(NS::porosity);
      63             : 
      64        6133 :   for (const auto tid : make_range(libMesh::n_threads()))
      65             :   {
      66        3189 :     _epss[tid] = &UserObject::_subproblem.getFunctor<ADReal>(porosity_name, tid, name(), true);
      67             : 
      68        3189 :     if (_smoothing_layers > 0)
      69             :     {
      70         593 :       if (!UserObject::_subproblem.hasFunctor(NS::smoothed_porosity, tid))
      71             :         // Smoothed porosity is only an envelope at this point, it will be set during pinsfvSetup()
      72         593 :         UserObject::_subproblem.addFunctor(NS::smoothed_porosity, _smoothed_eps, tid);
      73             : 
      74         593 :       _smoothed_epss[tid] =
      75         593 :           &UserObject::_subproblem.getFunctor<ADReal>(NS::smoothed_porosity, tid, name(), true);
      76             :     }
      77             :   }
      78        2944 : }
      79             : 
      80             : void
      81           0 : PINSFVRhieChowInterpolator::meshChanged()
      82             : {
      83           0 :   INSFVRhieChowInterpolator::meshChanged();
      84           0 :   pinsfvSetup();
      85           0 : }
      86             : 
      87             : void
      88        2566 : PINSFVRhieChowInterpolator::pinsfvSetup()
      89             : {
      90        2566 :   if (!_smoothing_layers)
      91             :     return;
      92             : 
      93         582 :   if (dynamic_cast<BernoulliPressureVariable *>(_p))
      94           0 :     paramError(
      95             :         NS::pressure,
      96             :         "If 'smoothing_layers' is non-zero, e.g. if the porosity is smooth(ed), "
      97             :         "then the pressure drop should be computed automatically. The "
      98             :         "'BernoulliPressureVariable' class enforces a pressure drop according to the Bernoulli "
      99             :         "equation on any face which has different pressure values on either side. This is "
     100             :         "undesirable when the porosity is varying smoothing and there may be pressure drops "
     101             :         "corresponding to viscous effects. Please just use the 'INSFVPressureVariable' "
     102             :         "class for pressure when the porosity is smooth.");
     103             : 
     104         582 :   const auto & all_fi = _moose_mesh.allFaceInfo();
     105         582 :   _geometric_fi.reserve(all_fi.size());
     106             : 
     107     3568618 :   for (const auto & fi : all_fi)
     108     3568036 :     if (isFaceGeometricallyRelevant(fi))
     109     3567904 :       _geometric_fi.push_back(&fi);
     110             : 
     111             :   _geometric_fi.shrink_to_fit();
     112             : 
     113         582 :   const auto saved_do_derivatives = ADReal::do_derivatives;
     114         582 :   ADReal::do_derivatives = true;
     115         582 :   Moose::FV::interpolateReconstruct(
     116         582 :       _smoothed_eps, _eps, _smoothing_layers, false, _geometric_fi, determineState());
     117         582 :   ADReal::do_derivatives = saved_do_derivatives;
     118             : 
     119             :   // Assign the new functor to all
     120         588 :   for (const auto tid : make_range((unsigned int)(1), libMesh::n_threads()))
     121             :   {
     122             :     auto & other_smoothed_epss = const_cast<Moose::Functor<ADReal> &>(
     123           6 :         UserObject::_subproblem.getFunctor<ADReal>(NS::smoothed_porosity, tid, name(), true));
     124           6 :     other_smoothed_epss.assign(_smoothed_eps);
     125             :   }
     126             : }
     127             : 
     128             : bool
     129     3568036 : PINSFVRhieChowInterpolator::isFaceGeometricallyRelevant(const FaceInfo & fi) const
     130             : {
     131     3568036 :   if (&fi.elem() == libMesh::remote_elem)
     132             :     return false;
     133             : 
     134     7136072 :   bool on_us = blockIDs().count(fi.elem().subdomain_id());
     135             : 
     136     3484572 :   if (fi.neighborPtr())
     137             :   {
     138     3484572 :     if (&fi.neighbor() == libMesh::remote_elem)
     139             :       return false;
     140             : 
     141     3484572 :     on_us = on_us || blockIDs().count(fi.neighbor().subdomain_id());
     142             :   }
     143             : 
     144       83464 :   if (!on_us)
     145             :     // Neither the element nor neighbor has a subdomain id on which we are active, so this face is
     146             :     // not relevant
     147           0 :     return false;
     148             : 
     149             :   //
     150             :   // Ok, we've established that either the element or neighbor is active on our subdomains and
     151             :   // neither of them are remote elements, so this face is still in the running to be considered
     152             :   // relevant. There is one more caveat to be considered. In the case that we are a boundary face,
     153             :   // we will generally need a two term expansion to compute our value, which will require a
     154             :   // cell-gradient evaluation. If that is the case, then all of our surrounding neighbors cannot be
     155             :   // remote. If we are not a boundary face, then at this point we're safe
     156             :   //
     157             : 
     158     3568036 :   if (!Moose::FV::onBoundary(blockIDs(), fi))
     159             :     return true;
     160             : 
     161           0 :   const auto & boundary_elem = (fi.neighborPtr() && blockIDs().count(fi.neighbor().subdomain_id()))
     162           0 :                                    ? fi.neighbor()
     163             :                                    : fi.elem();
     164             : 
     165      417056 :   for (auto * const neighbor : boundary_elem.neighbor_ptr_range())
     166             :   {
     167      333724 :     if (!neighbor)
     168       87982 :       continue;
     169             : 
     170      245742 :     if ((neighbor == libMesh::remote_elem))
     171         132 :       return false;
     172             :   }
     173             : 
     174             :   // We made it through all the tests!
     175       83332 :   return true;
     176             : }
     177             : 
     178             : void
     179       32019 : PINSFVRhieChowInterpolator::initialize()
     180             : {
     181       32019 :   if (_current_execute_flag == EXEC_INITIAL)
     182        2566 :     pinsfvSetup();
     183             :   else
     184             :     // Cannot compute Rhie Chow coefficients on initial
     185       29453 :     INSFVRhieChowInterpolator::initialize();
     186       32019 : }
     187             : 
     188             : void
     189       32019 : PINSFVRhieChowInterpolator::execute()
     190             : {
     191             :   // Cannot compute Rhie Chow coefficients on initial
     192       32019 :   if (_current_execute_flag != EXEC_INITIAL)
     193       29453 :     INSFVRhieChowInterpolator::execute();
     194       32019 : }
     195             : 
     196             : void
     197       32019 : PINSFVRhieChowInterpolator::finalize()
     198             : {
     199             :   // Cannot compute Rhie Chow coefficients on initial
     200       32019 :   if (_current_execute_flag != EXEC_INITIAL)
     201       29453 :     INSFVRhieChowInterpolator::finalize();
     202       32019 : }