Merge pull request #1786 from TomFischer/FixPETScNormComputation

Fix PETSc norm computation

NumLib/DOF/DOFTableUtil.cpp

@@ -14,33 +14,101 @@ namespace NumLib
 {
 namespace
 {
-
-template <class CalculateNorm>
+template <typename CalculateNorm>
 double norm(GlobalVector const& x, unsigned const global_component,
-             LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh, CalculateNorm calculate_norm)
+            LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
+            CalculateNorm calculate_norm)
 {
-    // TODO that also includes ghost nodes.
+#ifdef USE_PETSC
+    x.setLocalAccessibleVector();
+#endif
+
     double res = 0.0;
-    MeshLib::MeshSubsets const& mss = dof_table.getMeshSubsets(global_component);
-    for (unsigned i=0; i<mss.size(); i++)
+    MeshLib::MeshSubsets const& mss =
+        dof_table.getMeshSubsets(global_component);
+
+    for (unsigned i = 0; i < mss.size(); i++)
     {
         MeshLib::MeshSubset const& ms = mss.getMeshSubset(i);
         if (ms.getMeshID() != mesh.getID())
             continue;
         for (MeshLib::Node const* node : ms.getNodes())
         {
-            auto const value =
-                getNodalValue(x, mesh, dof_table, node->getID(), global_component);
-
+            auto const value = getNonGhostNodalValue(
+                x, mesh, dof_table, node->getID(), global_component);
             res = calculate_norm(res, value);
         }
     }
+    return res;
+}
+
+double norm1(GlobalVector const& x, unsigned const global_component,
+             LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
+{
+    double res =
+        norm(x, global_component, dof_table, mesh,
+             [](double res, double value) { return res + std::abs(value); });
+
+#ifdef USE_PETSC
+    double global_result = 0.0;
+    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
+                  PETSC_COMM_WORLD);
+    res = global_result;
+#endif
+    return res;
+}
+
+double norm2(GlobalVector const& x, unsigned const global_component,
+             LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
+{
+    double res =
+        norm(x, global_component, dof_table, mesh,
+             [](double res, double value) { return res + value * value; });
+
+#ifdef USE_PETSC
+    double global_result = 0.0;
+    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_SUM,
+                  PETSC_COMM_WORLD);
+    res = global_result;
+#endif
+    return std::sqrt(res);
+}
 
-    // TODO for PETSc some global accumulation is necessary.
+double normInfinity(GlobalVector const& x, unsigned const global_component,
+                    LocalToGlobalIndexMap const& dof_table,
+                    MeshLib::Mesh const& mesh)
+{
+    double res = norm(x, global_component, dof_table, mesh,
+                      [](double res, double value) {
+                          return std::max(res, std::abs(value));
+                      });
+
+#ifdef USE_PETSC
+    double global_result = 0.0;
+    MPI_Allreduce(&res, &global_result, 1, MPI_DOUBLE, MPI_MAX,
+                  PETSC_COMM_WORLD);
+    res = global_result;
+#endif
     return res;
 }
+}  // anonymous namespace
+
+double getNonGhostNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
+                             NumLib::LocalToGlobalIndexMap const& dof_table,
+                             std::size_t const node_id,
+                             std::size_t const global_component_id)
+{
+    MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
+                              node_id};
 
-} // anonymous namespace
+    auto const index = dof_table.getGlobalIndex(l, global_component_id);
+    assert (index != NumLib::MeshComponentMap::nop);
+
+    if (index < 0)  // ghost node value
+        return 0.0;
+
+    return x.get(index);
+}
 
 double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
                      NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -50,8 +118,7 @@ double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
     MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
                               node_id};
 
-    auto const index = dof_table.getLocalIndex(
-        l, global_component_id, x.getRangeBegin(), x.getRangeEnd());
+    auto const index = dof_table.getGlobalIndex(l, global_component_id);
     assert (index != NumLib::MeshComponentMap::nop);
 
     return x.get(index);
