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 "ComputeLinearFVGreenGaussGradientFaceThread.h"
      11             : #include "LinearFVBoundaryCondition.h"
      12             : #include "SystemBase.h"
      13             : #include "PetscVectorReader.h"
      14             : #include "FEProblemBase.h"
      15             : 
      16       43577 : ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFaceThread(
      17             :     FEProblemBase & fe_problem,
      18             :     SystemBase & system,
      19       43577 :     std::vector<std::unique_ptr<NumericVector<Number>>> & temporary_gradient)
      20       43577 :   : _fe_problem(fe_problem),
      21       43577 :     _dim(_fe_problem.mesh().dimension()),
      22       43577 :     _system(system),
      23       43577 :     _libmesh_system(system.system()),
      24       43577 :     _system_number(_libmesh_system.number()),
      25       43577 :     _temporary_gradient(temporary_gradient)
      26             : {
      27       43577 : }
      28             : 
      29           0 : ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFaceThread(
      30           0 :     ComputeLinearFVGreenGaussGradientFaceThread & x, Threads::split /*split*/)
      31           0 :   : _fe_problem(x._fe_problem),
      32           0 :     _dim(x._dim),
      33           0 :     _system(x._system),
      34           0 :     _libmesh_system(x._libmesh_system),
      35           0 :     _system_number(x._system_number),
      36             :     // This will be the vector we work on since the old gradient might still be needed
      37             :     // to compute extrapolated boundary conditions for example.
      38           0 :     _temporary_gradient(x._temporary_gradient)
      39             : {
      40           0 : }
      41             : 
      42             : void
      43       43577 : ComputeLinearFVGreenGaussGradientFaceThread::operator()(const FaceInfoRange & range)
      44             : {
      45       43577 :   ParallelUniqueId puid;
      46       43577 :   _tid = puid.id;
      47             : 
      48       43577 :   unsigned int size = 0;
      49             : 
      50       87154 :   for (const auto & variable : _system.getVariables(_tid))
      51             :   {
      52       43577 :     _current_var = dynamic_cast<MooseLinearVariableFV<Real> *>(variable);
      53       43577 :     if (!_current_var)
      54           0 :       continue;
      55             : 
      56       43577 :     if (_current_var->needsGradientVectorStorage())
      57             :     {
      58       43577 :       if (!size)
      59       43577 :         size = range.size();
      60             : 
      61       43577 :       std::vector<std::vector<Real>> temporary_values_elem(_temporary_gradient.size(),
      62       43577 :                                                            std::vector<Real>(size, 0.0));
      63       43577 :       std::vector<std::vector<Real>> temporary_values_neighbor(_temporary_gradient.size(),
      64       43577 :                                                                std::vector<Real>(size, 0.0));
      65       43577 :       std::vector<dof_id_type> dof_indices_elem(size, 0);
      66       43577 :       std::vector<dof_id_type> dof_indices_neighbor(size, 0);
      67             : 
      68             :       {
      69       43577 :         PetscVectorReader solution_reader(*_libmesh_system.current_local_solution);
      70             : 
      71             :         // Iterate over all the face infos in the range
      72       43577 :         auto face_iterator = range.begin();
      73    89531106 :         for (const auto & face_i : make_range(size))
      74             :         {
      75    89487529 :           const auto & face_info = *face_iterator;
      76             : 
      77             :           const auto current_face_type =
      78    89487529 :               face_info->faceType(std::make_pair(_current_var->number(), _system_number));
      79             : 
      80             :           // First we check if this face is internal to the variable, if yes, contribute to both
      81             :           // sides
      82    89487529 :           if (current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      83             :           {
      84    86972563 :             dof_indices_elem[face_i] =
      85    86972563 :                 face_info->elemInfo()->dofIndices()[_system_number][_current_var->number()];
      86    86972563 :             dof_indices_neighbor[face_i] =
      87    86972563 :                 face_info->neighborInfo()->dofIndices()[_system_number][_current_var->number()];
      88             : 
      89             :             const auto face_value =
      90    86972563 :                 Moose::FV::linearInterpolation(solution_reader(dof_indices_elem[face_i]),
      91    86972563 :                                                solution_reader(dof_indices_neighbor[face_i]),
      92             :                                                *face_info,
      93    86972563 :                                                true);
      94             : 
      95             :             const auto contribution =
      96    86972563 :                 face_info->normal() * face_info->faceArea() * face_info->faceCoord() * face_value;
      97             : 
      98   260816096 :             for (const auto i : make_range(_dim))
      99             :             {
     100   173843533 :               temporary_values_elem[i][face_i] = contribution(i);
     101   173843533 :               temporary_values_neighbor[i][face_i] = -contribution(i);
     102             :             }
     103             :           }
     104             :           // If this face is on the boundary of the block where the variable is defined, we
     105             :           // check for boundary conditions. If we don't find any we use an automatic one-term
     106             :           // expansion to compute the face value.
     107     2514966 :           else if (current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     108             :                    current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     109             :           {
     110             :             auto * bc_pointer =
     111     2514966 :                 _current_var->getBoundaryCondition(*face_info->boundaryIDs().begin());
     112             : 
     113     2514966 :             if (bc_pointer)
     114     2503286 :               bc_pointer->setupFaceData(face_info, current_face_type);
     115             : 
     116             :             const auto * const elem_info = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     117     2514966 :                                                ? face_info->elemInfo()
     118           0 :                                                : face_info->neighborInfo();
     119             : 
     120             :             // We have to account for cases when this face is an internal boundary and the normal
     121             :             // points in the wrong direction
     122     2514966 :             const auto multiplier =
     123     2514966 :                 current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0;
     124     2514966 :             auto & dof_id_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     125             :                                           ? dof_indices_elem
     126             :                                           : dof_indices_neighbor;
     127     2514966 :             auto & contribution_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     128             :                                                 ? temporary_values_elem
     129             :                                                 : temporary_values_neighbor;
     130             : 
     131     2514966 :             dof_id_container[face_i] =
     132     2514966 :                 elem_info->dofIndices()[_system_number][_current_var->number()];
     133             : 
     134             :             // If we don't have a boundary condition, then it's a natural condition. We'll use a
     135             :             // one-term expansion approximation in that case
     136     2514966 :             const auto contribution = multiplier * face_info->normal() * face_info->faceArea() *
     137     5029932 :                                       face_info->faceCoord() *
     138     2526646 :                                       (bc_pointer ? bc_pointer->computeBoundaryValue()
     139     2526646 :                                                   : solution_reader(dof_id_container[face_i]));
     140     7526680 :             for (const auto i : make_range(_dim))
     141     5011714 :               contribution_container[i][face_i] = contribution(i);
     142             :           }
     143    89487529 :           face_iterator++;
     144             :         }
     145       43577 :       }
     146      121291 :       for (const auto i : make_range(_dim))
     147             :       {
     148       77714 :         _temporary_gradient[i]->add_vector(temporary_values_elem[i].data(), dof_indices_elem);
     149       77714 :         _temporary_gradient[i]->add_vector(temporary_values_neighbor[i].data(),
     150             :                                            dof_indices_neighbor);
     151             :       }
     152       43577 :     }
     153             :   }
     154       43577 : }
     155             : 
     156             : void
     157           0 : ComputeLinearFVGreenGaussGradientFaceThread::join(
     158             :     const ComputeLinearFVGreenGaussGradientFaceThread & /*y*/)
     159             : {
     160           0 : }