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 "ComputeLinearFVGreenGaussGradientVolumeThread.h"
      11             : #include "LinearSystem.h"
      12             : #include "LinearFVBoundaryCondition.h"
      13             : #include "PetscVectorReader.h"
      14             : #include "FEProblemBase.h"
      15             : 
      16        9865 : ComputeLinearFVGreenGaussGradientVolumeThread::ComputeLinearFVGreenGaussGradientVolumeThread(
      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             : {
      25        9865 : }
      26             : 
      27           0 : ComputeLinearFVGreenGaussGradientVolumeThread::ComputeLinearFVGreenGaussGradientVolumeThread(
      28           0 :     ComputeLinearFVGreenGaussGradientVolumeThread & x, Threads::split /*split*/)
      29           0 :   : _fe_problem(x._fe_problem),
      30           0 :     _dim(x._dim),
      31           0 :     _linear_system_number(x._linear_system_number),
      32           0 :     _linear_system(x._linear_system),
      33           0 :     _system_number(x._system_number)
      34             : {
      35           0 : }
      36             : 
      37             : void
      38        9865 : ComputeLinearFVGreenGaussGradientVolumeThread::operator()(const ElemInfoRange & range)
      39             : {
      40        9865 :   ParallelUniqueId puid;
      41        9865 :   _tid = puid.id;
      42             : 
      43        9865 :   auto & linear_system = _fe_problem.getLinearSystem(_linear_system_number);
      44             : 
      45             :   // This will be the vector we work on since the old gradient might still be needed
      46             :   // to compute extrapolated boundary conditions for example.
      47        9865 :   auto & grad_container = linear_system.newGradientContainer();
      48             : 
      49             :   // Computing this size can be very expensive so we only want to do it once
      50        9865 :   unsigned int size = 0;
      51             : 
      52       19730 :   for (const auto & variable : linear_system.getVariables(_tid))
      53             :   {
      54        9865 :     _current_var = dynamic_cast<MooseLinearVariableFV<Real> *>(variable);
      55             :     mooseAssert(_current_var,
      56             :                 "This should be a linear FV variable, did we somehow add a nonlinear variable to "
      57             :                 "the linear system?");
      58        9865 :     if (_current_var->needsGradientVectorStorage())
      59             :     {
      60        8038 :       if (!size)
      61        8038 :         size = range.size();
      62             : 
      63        8038 :       const auto rz_radial_coord = _fe_problem.mesh().getAxisymmetricRadialCoord();
      64        8038 :       const auto state = Moose::currentState();
      65             : 
      66             :       std::vector<std::vector<Real>> new_values(grad_container.size(),
      67        8038 :                                                 std::vector<Real>(size, 0.0));
      68        8038 :       std::vector<dof_id_type> dof_indices(size, 0);
      69             :       {
      70        8038 :         std::vector<PetscVectorReader> grad_reader;
      71       21914 :         for (const auto dim_index : index_range(grad_container))
      72       13876 :           grad_reader.emplace_back(*grad_container[dim_index]);
      73             : 
      74             :         // Iterate over all the elements in the range
      75        8038 :         auto elem_iterator = range.begin();
      76     2495694 :         for (const auto elem_i : make_range(size))
      77             :         {
      78     2487656 :           const auto & elem_info = *elem_iterator;
      79     2487656 :           if (_current_var->hasBlocks(elem_info->subdomain_id()))
      80             :           {
      81     2487656 :             const auto coord_type = _fe_problem.mesh().getCoordSystem(elem_info->subdomain_id());
      82             : 
      83             :             mooseAssert(coord_type != Moose::CoordinateSystemType::COORD_RSPHERICAL,
      84             :                         "We have not yet implemented the correct translation from gradient to "
      85             :                         "divergence for "
      86             :                         "spherical coordinates yet.");
      87             : 
      88     2487656 :             dof_indices[elem_i] = elem_info->dofIndices()[_system_number][_current_var->number()];
      89     2487656 :             const auto volume = elem_info->volume() * elem_info->coordFactor();
      90             : 
      91     7424192 :             for (const auto dim_index : index_range(grad_container))
      92     4936536 :               new_values[dim_index][elem_i] = grad_reader[dim_index](dof_indices[elem_i]) / volume;
      93             : 
      94     2487656 :             if (coord_type == Moose::CoordinateSystemType::COORD_RZ)
      95             :             {
      96      439516 :               const auto radial_contrib = _current_var->getElemValue(*elem_info, state) /
      97      439516 :                                           elem_info->centroid()(rz_radial_coord);
      98      439516 :               new_values[rz_radial_coord][elem_i] += radial_contrib;
      99             :             }
     100             :           }
     101     2487656 :           elem_iterator++;
     102             :         }
     103        8038 :       }
     104       21914 :       for (const auto dim_index : index_range(grad_container))
     105             :       {
     106       13876 :         grad_container[dim_index]->zero();
     107       13876 :         grad_container[dim_index]->add_vector(new_values[dim_index].data(), dof_indices);
     108             :       }
     109        8038 :     }
     110             :   }
     111        9865 : }
     112             : 
     113             : void
     114           0 : ComputeLinearFVGreenGaussGradientVolumeThread::join(
     115             :     const ComputeLinearFVGreenGaussGradientVolumeThread & /*y*/)
     116             : {
     117           0 : }