Merge pull request #2082 from endJunction/EnableExtrapolationInPetsc

Enable extrapolation in PETSc

MathLib/LinAlg/LinAlg.cpp

@@ -308,6 +308,8 @@ void finalizeAssembly(EigenMatrix& x)
     x.getRawMatrix().makeCompressed();
 }
 
+void finalizeAssembly(EigenVector& x) {}
+
 } // namespace LinAlg
 
 } // namespace MathLib

MathLib/LinAlg/LinAlg.h

@@ -98,11 +98,8 @@ double norm(MatrixOrVector const& x, MathLib::VecNormType type)
     }
 }
 
-template<typename Matrix>
-void finalizeAssembly(Matrix& /*A*/)
-{
-    // By default do nothing.
-}
+template <typename Matrix>
+void finalizeAssembly(Matrix& /*A*/);
 
 // Matrix and Vector
 
@@ -260,6 +257,7 @@ void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
                 EigenVector const& v2, EigenVector& v3);
 
 void finalizeAssembly(EigenMatrix& A);
+void finalizeAssembly(EigenVector& A);
 
 } // namespace LinAlg
 

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp

@@ -47,12 +47,39 @@ void LocalLinearLeastSquaresExtrapolator::extrapolate(
 {
     auto const num_nodal_dof_result =
         _dof_table_single_component.dofSizeWithoutGhosts() * num_components;
+
+    std::vector<GlobalIndexType> ghost_indices;
+    {  // Create num_components times version of ghost_indices arranged by
+       // location. For example for 3 components and ghost_indices {5,6,10} we
+       // compute {15, 16, 17,  18, 19, 20,  30, 31, 32}.
+        auto const& single_component_ghost_indices =
+            _dof_table_single_component.getGhostIndices();
+        auto const single_component_ghost_indices_size =
+            single_component_ghost_indices.size();
+        ghost_indices.reserve(single_component_ghost_indices_size *
+                              num_components);
+        for (unsigned i = 0; i < single_component_ghost_indices_size; ++i)
+        {
+            for (unsigned c = 0; c < num_components; ++c)
+            {
+                ghost_indices.push_back(
+                    single_component_ghost_indices[i] * num_components + c);
+            }
+        }
+    }
+
     if (!_nodal_values ||
-        static_cast<std::size_t>(_nodal_values->size()) != num_nodal_dof_result)
+#ifdef USE_PETSC
+        static_cast<std::size_t>(_nodal_values->getLocalSize() +
+                                 _nodal_values->getGhostSize())
+#else
+        _nodal_values->size()
+#endif
+            != num_nodal_dof_result)
     {
         _nodal_values = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
-            MathLib::MatrixSpecifications{
-                num_nodal_dof_result, num_nodal_dof_result, nullptr, nullptr});
+            {num_nodal_dof_result, num_nodal_dof_result, &ghost_indices,
+             nullptr});
     }
     _nodal_values->setZero();
 
@@ -68,6 +95,7 @@ void LocalLinearLeastSquaresExtrapolator::extrapolate(
         extrapolateElement(i, num_components, extrapolatables, t,
                            current_solution, dof_table, *counts);
     }
+    MathLib::LinAlg::finalizeAssembly(*_nodal_values);
 
     MathLib::LinAlg::componentwiseDivide(*_nodal_values, *_nodal_values,
                                          *counts);
@@ -104,6 +132,7 @@ void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
         calculateResidualElement(i, num_components, extrapolatables, t,
                                  current_solution, dof_table);
     }
+    MathLib::LinAlg::finalizeAssembly(*_residuals);
 }
 
 void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
@@ -272,6 +301,8 @@ void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
     auto const int_pt_vals_mat =
         MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
 
