Move more vector operations out of the threaded loop. (idaholab#25722)

framework/include/loops/ComputeLinearFVGreenGaussGradientFaceThread.h

@@ -93,5 +93,5 @@ class ComputeLinearFVGreenGaussGradientFaceThread
 
   /// Cache for the new gradient which is being built. It is needed because in certain scenarios the
   /// old gradient is used for computing the new gradient.
-  std::vector<std::unique_ptr<NumericVector<Number>>> _new_gradient;
+  std::vector<std::unique_ptr<NumericVector<Number>>> & _new_gradient;
 };

framework/include/systems/LinearSystem.h

@@ -101,6 +101,14 @@ class LinearSystem : public SolverSystem, public PerfGraphInterface
    */
   void computeGradients();
 
+  /**
+   * Return a reference to the new (temporary) gradient container vectors
+   */
+  std::vector<std::unique_ptr<NumericVector<Number>>> & newGradientContainer()
+  {
+    return _new_gradient;
+  }
+
 protected:
   /**
    * Compute the right hand side and system matrix for given tags
@@ -152,6 +160,12 @@ class LinearSystem : public SolverSystem, public PerfGraphInterface
   /// Base class reference to the linear implicit system in libmesh
   LinearImplicitSystem & _linear_implicit_system;
 
+  /// Vectors to store the new gradients during the computation. This is needed
+  /// because the old gradients might still be needed to determine boundary values
+  /// (for extrapolated boundary conditions). Once the computation is done, we
+  /// move the nev vectors to the original containers.
+  std::vector<std::unique_ptr<NumericVector<Number>>> _new_gradient;
+
 private:
   /// The current states of the solution (0 = current, 1 = old, etc)
   std::vector<NumericVector<Number> *> _solution_state;

framework/include/variables/MooseVariableField.h

@@ -392,10 +392,6 @@ class MooseVariableField : public MooseVariableFieldBase,
 
   /// A dummy ADReal variable
   mutable ADReal _ad_real_dummy = 0;
-
-  /// Boolean to check if this variable needs a raw cell-gradient storage vector from the
-  /// system.
-  bool _needs_cell_gradients;
 };
 
 #define usingMooseVariableFieldMembers                                                             \

framework/src/loops/ComputeLinearFVGreenGaussGradientFaceThread.C

@@ -17,7 +17,8 @@ ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFa
     _dim(_fe_problem.mesh().dimension()),
     _linear_system_number(linear_system_num),
     _linear_system(libMesh::cast_ref<libMesh::LinearImplicitSystem &>(
-        _fe_problem.getLinearSystem(_linear_system_number).system()))
+        _fe_problem.getLinearSystem(_linear_system_number).system())),
+    _new_gradient(_fe_problem.getLinearSystem(_linear_system_number).newGradientContainer())
 {
 }
 
@@ -26,7 +27,10 @@ ComputeLinearFVGreenGaussGradientFaceThread::ComputeLinearFVGreenGaussGradientFa
   : _fe_problem(x._fe_problem),
     _dim(x._dim),
     _linear_system_number(x._linear_system_number),
-    _linear_system(x._linear_system)
+    _linear_system(x._linear_system),
+    // This will be the vector we work on since the old gradient might still be needed
+    // to compute extrapolated boundary conditions for example.
+    _new_gradient(x._new_gradient)
 {
 }
 
@@ -43,13 +47,7 @@ ComputeLinearFVGreenGaussGradientFaceThread::operator()(const FaceInfoRange & ra
     mooseAssert(_current_var,
                 "This should be a linear FV variable, did we somehow add a nonlinear variable to "
                 "the linear system?");
-    auto & grad_container = _fe_problem.getLinearSystem(_linear_system_number).gradientContainer();
     if (_current_var->needsGradientVectorStorage())
-    {
-      _new_gradient.clear();
-      for (auto & vec : grad_container)
-        _new_gradient.push_back(vec->zero_clone());
-
       // Iterate over all the elements in the range
       for (const auto & face_info : range)
       {
@@ -62,10 +60,6 @@ ComputeLinearFVGreenGaussGradientFaceThread::operator()(const FaceInfoRange & ra
                  current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
           onBoundaryFace(*face_info);
       }
-
-      for (const auto i : index_range(grad_container))
-        grad_container[i] = std::move(_new_gradient[i]);
-    }
   }
 }
 

framework/src/loops/ComputeLinearFVGreenGaussGradientVolumeThread.C

@@ -36,20 +36,26 @@ ComputeLinearFVGreenGaussGradientVolumeThread::operator()(const ElemInfoRange &
   ParallelUniqueId puid;
   _tid = puid.id;
 
-  for (const auto & variable :
-       _fe_problem.getLinearSystem(_linear_system_number).getVariables(_tid))
+  auto & linear_system = _fe_problem.getLinearSystem(_linear_system_number);
+
+  // This will be the vector we work on since the old gradient might still be needed
+  // to compute extrapolated boundary conditions for example.
+  auto & grad_container = linear_system.newGradientContainer();
+
+  for (const auto & variable : linear_system.getVariables(_tid))
   {
     _current_var = dynamic_cast<MooseLinearVariableFV<Real> *>(variable);
     mooseAssert(_current_var,
                 "This should be a linear FV variable, did we somehow add a nonlinear variable to "
                 "the linear system?");
-    auto & grad_container = _fe_problem.getLinearSystem(_linear_system_number).gradientContainer();
     if (_current_var->needsGradientVectorStorage())
     {
       const auto rz_radial_coord = _fe_problem.mesh().getAxisymmetricRadialCoord();
       const auto state = Moose::currentState();
+
       // Iterate over all the elements in the range
       for (const auto & elem_info : range)
+      {
         if (_current_var->hasBlocks(elem_info->subdomain_id()))
         {
           const auto coord_type = _fe_problem.mesh().getCoordSystem(elem_info->subdomain_id());
@@ -58,8 +64,10 @@ ComputeLinearFVGreenGaussGradientVolumeThread::operator()(const ElemInfoRange &
           const auto volume = elem_info->volume() * elem_info->coordFactor();
 
           for (const auto dim_index : index_range(grad_container))
-            grad_container[dim_index]->set(dof_id_elem,
-                                           (*grad_container[dim_index])(dof_id_elem) / volume);
+          {
+            const auto normalized_value = (*grad_container[dim_index])(dof_id_elem) / volume;
+            grad_container[dim_index]->set(dof_id_elem, normalized_value);
+          }
 
           if (coord_type == Moose::CoordinateSystemType::COORD_RZ)
           {
@@ -68,6 +76,7 @@ ComputeLinearFVGreenGaussGradientVolumeThread::operator()(const ElemInfoRange &
             grad_container[rz_radial_coord]->add(dof_id_elem, radial_contrib);
           }
         }
+      }
     }
   }
 }

framework/src/systems/LinearSystem.C

@@ -140,37 +140,41 @@ LinearSystem::computeLinearSystemTags(const std::set<TagID> & vector_tags,
 void
 LinearSystem::computeGradients()
 {
-  TIME_SECTION("LinearVariableFV_Gradients", 3 /*, "Computing Linear FV variable gradients"*/);
-
-  using ElemInfoRange = StoredRange<MooseMesh::const_elem_info_iterator, const ElemInfo *>;
-  ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
-                                _fe_problem.mesh().ownedElemInfoEnd());
+  _new_gradient.clear();
+  for (auto & vec : _raw_grad_container)
+    _new_gradient.push_back(vec->zero_clone());
 
