From ebf78a37af3a3b36492e156278fbcd8ca9866349 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Mon, 20 Aug 2018 13:58:03 +0200
Subject: [PATCH] Executable analyzing facedir/facemod variability for a grid

---
 bin/CMakeLists.txt                       |   4 +-
 bin/analyzegrid/CMakeLists.txt           |   2 +
 bin/analyzegrid/analyze_grid.cc          | 151 +++++++++++++++++++++++
 bin/analyzegrid/test_2d_structured.ini   |  26 ++++
 bin/analyzegrid/test_2d_unstructured.ini |  29 +++++
 bin/analyzegrid/test_3d_structured.ini   |  26 ++++
 bin/analyzegrid/test_3d_unstructured.ini |  30 +++++
 dune/perftool/sumfact/analyzegrid.hh     | 139 +++++++++++++++++++++
 python/dune/perftool/options.py          |   1 +
 python/dune/perftool/sumfact/switch.py   |  16 ++-
 10 files changed, 422 insertions(+), 2 deletions(-)
 create mode 100644 bin/analyzegrid/CMakeLists.txt
 create mode 100644 bin/analyzegrid/analyze_grid.cc
 create mode 100644 bin/analyzegrid/test_2d_structured.ini
 create mode 100644 bin/analyzegrid/test_2d_unstructured.ini
 create mode 100644 bin/analyzegrid/test_3d_structured.ini
 create mode 100644 bin/analyzegrid/test_3d_unstructured.ini
 create mode 100644 dune/perftool/sumfact/analyzegrid.hh

