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        9865 : ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFaceThread(
      17        9865 :     FEProblemBase & fe_problem, const unsigned int linear_system_num)
      18        9865 :   : _fe_problem(fe_problem),
      19        9865 :     _dim(_fe_problem.mesh().dimension()),
      20        9865 :     _linear_system_number(linear_system_num),
      21        9865 :     _linear_system(libMesh::cast_ref<libMesh::LinearImplicitSystem &>(
      22        9865 :         _fe_problem.getLinearSystem(_linear_system_number).system())),
      23        9865 :     _system_number(_linear_system.number()),
      24        9865 :     _new_gradient(_fe_problem.getLinearSystem(_linear_system_number).newGradientContainer())
      25             : {
      26        9865 : }
      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        9865 : ComputeLinearFVGreenGaussGradientFaceThread::operator()(const FaceInfoRange & range)
      43             : {
      44        9865 :   ParallelUniqueId puid;
      45        9865 :   _tid = puid.id;
      46             : 
      47        9865 :   unsigned int size = 0;
      48             : 
      49        9865 :   for (const auto & variable :
      50       29595 :        _fe_problem.getLinearSystem(_linear_system_number).getVariables(_tid))
      51             :   {
      52        9865 :     _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        9865 :     if (_current_var->needsGradientVectorStorage())
      57             :     {
      58        8038 :       if (!size)
      59        8038 :         size = range.size();
      60             : 
      61        8038 :       std::vector<std::vector<Real>> new_values_elem(_new_gradient.size(),
      62        8038 :                                                      std::vector<Real>(size, 0.0));
      63        8038 :       std::vector<std::vector<Real>> new_values_neighbor(_new_gradient.size(),
      64        8038 :                                                          std::vector<Real>(size, 0.0));
      65        8038 :       std::vector<dof_id_type> dof_indices_elem(size, 0);
      66        8038 :       std::vector<dof_id_type> dof_indices_neighbor(size, 0);
      67             : 
      68             :       {
      69        8038 :         PetscVectorReader solution_reader(*_linear_system.current_local_solution);
      70             : 
      71             :         // Iterate over all the face infos in the range
      72        8038 :         auto face_iterator = range.begin();
      73     4360482 :         for (const auto & face_i : make_range(size))
      74             :         {
      75     4352444 :           const auto & face_info = *face_iterator;
      76             : 
      77             :           const auto current_face_type =
      78     4352444 :               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     4352444 :           if (current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      83             :           {
      84     4049620 :             dof_indices_elem[face_i] =
      85     4049620 :                 face_info->elemInfo()->dofIndices()[_system_number][_current_var->number()];
      86     4049620 :             dof_indices_neighbor[face_i] =
      87     4049620 :                 face_info->neighborInfo()->dofIndices()[_system_number][_current_var->number()];
      88             : 
      89             :             const auto face_value =
      90     4049620 :                 Moose::FV::linearInterpolation(solution_reader(dof_indices_elem[face_i]),
      91     4049620 :                                                solution_reader(dof_indices_neighbor[face_i]),
      92             :                                                *face_info,
      93     4049620 :                                                true);
      94             : 
      95             :             const auto contribution =
      96     4049620 :                 face_info->normal() * face_info->faceArea() * face_info->faceCoord() * face_value;
      97             : 
      98    12112136 :             for (const auto i : make_range(_dim))
      99             :             {
     100     8062516 :               new_values_elem[i][face_i] = contribution(i);
     101     8062516 :               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      302824 :           else if (current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     108             :                    current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     109             :           {
     110             :             auto * bc_pointer =
     111      302824 :                 _current_var->getBoundaryCondition(*face_info->boundaryIDs().begin());
     112             : 
     113      302824 :             if (bc_pointer)
     114      289704 :               bc_pointer->setupFaceData(face_info, current_face_type);
     115             : 
     116             :             const auto * const elem_info = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     117      302824 :                                                ? 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      302824 :             const auto multiplier =
     123      302824 :                 current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0;
     124      302824 :             auto & dof_id_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     125             :                                           ? dof_indices_elem
     126             :                                           : dof_indices_neighbor;
     127      302824 :             auto & contribution_container = current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     128             :                                                 ? new_values_elem
     129             :                                                 : new_values_neighbor;
     130             : 
     131      302824 :             dof_id_container[face_i] =
     132      302824 :                 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      302824 :             const auto contribution = multiplier * face_info->normal() * face_info->faceArea() *
     137      605648 :                                       face_info->faceCoord() *
     138      315944 :                                       (bc_pointer ? bc_pointer->computeBoundaryValue()
     139      315944 :                                                   : solution_reader(dof_id_container[face_i]));
     140      904368 :             for (const auto i : make_range(_dim))
     141      601544 :               contribution_container[i][face_i] = contribution(i);
     142             :           }
     143     4352444 :           face_iterator++;
     144             :         }
     145        8038 :       }
     146       21914 :       for (const auto i : make_range(_dim))
     147             :       {
     148       13876 :         _new_gradient[i]->add_vector(new_values_elem[i].data(), dof_indices_elem);
     149       13876 :         _new_gradient[i]->add_vector(new_values_neighbor[i].data(), dof_indices_neighbor);
     150             :       }
     151        8038 :     }
     152             :   }
     153        9865 : }
     154             : 
     155             : void
     156           0 : ComputeLinearFVGreenGaussGradientFaceThread::join(
     157             :     const ComputeLinearFVGreenGaussGradientFaceThread & /*y*/)
     158             : {
     159           0 : }