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        3594 : PINSFVRhieChowInterpolator::validParams()
      20             : {
      21        3594 :   auto params = INSFVRhieChowInterpolator::validParams();
      22        3594 :   params.addClassDescription("Performs interpolations and reconstructions of porosity and computes "
      23             :                              "the Rhie-Chow face velocities.");
      24        3594 :   ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", /*quiet=*/true);
      25        3594 :   exec_enum.setAdditionalValue(EXEC_INITIAL);
      26        3594 :   params.addRequiredParam<MooseFunctorName>(NS::porosity, "The porosity");
      27       10782 :   params.addParam<unsigned short>(
      28             :       "smoothing_layers",
      29        7188 :       0,
      30             :       "The number of interpolation-reconstruction operations to perform on the porosity");
      31        7188 :   params.addRelationshipManager("ElementSideNeighborLayers",
      32             :                                 Moose::RelationshipManagerType::GEOMETRIC,
      33        4061 :                                 [](const InputParameters & obj_params, InputParameters & rm_params)
      34             :                                 {
      35        4061 :                                   rm_params.set<unsigned short>("layers") =
      36        8122 :                                       obj_params.get<unsigned short>("smoothing_layers");
      37        4061 :                                   rm_params.set<bool>("use_displaced_mesh") =
      38        4061 :                                       obj_params.get<bool>("use_displaced_mesh");
      39        4061 :                                 });
      40        3594 :   return params;
      41           0 : }
      42             : 
      43        1798 : PINSFVRhieChowInterpolator::PINSFVRhieChowInterpolator(const InputParameters & params)
      44             :   : INSFVRhieChowInterpolator(params),
      45        3596 :     _eps(getFunctor<ADReal>(NS::porosity)),
      46        1798 :     _smoothed_eps(_moose_mesh, NS::smoothed_porosity, /*extrapolated_boundary*/ true),
      47        1798 :     _epss(libMesh::n_threads(), nullptr),
      48        1798 :     _smoothed_epss(libMesh::n_threads(), nullptr),
      49        5394 :     _smoothing_layers(getParam<unsigned short>("smoothing_layers"))
      50             : {
      51        1798 :   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        1796 :   const auto porosity_name = deduceFunctorName(NS::porosity);
      63             : 
      64        3679 :   for (const auto tid : make_range(libMesh::n_threads()))
      65             :   {
      66        1883 :     _epss[tid] = &UserObject::_subproblem.getFunctor<ADReal>(porosity_name, tid, name(), true);
      67             : 
      68        1883 :     if (_smoothing_layers > 0)
      69             :     {
      70         353 :       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         353 :         UserObject::_subproblem.addFunctor(NS::smoothed_porosity, _smoothed_eps, tid);
      73             : 
      74         353 :       _smoothed_epss[tid] =
      75         353 :           &UserObject::_subproblem.getFunctor<ADReal>(NS::smoothed_porosity, tid, name(), true);
      76             :     }
      77             :   }
      78        1796 : }
      79             : 
      80             : void
      81           0 : PINSFVRhieChowInterpolator::meshChanged()
      82             : {
      83           0 :   INSFVRhieChowInterpolator::meshChanged();
      84           0 :   pinsfvSetup();
      85           0 : }
      86             : 
      87             : void
      88        1716 : PINSFVRhieChowInterpolator::pinsfvSetup()
      89             : {
      90        1716 :   if (!_smoothing_layers)
      91             :     return;
      92             : 
      93         350 :   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         350 :   const auto & all_fi = _moose_mesh.allFaceInfo();
     105         350 :   _geometric_fi.reserve(all_fi.size());
     106             : 
     107     2140850 :   for (const auto & fi : all_fi)
     108     2140500 :     if (isFaceGeometricallyRelevant(fi))
     109     2140456 :       _geometric_fi.push_back(&fi);
     110             : 
     111             :   _geometric_fi.shrink_to_fit();
     112             : 
     113         350 :   const auto saved_do_derivatives = ADReal::do_derivatives;
     114         350 :   ADReal::do_derivatives = true;
     115         350 :   Moose::FV::interpolateReconstruct(
     116         350 :       _smoothed_eps, _eps, _smoothing_layers, false, _geometric_fi, determineState());
     117         350 :   ADReal::do_derivatives = saved_do_derivatives;
     118             : 
     119             :   // Assign the new functor to all
     120         352 :   for (const auto tid : make_range((unsigned int)(1), libMesh::n_threads()))
     121             :   {
     122             :     auto & other_smoothed_epss = const_cast<Moose::Functor<ADReal> &>(
     123           2 :         UserObject::_subproblem.getFunctor<ADReal>(NS::smoothed_porosity, tid, name(), true));
     124           2 :     other_smoothed_epss.assign(_smoothed_eps);
     125             :   }
     126             : }
     127             : 
     128             : bool
     129     2140500 : PINSFVRhieChowInterpolator::isFaceGeometricallyRelevant(const FaceInfo & fi) const
     130             : {
     131     2140500 :   if (&fi.elem() == libMesh::remote_elem)
     132             :     return false;
     133             : 
     134     4281000 :   bool on_us = blockIDs().count(fi.elem().subdomain_id());
     135             : 
     136     2090444 :   if (fi.neighborPtr())
     137             :   {
     138     2090444 :     if (&fi.neighbor() == libMesh::remote_elem)
     139             :       return false;
     140             : 
     141     2090444 :     on_us = on_us || blockIDs().count(fi.neighbor().subdomain_id());
     142             :   }
     143             : 
     144       50056 :   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     2140500 :   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      250192 :   for (auto * const neighbor : boundary_elem.neighbor_ptr_range())
     166             :   {
     167      200180 :     if (!neighbor)
     168       52810 :       continue;
     169             : 
     170      147370 :     if ((neighbor == libMesh::remote_elem))
     171          44 :       return false;
     172             :   }
     173             : 
     174             :   // We made it through all the tests!
     175       50012 :   return true;
     176             : }
     177             : 
     178             : void
     179       21312 : PINSFVRhieChowInterpolator::initialize()
     180             : {
     181       21312 :   if (_current_execute_flag == EXEC_INITIAL)
     182        1716 :     pinsfvSetup();
     183             :   else
     184             :     // Cannot compute Rhie Chow coefficients on initial
     185       19596 :     INSFVRhieChowInterpolator::initialize();
     186       21312 : }
     187             : 
     188             : void
     189       21312 : PINSFVRhieChowInterpolator::execute()
     190             : {
     191             :   // Cannot compute Rhie Chow coefficients on initial
     192       21312 :   if (_current_execute_flag != EXEC_INITIAL)
     193       19596 :     INSFVRhieChowInterpolator::execute();
     194       21312 : }
     195             : 
     196             : void
     197       21312 : PINSFVRhieChowInterpolator::finalize()
     198             : {
     199             :   // Cannot compute Rhie Chow coefficients on initial
     200       21312 :   if (_current_execute_flag != EXEC_INITIAL)
     201       19596 :     INSFVRhieChowInterpolator::finalize();
     202       21312 : }