@@ -95,28 +162,20 @@ NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(
     return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
 }
 
-
 double norm(GlobalVector const& x, unsigned const global_component,
             MathLib::VecNormType norm_type,
             LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
 {
-    switch(norm_type)
+    switch (norm_type)
     {
-    case MathLib::VecNormType::NORM1:
-        return norm(
-            x, global_component, dof_table, mesh,
-            [](double res, double value) { return res + std::abs(value); });
-    case MathLib::VecNormType::NORM2:
-        return std::sqrt(
-            norm(x, global_component, dof_table, mesh,
-                 [](double res, double value) { return res + value * value; }));
-    case MathLib::VecNormType::INFINITY_N:
-        return norm(x, global_component, dof_table, mesh,
-                    [](double res, double value) {
-                        return std::max(res, std::abs(value));
-                    });
-    default:
-        OGS_FATAL("An invalid norm type has been passed.");
+        case MathLib::VecNormType::NORM1:
+            return norm1(x, global_component, dof_table, mesh);
+        case MathLib::VecNormType::NORM2:
+            return norm2(x, global_component, dof_table, mesh);
+        case MathLib::VecNormType::INFINITY_N:
+            return normInfinity(x, global_component, dof_table, mesh);
+        default:
+            OGS_FATAL("An invalid norm type has been passed.");
     }
 }
 

NumLib/DOF/DOFTableUtil.h

@@ -14,6 +14,13 @@
 
 namespace NumLib
 {
+//! Returns the value for the given \c node_id and \c global_component_id from
+//! the vector \c x in case the node is not a ghost node. Else 0.0 is returned.
+double getNonGhostNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
+                             NumLib::LocalToGlobalIndexMap const& dof_table,
+                             std::size_t const node_id,
+                             std::size_t const global_component_id);
+
 //! Returns the value for the given \c node_id and \c global_component_id from
 //! the vector \c x.
 double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,

NumLib/DOF/LocalToGlobalIndexMap.h

@@ -129,6 +129,12 @@ public:
         return _mesh_component_map.getGlobalIndex(l, c);
     }
 