diff --git a/bin/CMakeLists.txt b/bin/CMakeLists.txt
index 5dd2d517..d334abb7 100644
--- a/bin/CMakeLists.txt
+++ b/bin/CMakeLists.txt
@@ -2,4 +2,6 @@ dune_python_install_package(PATH .)
 
 dune_symlink_to_source_files(FILES make_graph.sh
                                    knltimings.sh
-                             )
+                                   )
+
+add_subdirectory(analyzegrid)
diff --git a/bin/analyzegrid/CMakeLists.txt b/bin/analyzegrid/CMakeLists.txt
new file mode 100644
index 00000000..2e96704b
--- /dev/null
+++ b/bin/analyzegrid/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(analyze_grid analyze_grid.cc)
+dune_symlink_to_source_files(FILES test_2d_structured.ini test_3d_structured.ini test_2d_unstructured.ini test_3d_unstructured.ini)
diff --git a/bin/analyzegrid/analyze_grid.cc b/bin/analyzegrid/analyze_grid.cc
new file mode 100644
index 00000000..c4bf9eb6
--- /dev/null
+++ b/bin/analyzegrid/analyze_grid.cc
@@ -0,0 +1,151 @@
+#include "config.h"
+
+#include <algorithm>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include "dune/common/parametertreeparser.hh"
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include <dune/testtools/gridconstruction.hh>
+
+#include <dune/consistent-edge-orientation/createconsistentgrid.hh>
+
+#include <dune/perftool/sumfact/analyzegrid.hh>
+
+int main(int argc, char** argv){
+  try
+  {
+    if (argc != 2){
+      std::cout << "Need ini file as command line argument." << std::endl;
+      return 1;
+    }
+
+    // Initialize basic stuff...
+    Dune::MPIHelper& mpihelper = Dune::MPIHelper::instance(argc, argv);
+    using RangeType = double;
+    Dune::ParameterTree initree;
+    Dune::ParameterTreeParser::readINITree(argv[1], initree);
+
+    bool quadrilateral = false;
+    std::string elementType = initree.get("elementType", "");
+    if (elementType == "quadrilateral"){
+      quadrilateral = true;
+    }
+    quadrilateral = quadrilateral || initree.hasKey("cells");
+    std::cout << "quadrilateral: " << quadrilateral << std::endl;
+
+    bool unstructured = initree.get("formcompiler.grid_unstructured", false);
+    std::cout << "unstructured: " << unstructured << std::endl;
+
+    std::string filename = initree.get("analyze_grid_filename", "analyze_grid.txt");
+    std::cout << "writing into: " << filename << std::endl;
+
+    if (quadrilateral){
+      if (unstructured){
+        std::vector<int> tmp;
+        std::vector<int> dim_counter = initree.get("elements", tmp);
+        const int dim = dim_counter.size();
+        std::cout << "dim: " << dim << std::endl;
+
+        if (dim == 2){
+          std::cout << "Analyse 2d unstructured quadrilateral grid" << std::endl;
+
+          // Setup grid (view)...
+          using Grid = Dune::UGGrid<2>;
+          using GV = Grid::LeafGridView;
+          IniGridFactory<Grid> factory(initree);
+          std::shared_ptr<Grid> grid_nonconsistent = factory.getGrid();
+          std::shared_ptr<Grid> grid = createConsistentGrid(grid_nonconsistent);
+          GV gv = grid->leafGridView();
+
+          // Extract facemod/facedir intersection variation
+          using ES = Dune::PDELab::AllEntitySet<GV>;
+          ES es(gv);
+          analyze_grid(es, filename);
+        }
+        else if (dim == 3){
+          std::cout << "Analyse 3d unstructured quadrilateral grid" << std::endl;
+
+          // Setup grid (view)...
+          using Grid = Dune::UGGrid<3>;
+          using GV = Grid::LeafGridView;
+          IniGridFactory<Grid> factory(initree);
+          std::shared_ptr<Grid> grid_nonconsistent = factory.getGrid();
+          std::shared_ptr<Grid> grid = createConsistentGrid(grid_nonconsistent);
+          GV gv = grid->leafGridView();
+
+          // Extract facemod/facedir intersection variation
+          using ES = Dune::PDELab::AllEntitySet<GV>;
+          ES es(gv);
+          analyze_grid(es, filename);
+        }
+        else{
+          assert(false);
+        }
+      }
+      else {
+        std::vector<int> tmp;
+        std::vector<int> dim_counter = initree.get("cells", tmp);
+        const int dim = dim_counter.size();
+        std::cout << "dim: " << dim << std::endl;
+
+        if (dim == 2){
+          std::cout << "Analyse 2d structured quadrilateral grid" << std::endl;
+
+          // For structured grids we already know the variation beforehand. This
+          // is only implemented to cover all the cases.
+
+          // Setup grid (view)...
+          using Grid = Dune::YaspGrid<2>;
+          using GV = Grid::LeafGridView;
+          IniGridFactory<Grid> factory(initree);
+          std::shared_ptr<Grid> grid = factory.getGrid();
+          GV gv = grid->leafGridView();
+
+          // Extract facemod/facedir intersection variation
+          using ES = Dune::PDELab::AllEntitySet<GV>;
+          ES es(gv);
+          analyze_grid(es, filename);
+        }
+        if (dim == 3){
+          std::cout << "Analyse 3d structured quadrilateral grid" << std::endl;
+
+          // For structured grids we already know the variation beforehand. This
+          // is only implemented to cover all the cases.
+
+          // Setup grid (view)...
+          using Grid = Dune::YaspGrid<3>;
+          using GV = Grid::LeafGridView;
+          IniGridFactory<Grid> factory(initree);
+          std::shared_ptr<Grid> grid = factory.getGrid();
+          GV gv = grid->leafGridView();
+
+          // Extract facemod/facedir intersection variation
+          using ES = Dune::PDELab::AllEntitySet<GV>;
+          ES es(gv);
+          analyze_grid(es, filename);
+        }
+        else{
+          assert (false);
+        }
+      }
+    }
+    else{
+      // We can't do sum factorization on simplicial meshes
+      assert (false);
+    }
+
+    return 0;
+
+  }
+  catch (Dune::Exception& e)
+  {    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (std::exception& e)
+  {    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}
diff --git a/bin/analyzegrid/test_2d_structured.ini b/bin/analyzegrid/test_2d_structured.ini
new file mode 100644
index 00000000..29eca446
--- /dev/null
+++ b/bin/analyzegrid/test_2d_structured.ini
@@ -0,0 +1,26 @@
+__exec_suffix = deg2_symdiff_nonquadvec_nongradvec
+__name = /home/rene/phd/dune-les/build-optimized-clang60-python3/dune-perftool/test/sumfact/poisson/sumfact_poisson_dg_2d_deg2_symdiff_nonquadvec_nongradvec
+cells = 16 16
+deg_suffix = deg2
+diff_suffix = symdiff
+extension = 1. 1.
+gradvec_suffix = nongradvec
+quadvec_suffix = nonquadvec
+
+[formcompiler]
+compare_l2errorsquared = 5e-7
+
+[formcompiler.r]
+numerical_jacobian = 0
+sumfact = 1
+vectorization_quadloop = 0
+vectorization_strategy = none
+
+[formcompiler.ufl_variants]
+degree = 2
+
+[wrapper]
+
+[wrapper.vtkcompare]
+extension = vtu
+name = sumfact_poisson_dg_2d_deg2_symdiff_nonquadvec_nongradvec
diff --git a/bin/analyzegrid/test_2d_unstructured.ini b/bin/analyzegrid/test_2d_unstructured.ini
new file mode 100644
index 00000000..7e9d460d
--- /dev/null
+++ b/bin/analyzegrid/test_2d_unstructured.ini
@@ -0,0 +1,29 @@
+__exec_suffix = deg2_symdiff_nonquadvec_nongradvec
+__name = /home/rene/phd/dune-les/build-optimized-clang60-python3/dune-perftool/test/sumfact/poisson/sumfact_poisson_dg_2d_unstructured_deg2_symdiff_nonquadvec_nongradvec
+deg_suffix = deg2
+diff_suffix = symdiff
+elementType = quadrilateral
+elements = 16 16
+gradvec_suffix = nongradvec
+lowerleft = 0.0 0.0
+quadvec_suffix = nonquadvec
+upperright = 1.0 1.0
+
+[formcompiler]
+compare_l2errorsquared = 5e-5
+grid_unstructured = 1
+
+[formcompiler.r]
+numerical_jacobian = 0
+sumfact = 1
+vectorization_quadloop = 0
+vectorization_strategy = none
+
+[formcompiler.ufl_variants]
+degree = 2
+
+[wrapper]
+
+[wrapper.vtkcompare]
+extension = vtu
+name = sumfact_poisson_dg_2d_unstructured_deg2_symdiff_nonquadvec_nongradvec
diff --git a/bin/analyzegrid/test_3d_structured.ini b/bin/analyzegrid/test_3d_structured.ini
new file mode 100644
index 00000000..575f8877
--- /dev/null
+++ b/bin/analyzegrid/test_3d_structured.ini
@@ -0,0 +1,26 @@
+__exec_suffix = deg2_symdiff_nonquadvec_nongradvec
+__name = /home/rene/phd/dune-les/build-optimized-clang60-python3/dune-perftool/test/sumfact/poisson/sumfact_poisson_dg_3d_deg2_symdiff_nonquadvec_nongradvec
+cells = 8 8 8
+deg_suffix = deg2
+diff_suffix = symdiff
+extension = 1. 1. 1.
+gradvec_suffix = nongradvec
+quadvec_suffix = nonquadvec
+
+[formcompiler]
+compare_l2errorsquared = 5e-6
+
+[formcompiler.r]
+numerical_jacobian = 0
+sumfact = 1
+vectorization_quadloop = 0
+vectorization_strategy = none
+
+[formcompiler.ufl_variants]
+degree = 2
+
+[wrapper]
+
+[wrapper.vtkcompare]
+extension = vtu
+name = sumfact_poisson_dg_3d_deg2_symdiff_nonquadvec_nongradvec
diff --git a/bin/analyzegrid/test_3d_unstructured.ini b/bin/analyzegrid/test_3d_unstructured.ini
new file mode 100644
index 00000000..900efe54
--- /dev/null
+++ b/bin/analyzegrid/test_3d_unstructured.ini
@@ -0,0 +1,30 @@
+__exec_suffix = deg2_symdiff_nonquadvec_nongradvec
+__name = /home/rene/phd/dune-les/build-optimized-clang60-python3/dune-perftool/test/sumfact/poisson/sumfact_poisson_dg_3d_unstructured_deg2_symdiff_nonquadvec_nongradvec
+deg_suffix = deg2
+diff_suffix = symdiff
+elementType = quadrilateral
+elements = 8 8 8
+gradvec_suffix = nongradvec
+lowerleft = 0.0 0.0 0.0
+quadvec_suffix = nonquadvec
+upperright = 1.0 1.0 1.0
+
+[formcompiler]
+compare_l2errorsquared = 5e-6
+grid_unstructured = 1
+
+[formcompiler.r]
+numerical_jacobian = 0
+sumfact = 1
+sumfact_regular_jacobians = 1
+vectorization_quadloop = 0
+vectorization_strategy = none
+
+[formcompiler.ufl_variants]
+degree = 2
+
+[wrapper]
+
+[wrapper.vtkcompare]
+extension = vtu
+name = sumfact_poisson_dg_3d_unstructured_deg2_symdiff_nonquadvec_nongradvec
diff --git a/dune/perftool/sumfact/analyzegrid.hh b/dune/perftool/sumfact/analyzegrid.hh
new file mode 100644
index 00000000..c6466161
--- /dev/null
+++ b/dune/perftool/sumfact/analyzegrid.hh
@@ -0,0 +1,139 @@
+#ifndef DUNE_PERFTOOL_SUMFACT_ANALYZEGRID_HH
+#define DUNE_PERFTOOL_SUMFACT_ANALYZEGRID_HH
+
+
+#include <dune/pdelab/common/intersectiontype.hh>
+
+
+struct GridInfo{
+  int dim;
+  int number_of_elements;
+  int number_of_skeleton_faces;
+  int number_of_boundary_faces;
+  std::set<std::size_t> skeleton_variants;
+  std::set<std::size_t> boundary_variants;
+
+  void export_grid_info(std::string);
+
+  friend std::ostream & operator<<(std::ostream &os, const GridInfo& grid_info);
+};
+
+void GridInfo::export_grid_info(std::string filename){
+  std::ofstream output_file;
+  output_file.open(filename);
+  for (std::size_t var : skeleton_variants){
+    std::size_t index_s = var / (2 * dim);
+    std::size_t facedir_s = index_s / 2;
+    std::size_t facemod_s = index_s % 2;
+
+    std::size_t index_n = var % (2 * dim);
+    std::size_t facedir_n = index_n / 2;
+    std::size_t facemod_n = index_n % 2;
+
+    output_file << "skeleton" << " "
+                << facedir_s << " "
+                << facemod_s << " "
+                << facedir_n << " "
+                << facemod_n << std::endl;
+  }
+  for (std::size_t var : boundary_variants){
+    std::size_t index_s = var;
+    std::size_t facedir_s = index_s / 2;
+    std::size_t facemod_s = index_s % 2;
+
+    output_file << "boundary" << " "
+                << facedir_s << " "
+                << facemod_s << " "
+                << -1 << " "
+                << -1 << std::endl;
+  }
+  output_file.close();
+}
+
+std::ostream & operator<<(std::ostream &os, const GridInfo& grid_info){
+  return os << "== Print GridInfo" << std::endl
+            << "Dimension: " << grid_info.dim << std::endl
+            << "Number of elements: " << grid_info.number_of_elements << std::endl
+            << std::endl
+            << "Number of skeleton faces: " << grid_info.number_of_skeleton_faces << std::endl
+            << "Number of skeleton variants: " << grid_info.skeleton_variants.size() << std::endl
+            << std::endl
+            << "Number of boundary faces: " << grid_info.number_of_boundary_faces << std::endl
+            << "Number of boundary variants: " << grid_info.boundary_variants.size() << std::endl;
+}
+
+
+template <typename ES>
+void analyze_grid(const ES& es, std::string filename){
+  const int dim = ES::dimension;
+
+  GridInfo grid_info;
+  grid_info.dim = dim;
+  grid_info.number_of_elements = 0;
+  grid_info.number_of_skeleton_faces = 0;
+  grid_info.number_of_boundary_faces = 0;
+
+
+  auto& index_set = es.indexSet();
+  for (const auto& element : elements(es)){
+    grid_info.number_of_elements += 1;
+
+    auto ids = index_set.index(element);
+
+    for (const auto& intersection : intersections(es, element)){
+      auto intersection_data = Dune::PDELab::classifyIntersection(es, intersection);
+      auto intersection_type = std::get<0>(intersection_data);
+      auto& outside_element = std::get<1>(intersection_data);
+
+      switch(intersection_type){
+      case Dune::PDELab::IntersectionType::skeleton :
+      {
+        auto idn = index_set.index(outside_element);
+        bool first_visit = ids > idn;
+        if (first_visit){
+          grid_info.number_of_skeleton_faces += 1;
+
+          size_t facedir_s = intersection.indexInInside() / 2 ;
+          size_t facemod_s = intersection.indexInInside() % 2 ;
+          size_t facedir_n = intersection.indexInOutside() / 2 ;
+          size_t facemod_n = intersection.indexInOutside() % 2 ;
+
+          // std::cout << facedir_s << " "
+          //           << facemod_s << " "
+          //           << facedir_n << " "
+          //           << facemod_n << std::endl;
+
+          size_t variant = intersection.indexInOutside() + 2 * dim * intersection.indexInInside();
+          grid_info.skeleton_variants.insert(variant);
+        }
+
+        break;
+      }
+      case Dune::PDELab::IntersectionType::periodic :
+      {
+        assert (false);
+        break;
+      }
+      case Dune::PDELab::IntersectionType::boundary :
+      {
+        grid_info.number_of_boundary_faces += 1;
+
+        size_t variant = intersection.indexInInside();
+        grid_info.boundary_variants.insert(variant);
+
+        break;
+      }
+      case Dune::PDELab::IntersectionType::processor :
+      {
+        assert (false);
+        break;
+      }
+      }
+    }
+  }
+
+  grid_info.export_grid_info(filename);
+  std::cout << grid_info << std::endl;
+}
+
+#endif
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index e73ff9ef..b909e5b1 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -80,6 +80,7 @@ class PerftoolFormOptionsArray(ImmutableRecord):
     fastdg = PerftoolOption(default=False, helpstr="Use FastDGGridOperator from PDELab.")
     sumfact = PerftoolOption(default=False, helpstr="Use sumfactorization")
     sumfact_regular_jacobians = PerftoolOption(default=False, helpstr="Generate non sum-factorized jacobians (only useful if sumfact is set)")
+    sumfact_grid_info = PerftoolOption(default=None, helpstr="Path to file with information about facedir and facemod variations. This can be used to limit the generation of skeleton kernels.")
     vectorization_quadloop = PerftoolOption(default=False, helpstr="whether to generate code with explicit vectorization")
     vectorization_strategy = PerftoolOption(default="none", helpstr="The identifier of the vectorization cost model. Possible values: none|explicit|model")
     vectorization_not_fully_vectorized_error = PerftoolOption(default=False, helpstr="throw an error if nonquadloop vectorization did not fully vectorize")
diff --git a/python/dune/perftool/sumfact/switch.py b/python/dune/perftool/sumfact/switch.py
index aee3c5a1..5d9e971f 100644
--- a/python/dune/perftool/sumfact/switch.py
+++ b/python/dune/perftool/sumfact/switch.py
@@ -1,5 +1,7 @@
 """ boundary and skeleton integrals come in variants in sum factorization - implement the switch! """
 
+import csv
+
 from dune.perftool.generation import (backend,
                                       get_global_context_value,
                                       global_context,
@@ -54,7 +56,17 @@ def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=No
 def decide_if_kernel_is_necessary(facedir_s, facemod_s, facedir_n, facemod_n):
     # If we are not using YaspGrid, all variants need to be realized
     if not get_form_option("diagonal_transformation_matrix"):
-        return True
+        # Reduce the variability according to grid info file
+        if get_form_option("sumfact_grid_info") is not None:
+            filename = get_form_option("sumfact_grid_info")
+            with open(filename) as csv_file:
+                csv_reader = csv.reader(csv_file, delimiter=" ")
+                for row in csv_reader:
+                    if (facedir_s == int(row[1])) and (facemod_s == int(row[2])) and (facedir_n == int(row[3])) and (facemod_n == int(row[4])):
+                        return True
+                return False
+        else:
+            return True
 
     # The PDELab machineries visit-once policy combined with Yasp avoids any visits
     # with facemod_s being True
@@ -98,6 +110,7 @@ def generate_exterior_facet_switch():
                                                                               ),
                                                               args))
 
+    block.append("    default: assert(false);")
     block.append("  }")
     block.append("}")
 
@@ -130,6 +143,7 @@ def generate_interior_facet_switch():
                                                                                           ),
                                                                           args))
 
+    block.append("    default: assert(false);")
     block.append("  }")
     block.append("}")
 
-- 
GitLab