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 "LinearSystem.h"
      12             : #include "LinearFVBoundaryCondition.h"
      13             : #include "PetscVectorReader.h"
      14             : #include "FEProblemBase.h"
      15             : 
      16       10678 : ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFaceThread(
      17       10678 :     FEProblemBase & fe_problem, const unsigned int linear_system_num)
      18       10678 :   : _fe_problem(fe_problem),
      19       10678 :     _dim(_fe_problem.mesh().dimension()),
      20       10678 :     _linear_system_number(linear_system_num),
      21       10678 :     _linear_system(libMesh::cast_ref<libMesh::LinearImplicitSystem &>(
      22       10678 :         _fe_problem.getLinearSystem(_linear_system_number).system())),
      23       10678 :     _system_number(_linear_system.number()),
      24       10678 :     _new_gradient(_fe_problem.getLinearSystem(_linear_system_number).newGradientContainer())
      25             : {
      26       10678 : }
      27             : 
      28           0 : ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFaceThread(
      29           0 :     ComputeLinearFVGreenGaussGradientFaceThread & x, Threads::split /*split*/)
      30           0 :   : _fe_problem(x._fe_problem),
      31           0 :     _dim(x._dim),
      32           0 :     _linear_system_number(x._linear_system_number),
      33           0 :     _linear_system(x._linear_system),
      34           0 :     _system_number(x._system_number),
      35             :     // This will be the vector we work on since the old gradient might still be needed
      36             :     // to compute extrapolated boundary conditions for example.
      37           0 :     _new_gradient(x._new_gradient)
      38             : {
      39           0 : }
      40             : 
      41             : void
      42       10678 : ComputeLinearFVGreenGaussGradientFaceThread::operator()(const FaceInfoRange & range)
      43             : {
      44       10678 :   ParallelUniqueId puid;
      45       10678 :   _tid = puid.id;
      46             : 
      47       10678 :   unsigned int size = 0;
      48             : 
      49       10678 :   for (const auto & variable :
      50       32034 :        _fe_problem.getLinearSystem(_linear_system_number).getVariables(_tid))
      51             :   {
      52       10678 :     _current_var = dynamic_cast<MooseLinearVariableFV<Real> *>(variable);
      53             :     mooseAssert(_current_var,
      54             :                 "This should be a linear FV variable, did we somehow add a nonlinear variable to "
      55             :                 "the linear system?");
      56       10678 :     if (_current_var->needsGradientVectorStorage())
      57             :     {
      58        9320 :       if (!size)
      59        9320 :         size = range.size();
      60             : 
      61        9320 :       std::vector<std::vector<Real>> new_values_elem(_new_gradient.size(),
      62        9320 :                                                      std::vector<Real>(size, 0.0));
      63        9320 :       std::vector<std::vector<Real>> new_values_neighbor(_new_gradient.size(),
      64        9320 :                                                          std::vector<Real>(size, 0.0));
      65        9320 :       std::vector<dof_id_type> dof_indices_elem(size, 0);
      66        9320 :       std::vector<dof_id_type> dof_indices_neighbor(size, 0);
      67             : 
      68             :       {
      69        9320 :         PetscVectorReader solution_reader(*_linear_system.current_local_solution);
      70             : 
      71             :         // Iterate over all the face infos in the range
      72        9320 :         auto face_iterator = range.begin();
      73     4700710 :         for (const auto & face_i : make_range(size))
      74             :         {
      75     4691390 :           const auto & face_info = *face_iterator;
      76             : 
      77             :           const auto current_face_type =
      78     4691390 :               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     4691390 :           if (current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      83             :           {
      84     4349474 :             dof_indices_elem[face_i] =
      85     4349474 :                 face_info->elemInfo()->dofIndices()[_system_number][_current_var->number()];
      86     4349474 :             dof_indices_neighbor[face_i] =
      87     4349474 :                 face_info->neighborInfo()->dofIndices()[_system_number][_current_var->number()];
      88             : 
      89             :             const auto face_value =
      90     4349474 :                 Moose::FV::linearInterpolation(solution_reader(dof_indices_elem[face_i]),
      91     4349474 :                                                solution_reader(dof_indices_neighbor[face_i]),
      92             :                                                *face_info,
      93     4349474 :                                                true);
      94             : 
      95             :             const auto contribution =
      96     4349474 :                 face_info->normal() * face_info->faceArea() * face_info->faceCoord() * face_value;
      97             : 
      98    13007268 :             for (const auto i : make_range(_dim))
      99             :             {
     100     8657794 :               new_values_elem[i][face_i] = contribution(i);
     101     8657794 :               new_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      341916 :           else if (current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     108             :                    current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     109             :           {
     110             :             auto * bc_pointer =
     111      341916 :                 _current_var->getBoundaryCondition(*face_info->boundaryIDs().begin());
     112             : 
     113      341916 :             if (bc_pointer)
     114      326956 :               bc_pointer->setupFaceData(face_info, current_face_type);
     115             : 
     116             :             const auto * const elem_info = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     117      341916 :                                                ? 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      341916 :             const auto multiplier =
     123      341916 :                 current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0;
     124      341916 :             auto & dof_id_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     125             :                                           ? dof_indices_elem
     126             :                                           : dof_indices_neighbor;
     127      341916 :             auto & contribution_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     128             :                                                 ? new_values_elem
     129             :                                                 : new_values_neighbor;
     130             : 
     131      341916 :             dof_id_container[face_i] =
     132      341916 :                 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      341916 :             const auto contribution = multiplier * face_info->normal() * face_info->faceArea() *
     137      683832 :                                       face_info->faceCoord() *
     138      356876 :                                       (bc_pointer ? bc_pointer->computeBoundaryValue()
     139      356876 :                                                   : solution_reader(dof_id_container[face_i]));
     140     1020984 :             for (const auto i : make_range(_dim))
     141      679068 :               contribution_container[i][face_i] = contribution(i);
     142             :           }
     143     4691390 :           face_iterator++;
     144             :         }
     145        9320 :       }
     146       25430 :       for (const auto i : make_range(_dim))
     147             :       {
     148       16110 :         _new_gradient[i]->add_vector(new_values_elem[i].data(), dof_indices_elem);
     149       16110 :         _new_gradient[i]->add_vector(new_values_neighbor[i].data(), dof_indices_neighbor);
     150             :       }
     151        9320 :     }
     152             :   }
     153       10678 : }
     154             : 
     155             : void
     156           0 : ComputeLinearFVGreenGaussGradientFaceThread::join(
     157             :     const ComputeLinearFVGreenGaussGradientFaceThread & /*y*/)
     158             : {
     159           0 : }