+    GlobalIndexType getGlobalIndex(MeshLib::Location const& l,
+                                   int const global_component_id) const
+    {
+        return _mesh_component_map.getGlobalIndex(l, global_component_id);
+    }
+
     /// Forwards the respective method from MeshComponentMap.
     std::vector<GlobalIndexType> getGlobalIndices(const MeshLib::Location &l) const
     {

NumLib/DOF/MeshComponentMap.cpp

@@ -77,7 +77,7 @@ MeshComponentMap::MeshComponentMap(
                     // arrange non ghost entries of a partition within
                     // a rank in the parallel computing.
                     OGS_FATAL("Global index in the system of equations"
-                              " can only be numbered by the oder type"
+                              " can only be numbered by the order type"
                               " of ComponentOrder::BY_LOCATION");
                 }
                 global_id = static_cast<GlobalIndexType>(

Tests/CMakeLists.txt

@@ -93,7 +93,7 @@ add_custom_target(tests-cleanup ${CMAKE_COMMAND} -E remove testrunner.xml)
 set_target_properties(tests-cleanup PROPERTIES FOLDER Testing)
 
 if(OGS_USE_PETSC)
-    set(TEST_FILTER_MPI --gtest_filter=-MPITest_Math.*)
+    set(TEST_FILTER_MPI --gtest_filter=-MPITest*)
     add_custom_target(tests
         mpirun -np 1 $<TARGET_FILE:testrunner> ${TESTRUNNER_ADDITIONAL_ARGUMENTS} ${TEST_FILTER_MPI}
         DEPENDS testrunner tests-cleanup

Tests/NumLib/TestComponentNorms.cpp

@@ -11,10 +11,17 @@
 #include <limits>
 #include <logog/include/logog.hpp>
 
+#ifdef USE_PETSC
+#include "BaseLib/BuildInfo.h"
+#endif
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/MatrixSpecifications.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
+#ifdef USE_PETSC
+#include "MeshLib/IO/readMeshFromFile.h"
+#else
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
+#endif
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "NumLib/NumericsConfig.h"
 #include "NumLib/DOF/DOFTableUtil.h"
@@ -34,8 +41,14 @@ struct DOFTableData
 {
     DOFTableData(unsigned const num_components,
                  NumLib::ComponentOrder const order)
+#ifdef USE_PETSC
+        : mesh(MeshLib::IO::readMeshFromFile(
+              BaseLib::BuildInfo::data_path +
+              "/EllipticPETSc/quad_20x10_GroundWaterFlow.")),
+#else
         : mesh(MeshLib::MeshGenerator::generateRegularHexMesh(
               mesh_length, mesh_elements_in_each_direction)),
+#endif
           mesh_subset_all_nodes(*mesh, &mesh->getNodes()),
           dof_table(createMeshSubsets(mesh_subset_all_nodes, num_components),
                     order)
@@ -46,12 +59,16 @@ struct DOFTableData
     MeshLib::MeshSubset const mesh_subset_all_nodes;
     NumLib::LocalToGlobalIndexMap const dof_table;
 
+#ifndef USE_PETSC
 private:
     static const double mesh_length;
     static const std::size_t mesh_elements_in_each_direction;
+#endif
 };
+#ifndef USE_PETSC
 const double DOFTableData::mesh_length = 1;
 const std::size_t DOFTableData::mesh_elements_in_each_direction = 5;
+#endif
 
 template <typename AccumulateCallback, typename AccumulateFinishCallback>
 void do_test(unsigned const num_components,
@@ -65,7 +82,11 @@ void do_test(unsigned const num_components,
      * vector.
      */
     using CO = NumLib::ComponentOrder;
+#ifdef USE_PETSC
+    for (auto order : {CO::BY_LOCATION})
+#else
     for (auto order : {CO::BY_COMPONENT, CO::BY_LOCATION})
+#endif
     {
         DOFTableData dtd(num_components, order);
 
@@ -97,7 +118,7 @@ void do_test(unsigned const num_components,
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormSingleComponent)
 #else
-TEST(NumLib, DISABLED_ComponentNormSingleComponent)
+TEST(MPITest_NumLib, ComponentNormSingleComponent)
 #endif
 {
     unsigned const num_components = 1;
@@ -116,7 +137,7 @@ TEST(NumLib, DISABLED_ComponentNormSingleComponent)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponent1)
 #else
-TEST(NumLib, DISABLED_ComponentNormMultiComponent1)
+TEST(MPITest_NumLib, ComponentNormMultiComponent1)
 #endif
 {
     unsigned const num_components = 3;
@@ -132,7 +153,7 @@ TEST(NumLib, DISABLED_ComponentNormMultiComponent1)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponent2)
 #else
-TEST(NumLib, DISABLED_ComponentNormMultiComponent2)
+TEST(MPITest_NumLib, ComponentNormMultiComponent2)
 #endif
 {
     unsigned const num_components = 3;
@@ -148,7 +169,7 @@ TEST(NumLib, DISABLED_ComponentNormMultiComponent2)
 #ifndef USE_PETSC
 TEST(NumLib, ComponentNormMultiComponentInfinity)
 #else
-TEST(NumLib, DISABLED_ComponentNormMultiComponentInfinity)
+TEST(MPITest_NumLib, ComponentNormMultiComponentInfinity)
 #endif
 {
     unsigned const num_components = 3;

Tests/VectorUtils.h

@@ -24,4 +24,7 @@ void fillVectorRandomly(Vector& x)
     for (Index i = 0; i < size; ++i) {
         MathLib::setVector(x, i, rnd(random_number_generator));
     }
+#ifdef USE_PETSC
+    finalizeVectorAssembly(x);
+#endif
 }