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 "ComputeLinearFVLimitedGradientThread.h"
      11             : 
      12             : #include "GradientLimiterType.h"
      13             : #include "SystemBase.h"
      14             : #include "PetscVectorReader.h"
      15             : #include "FEProblemBase.h"
      16             : #include "FVUtils.h"
      17             : 
      18             : #include "libmesh/dof_object.h"
      19             : 
      20             : #include <algorithm>
      21             : #include <cmath>
      22             : #include <limits>
      23             : 
      24       11870 : ComputeLinearFVLimitedGradientThread::ComputeLinearFVLimitedGradientThread(
      25             :     FEProblemBase & fe_problem,
      26             :     SystemBase & system,
      27             :     const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_gradient,
      28             :     std::vector<std::unique_ptr<NumericVector<Number>>> & temporary_limited_gradient,
      29             :     const Moose::FV::GradientLimiterType limiter_type,
      30       11870 :     const std::unordered_set<unsigned int> & requested_variables)
      31       11870 :   : _fe_problem(fe_problem),
      32       11870 :     _dim(_fe_problem.mesh().dimension()),
      33       11870 :     _system(system),
      34       11870 :     _libmesh_system(system.system()),
      35       11870 :     _system_number(_libmesh_system.number()),
      36       11870 :     _raw_gradient(raw_gradient),
      37       11870 :     _limiter_type(limiter_type),
      38       11870 :     _requested_variables(requested_variables),
      39       11870 :     _temporary_limited_gradient(temporary_limited_gradient)
      40             : {
      41       11870 : }
      42             : 
      43           0 : ComputeLinearFVLimitedGradientThread::ComputeLinearFVLimitedGradientThread(
      44           0 :     ComputeLinearFVLimitedGradientThread & x, Threads::split /*split*/)
      45           0 :   : _fe_problem(x._fe_problem),
      46           0 :     _dim(x._dim),
      47           0 :     _system(x._system),
      48           0 :     _libmesh_system(x._libmesh_system),
      49           0 :     _system_number(x._system_number),
      50           0 :     _raw_gradient(x._raw_gradient),
      51           0 :     _limiter_type(x._limiter_type),
      52           0 :     _requested_variables(x._requested_variables),
      53           0 :     _temporary_limited_gradient(x._temporary_limited_gradient)
      54             : {
      55           0 : }
      56             : 
      57             : void
      58       11870 : ComputeLinearFVLimitedGradientThread::operator()(const ElemInfoRange & range)
      59             : {
      60       11870 :   ParallelUniqueId puid;
      61       11870 :   _tid = puid.id;
      62             : 
      63       11870 :   const auto & raw_grad_container = _raw_gradient;
      64             : 
      65       11870 :   if (_limiter_type != Moose::FV::GradientLimiterType::Venkatakrishnan)
      66           0 :     mooseError("ComputeLinearFVLimitedGradientThread currently supports only the Venkatakrishnan "
      67             :                "limiter.");
      68             : 
      69       11870 :   unsigned int size = 0;
      70             : 
      71       23740 :   for (const auto & variable : _system.getVariables(_tid))
      72             :   {
      73       11870 :     _current_var = dynamic_cast<MooseLinearVariableFV<Real> *>(variable);
      74       11870 :     if (!_current_var)
      75           0 :       continue;
      76             : 
      77       11870 :     if (!_current_var->needsGradientVectorStorage())
      78           0 :       continue;
      79             : 
      80       11870 :     if (!_requested_variables.count(_current_var->number()))
      81           0 :       continue;
      82             : 
      83       11870 :     if (!size)
      84       11870 :       size = range.size();
      85             : 
      86       11870 :     std::vector<std::vector<Real>> temporary_values(_temporary_limited_gradient.size(),
      87       11870 :                                                     std::vector<Real>(size, 0.0));
      88       11870 :     std::vector<dof_id_type> dof_indices(size, 0);
      89             : 
      90       11870 :     PetscVectorReader solution_reader(*_libmesh_system.current_local_solution);
      91       11870 :     std::vector<PetscVectorReader> grad_reader;
      92       11870 :     grad_reader.reserve(raw_grad_container.size());
      93       32580 :     for (const auto dim_index : index_range(raw_grad_container))
      94       20710 :       grad_reader.emplace_back(*raw_grad_container[dim_index]);
      95             : 
      96             :     mooseAssert(raw_grad_container.size() >= _dim,
      97             :                 "Raw gradient container has fewer components than mesh dimension.");
      98             :     mooseAssert(_temporary_limited_gradient.size() >= _dim,
      99             :                 "Limited gradient container has fewer components than mesh dimension.");
     100             : 
     101       11870 :     auto elem_iterator = range.begin();
     102    16459982 :     for (const auto elem_i : make_range(size))
     103             :     {
     104    16448112 :       const auto & elem_info = *elem_iterator;
     105    16448112 :       elem_iterator++;
     106             : 
     107    16448112 :       if (!_current_var->hasBlocks(elem_info->subdomain_id()))
     108     1010592 :         continue;
     109             : 
     110    16448112 :       const dof_id_type dof = elem_info->dofIndices()[_system_number][_current_var->number()];
     111    16448112 :       if (dof == libMesh::DofObject::invalid_id)
     112           0 :         continue;
     113             : 
     114    16448112 :       dof_indices[elem_i] = dof;
     115             : 
     116    16448112 :       const Real phi_elem = solution_reader(dof);
     117    16448112 :       Real max_value = phi_elem;
     118    16448112 :       Real min_value = phi_elem;
     119             : 
     120             :       // Gather one-ring min/max values.
     121    16448112 :       const Elem * const elem = elem_info->elem();
     122    73216896 :       for (const auto side : make_range(elem->n_sides()))
     123             :       {
     124    56768784 :         const Elem * const neighbor = elem->neighbor_ptr(side);
     125    56768784 :         if (!neighbor)
     126      735804 :           continue;
     127             : 
     128    56032980 :         const auto & neighbor_info = _fe_problem.mesh().elemInfo(neighbor->id());
     129    56032980 :         if (!_current_var->hasBlocks(neighbor_info.subdomain_id()))
     130           0 :           continue;
     131             : 
     132             :         const dof_id_type neighbor_dof =
     133    56032980 :             neighbor_info.dofIndices()[_system_number][_current_var->number()];
     134    56032980 :         if (neighbor_dof == libMesh::DofObject::invalid_id)
     135           0 :           continue;
     136             : 
     137    56032980 :         const Real phi_neighbor = solution_reader(neighbor_dof);
     138    56032980 :         max_value = std::max(max_value, phi_neighbor);
     139    56032980 :         min_value = std::min(min_value, phi_neighbor);
     140             :       }
     141             : 
     142             :       // Read the raw cell gradient.
     143    16448112 :       VectorValue<Real> raw_grad;
     144    16448112 :       raw_grad.zero();
     145    49306764 :       for (const auto dim_index : make_range(_dim))
     146    32858652 :         raw_grad(dim_index) = grad_reader[dim_index](dof);
     147             : 
     148             :       // If the stencil is constant (or nearly constant), don't attempt to limit.
     149    16448112 :       if (std::abs(max_value - min_value) < 1e-14)
     150             :       {
     151     3031590 :         for (const auto dim_index : make_range(_dim))
     152     2020998 :           temporary_values[dim_index][elem_i] = raw_grad(dim_index);
     153     1010592 :         continue;
     154     1010592 :       }
     155             : 
     156    15437520 :       Real alpha = 1.0;
     157    15437520 :       const Point & elem_centroid = elem_info->centroid();
     158             : 
     159    69076382 :       for (const auto side : make_range(elem->n_sides()))
     160             :       {
     161    53638862 :         const Elem * const neighbor = elem->neighbor_ptr(side);
     162    53638862 :         if (!neighbor)
     163      670203 :           continue;
     164             : 
     165    52968659 :         const auto & neighbor_info = _fe_problem.mesh().elemInfo(neighbor->id());
     166    52968659 :         if (!_current_var->hasBlocks(neighbor_info.subdomain_id()))
     167           0 :           continue;
     168             : 
     169             :         const dof_id_type neighbor_dof =
     170    52968659 :             neighbor_info.dofIndices()[_system_number][_current_var->number()];
     171    52968659 :         if (neighbor_dof == libMesh::DofObject::invalid_id)
     172           0 :           continue;
     173             : 
     174    52968659 :         const bool elem_has_face_info = Moose::FV::elemHasFaceInfo(*elem, neighbor);
     175    52968659 :         const Elem * const fi_elem = elem_has_face_info ? elem : neighbor;
     176             :         const unsigned int fi_side =
     177    52968659 :             elem_has_face_info ? side : neighbor->which_neighbor_am_i(elem);
     178    52968659 :         const auto * fi = _fe_problem.mesh().faceInfo(fi_elem, fi_side);
     179             :         mooseAssert(fi,
     180             :                     "Missing FaceInfo for neighboring elements with centroid " +
     181             :                         Moose::stringify(elem_info->centroid()) + " and " +
     182             :                         Moose::stringify(neighbor->vertex_average()) +
     183             :                         " while computing limited gradients.");
     184             : 
     185    52968659 :         const Point face_point = fi->faceCentroid();
     186             : 
     187    52968659 :         const Real delta_face = raw_grad * (face_point - elem_centroid);
     188             : 
     189    52968659 :         Real h = elem->hmin();
     190    52968659 :         Real grad_mag = raw_grad.norm();
     191             : 
     192    52968659 :         Real eps = 0.1 * (grad_mag * h) * (grad_mag * h) + 1e-20;
     193             : 
     194    52968659 :         const Real delta_max = std::abs(max_value - phi_elem) + eps;
     195    52968659 :         const Real delta_min = std::abs(min_value - phi_elem) + eps;
     196             : 
     197    79652239 :         const Real rf = (delta_face >= 0.0) ? std::abs(delta_face) / delta_max
     198    26683580 :                                             : std::abs(delta_face) / delta_min;
     199             : 
     200    52968659 :         const Real beta = (2.0 * rf + 1.0) / (rf * (2.0 * rf + 1.0) + 1.0);
     201    52968659 :         alpha = std::min(alpha, beta);
     202             :       }
     203             : 
     204    15437520 :       const VectorValue<Real> limited_grad = alpha * raw_grad;
     205    46275174 :       for (const auto dim_index : make_range(_dim))
     206    30837654 :         temporary_values[dim_index][elem_i] = limited_grad(dim_index);
     207             :     }
     208             : 
     209       32580 :     for (const auto dim_index : make_range(_dim))
     210             :     {
     211       20710 :       _temporary_limited_gradient[dim_index]->add_vector(temporary_values[dim_index].data(),
     212             :                                                          dof_indices);
     213             :     }
     214       11870 :   }
     215       11870 : }
     216             : 
     217             : void
     218           0 : ComputeLinearFVLimitedGradientThread::join(const ComputeLinearFVLimitedGradientThread & /*y*/)
     219             : {
     220           0 : }