+    MathLib::LinAlg::setLocalAccessibleVector(
+        *_nodal_values);  // For access in the for-loop.
     for (unsigned comp = 0; comp < num_components; ++comp)
     {
         // filter nodal values of the current element
@@ -279,7 +310,7 @@ void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
         {
             // TODO PETSc negative indices?
             auto const idx = num_components * global_indices[i] + comp;
-            nodal_vals_element[i] = (*_nodal_values)[idx];
+            nodal_vals_element[i] = _nodal_values->get(idx);
         }
 
         double const residual = (interpolation_matrix * nodal_vals_element -

ProcessLib/Output/CachedSecondaryVariable.cpp

@@ -68,6 +68,9 @@ GlobalVector const& CachedSecondaryVariable::evalFieldNoArgs() const
         1, *_extrapolatables, _t, *_current_solution, *_dof_table);
     auto const& nodal_values = _extrapolator.getNodalValues();
     MathLib::LinAlg::copy(nodal_values, _cached_nodal_values);
+
+    MathLib::LinAlg::setLocalAccessibleVector(
+        _cached_nodal_values);  // For access in the getValue()
     _needs_recomputation = false;
     return nodal_values;
 }

ProcessLib/Output/ProcessOutput.cpp

@@ -10,6 +10,7 @@
 #include "ProcessOutput.h"
 
 #include "BaseLib/BuildInfo.h"
+#include "MathLib/LinAlg/LinAlg.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
@@ -27,7 +28,6 @@ static void addOgsVersion(MeshLib::Mesh& mesh)
                              BaseLib::BuildInfo::ogs_version.end());
 }
 
-#ifndef USE_PETSC  // Not used in PETSc case
 static void addSecondaryVariableNodes(
     double const t,
     GlobalVector const& x,
@@ -55,21 +55,22 @@ static void addSecondaryVariableNodes(
     std::unique_ptr<GlobalVector> result_cache;
     auto const& nodal_values =
         var.fcts.eval_field(t, x, dof_table, result_cache);
-    if (nodal_values_mesh.size() !=
-        static_cast<std::size_t>(nodal_values.size()))
+#ifdef USE_PETSC
+    std::size_t const global_vector_size =
+        nodal_values.getLocalSize() + nodal_values.getGhostSize();
+#else
+    std::size_t const global_vector_size = nodal_values.size();
+#endif
+    if (nodal_values_mesh.size() != global_vector_size)
     {
         OGS_FATAL(
             "Secondary variable `%s' did not evaluate to the right "
             "number of components. Expected: %d, actual: %d.",
-            var.name.c_str(), nodal_values_mesh.size(), nodal_values.size());
+            var.name.c_str(), nodal_values_mesh.size(), global_vector_size);
     }
 
     // Copy result
-    for (GlobalIndexType i = 0; i < nodal_values.size(); ++i)
-    {
-        assert(!std::isnan(nodal_values[i]));
-        nodal_values_mesh[i] = nodal_values[i];
-    }
+    nodal_values.copyValues(nodal_values_mesh);
 }
 
 static void addSecondaryVariableResiduals(
@@ -103,22 +104,23 @@ static void addSecondaryVariableResiduals(
     std::unique_ptr<GlobalVector> result_cache;
     auto const& residuals =
         var.fcts.eval_residuals(t, x, dof_table, result_cache);
-    if (residuals_mesh.size() != static_cast<std::size_t>(residuals.size()))
+#ifdef USE_PETSC
+    std::size_t const global_vector_size =
+        residuals.getLocalSize() + residuals.getGhostSize();
+#else
+    std::size_t const global_vector_size = residuals.size();
+#endif
+    if (residuals_mesh.size() != global_vector_size)
     {
         OGS_FATAL(
             "The residual of secondary variable `%s' did not evaluate to the "
             "right number of components. Expected: %d, actual: %d.",
-            var.name.c_str(), residuals_mesh.size(), residuals.size());
+            var.name.c_str(), residuals_mesh.size(), global_vector_size);
     }
 
     // Copy result
-    for (GlobalIndexType i = 0; i < residuals.size(); ++i)
-    {
-        assert(!std::isnan(residuals[i]));
-        residuals_mesh[i] = residuals[i];
-    }
+    residuals.copyValues(residuals_mesh);
 }
-#endif  // USE_PETSC
 
 namespace ProcessLib
 {
@@ -212,8 +214,6 @@ void processOutputData(
         }
     }
 
-#ifndef USE_PETSC
-
     // Secondary variables output
     for (auto const& external_variable_name : output_variables)
     {
@@ -234,11 +234,6 @@ void processOutputData(
                 external_variable_name, mesh);
         }
     }
-#else
-    (void)mesh;
-    (void)secondary_variables;
-    (void)t;
-#endif  // USE_PETSC
 
     addIntegrationPointWriter(mesh, integration_point_writer);
 }