-  using FaceInfoRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
-  FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
-                                _fe_problem.mesh().ownedFaceInfoEnd());
+  TIME_SECTION("LinearVariableFV_Gradients", 3 /*, "Computing Linear FV variable gradients"*/);
 
   PARALLEL_TRY
   {
+    using FaceInfoRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
+    FaceInfoRange face_info_range(_fe_problem.mesh().ownedFaceInfoBegin(),
+                                  _fe_problem.mesh().ownedFaceInfoEnd());
+
     ComputeLinearFVGreenGaussGradientFaceThread gradient_face_thread(
         _fe_problem, _fe_problem.linearSysNum(name()));
     Threads::parallel_reduce(face_info_range, gradient_face_thread);
   }
   PARALLEL_CATCH;
 
-  // We must assemble here since we may have added face contributions to cells owned by neighboring
-  // processes, and we must perform all our face summations before performing division by the cell
-  // volumes in our volume thread
-  for (auto & grad_component : _raw_grad_container)
-    grad_component->close();
+  for (auto & vec : _new_gradient)
+    vec->close();
 
   PARALLEL_TRY
   {
+    using ElemInfoRange = StoredRange<MooseMesh::const_elem_info_iterator, const ElemInfo *>;
+    ElemInfoRange elem_info_range(_fe_problem.mesh().ownedElemInfoBegin(),
+                                  _fe_problem.mesh().ownedElemInfoEnd());
+
     ComputeLinearFVGreenGaussGradientVolumeThread gradient_volume_thread(
         _fe_problem, _fe_problem.linearSysNum(name()));
     Threads::parallel_reduce(elem_info_range, gradient_volume_thread);
   }
   PARALLEL_CATCH;
+
+  for (const auto i : index_range(_raw_grad_container))
+    _raw_grad_container[i] = std::move(_new_gradient[i]);
 }
 
 void