From 0528733ecfc75c34d6ab914b8d256bfe960bcb3f Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 3 May 2022 15:47:21 +0200
Subject: [PATCH 1/2] [test][1p] Fix missing overload of gridView.grid() for
 mock grid

---
 test/porousmediumflow/1p/incompressible/properties.hh | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/test/porousmediumflow/1p/incompressible/properties.hh b/test/porousmediumflow/1p/incompressible/properties.hh
index 5c59e2dfa1..d68b73ba04 100644
--- a/test/porousmediumflow/1p/incompressible/properties.hh
+++ b/test/porousmediumflow/1p/incompressible/properties.hh
@@ -71,6 +71,9 @@ public:
 
     struct Traits : public ParentType::Traits
     { using Grid = NoCommunicateGrid<dim>; };
+
+    const typename Traits::Grid& grid() const
+    { return static_cast<const typename Traits::Grid&>(ParentType::grid()); }
 };
 
 template<int dim>
-- 
GitLab


From e23967f7575d26ccdaef40978730a8e5c82e880e Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 3 May 2022 14:55:35 +0200
Subject: [PATCH 2/2] [bugfix][gridgeometry] Always construct bboxtree and
 element map

A thread-safe implementation with double-checked lock needs atomic load/store
operations to update the pointer to be thread-safe. Alternatively, a lock
has to be acquired every time when the function is called. Both implementations
impair the scalability of the multihreaded assembly (in the case of tpfa, or
any simulation that uses element maps or bounding box tree on an element level).

This changes the beahviour to always build the tree and the element map. This
usually occurs little runtime overhead in cases where the features are unused.
It is the most performant version as soon as one of the features is used.
There is however some memory overhead in the order of one or more ints per
element (seeds) and a bounding box hierarchy.
---
 dumux/discretization/basegridgeometry.hh | 49 ++++++++----------------
 1 file changed, 15 insertions(+), 34 deletions(-)

diff --git a/dumux/discretization/basegridgeometry.hh b/dumux/discretization/basegridgeometry.hh
index c14853150c..ae992f926e 100644
--- a/dumux/discretization/basegridgeometry.hh
+++ b/dumux/discretization/basegridgeometry.hh
@@ -26,7 +26,6 @@
 
 #include <utility>
 #include <type_traits>
-#include <mutex>
 
 #include <dune/grid/common/mcmgmapper.hh>
 
@@ -149,34 +148,13 @@ public:
      * \brief Returns the bounding box tree of the grid
      */
     const BoundingBoxTree& boundingBoxTree() const
-    {
-        if (!boundingBoxTree_)
-        {
-            std::scoped_lock{ boundingBoxTreeMutex_ };
-            if (!boundingBoxTree_)
-                elementMap(); // make sure the element map is built
-                boundingBoxTree_ = std::make_unique<BoundingBoxTree>(
-                    std::make_shared<ElementSet>(gridView_, elementMapper(), elementMap_)
-                );
-        }
-
-        return *boundingBoxTree_;
-    }
+    { return *boundingBoxTree_; }
 
     /*!
      * \brief Returns the element index to element map
      */
     const ElementMap& elementMap() const
-    {
-        if (!elementMap_)
-        {
-            std::scoped_lock{ elementMapMutex_ };
-            if (!elementMap_)
-                elementMap_ = std::make_shared<ElementMap>(gridView_.grid(), elementMapper_);
-        }
-
-        return *elementMap_;
-    }
+    { return *elementMap_; }
 
     /*!
      * \brief Get an element from a global element index
@@ -258,7 +236,7 @@ private:
 
     void update_()
     {
-        //! Update the mappers
+        // Update the mappers
         if constexpr (Deprecated::hasUpdateGridView<ElementMapper, GridView>())
             elementMapper_.update(gridView_);
         else
@@ -269,12 +247,17 @@ private:
         else
             Deprecated::update(vertexMapper_);
 
-        //! Compute the bounding box of the entire domain, for e.g. setting boundary conditions
+        // Compute the bounding box of the entire domain, for e.g. setting boundary conditions
         computeGlobalBoundingBox_();
 
-        //! reset bounding box tree and the element map until requested the next time
-        boundingBoxTree_.release();
-        elementMap_.reset();
+        // update element map and bounding box tree
+        // always building these comes at a memory overhead but improved
+        // performance and thread-safe element level access (e.g. during assembly)
+        // for all simulation that use these features
+        elementMap_ = std::make_shared<ElementMap>(gridView_.grid(), elementMapper_);
+        boundingBoxTree_ = std::make_unique<BoundingBoxTree>(
+            std::make_shared<ElementSet>(gridView_, elementMapper(), elementMap_)
+        );
     }
 
     //! the process grid view
@@ -285,12 +268,10 @@ private:
     VertexMapper vertexMapper_;
 
     //! the bounding box tree of the grid view for efficient element intersections
-    mutable std::unique_ptr<BoundingBoxTree> boundingBoxTree_;
-    mutable std::mutex boundingBoxTreeMutex_;
+    std::unique_ptr<const BoundingBoxTree> boundingBoxTree_;
 
-    //! a map from element index to elements (needed in the bounding box tree and for assembling cell-centered discretization)
-    mutable std::shared_ptr<ElementMap> elementMap_;
-    mutable std::mutex elementMapMutex_;
+    //! a map from element index to elements
+    std::shared_ptr<const ElementMap> elementMap_;
 
     //! the bounding box of the whole domain
     GlobalCoordinate bBoxMin_;
-- 
GitLab