Initial cache implementation

examples/init/ellipt.dat.3d

@@ -4,7 +4,7 @@ ellipt->mesh:                   elliptMesh
 
 ellipt->solver->name:           bicgstab
 ellipt->solver->max iteration:  1000
-ellipt->solver->relative tolerance: 1.e-8
+ellipt->solver->relative tolerance: 1.e-9
 ellipt->solver->info:           10
 ellipt->solver->left precon:    diag
 ellipt->solver->right precon:   no

examples/init/stokes.dat.2d

@@ -9,7 +9,7 @@ stokes->mesh: stokesMesh
 stokes->solver->name: bicgstab_ell
 stokes->solver->ell: 3
 stokes->solver->max iteration: 1000
-stokes->solver->reduction: 1e-8
+stokes->solver->relative tolerance: 1e-8
 stokes->solver->restart: 50
 stokes->solver->info:	2
 

src/amdis/LocalBasisCache.hpp

+#pragma once
+
+#include <vector>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+
+#include <amdis/common/TupleUtility.hpp>
+#include <amdis/utility/ConcurrentCache.hpp>
+
+namespace AMDiS
+{
+  namespace Impl
+  {
+    struct HashableGeometryType
+    {
+      using type = std::tuple<unsigned int,unsigned int,bool>;
+      type value;
+
+      HashableGeometryType(Dune::GeometryType const& type)
+        : value{type.id(), type.dim(), type.isNone()}
+      {}
+
+      friend std::size_t hash_value(HashableGeometryType const& t)
+      {
+        std::size_t seed = 0;
+        Impl::HashTupleImpl<HashableGeometryType::type>::apply(seed, t.value);
+        return seed;
+      }
+
+      friend bool operator==(HashableGeometryType const& lhs, HashableGeometryType const& rhs)
+      {
+        return lhs.value == rhs.value;
+      }
+    };
+
+    /// Pair of GeometryType and quadrature order and size
+    template <class LocalBasisType, class Tag>
+    struct LocalBasisCacheKey
+    {
+      using type = std::tuple<HashableGeometryType,int,int>;
+      type value;
+
+      friend std::size_t hash_value(LocalBasisCacheKey const& t)
+      {
+        std::size_t seed = 0;
+        Impl::HashTupleImpl<LocalBasisCacheKey::type>::apply(seed, t.value);
+        return seed;
+      }
+
+      friend bool operator==(LocalBasisCacheKey const& lhs, LocalBasisCacheKey const& rhs)
+      {
+        return lhs.value == rhs.value;
+      }
+    };
+  }
+} // end namespace AMDiS
+
+DUNE_DEFINE_STD_HASH( , AMDiS::Impl::HashableGeometryType)
+DUNE_DEFINE_HASH(DUNE_HASH_TEMPLATE_ARGS(class LocalBasisType, class Tag),DUNE_HASH_TYPE(AMDiS::Impl::LocalBasisCacheKey<LocalBasisType,Tag>))
+
+namespace AMDiS
+{
+  template <class LocalBasisType>
+  class LocalBasisCache
+  {
+  public:
+    using Traits = typename LocalBasisType::Traits;
+
+    using RangeType = typename Traits::RangeType;
+    using RangeFieldType = typename Traits::RangeFieldType;
+    using JacobianType = typename Traits::JacobianType;
+
+    using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
+
+  private:
+    using ValuesKey = Impl::LocalBasisCacheKey<LocalBasisType, struct ValuesTag>;
+    using GradientsKey = Impl::LocalBasisCacheKey<LocalBasisType, struct GradientsTag>;
+    using Policy = tag::thread_local_policy;
+
+  public:
+    using ShapeValuesContainer = std::vector<std::vector<RangeType>>;
+    using ShapeGradientsContainer = std::vector<std::vector<JacobianType>>;
+
+    LocalBasisCache(LocalBasisType const& localBasis)
+      : localBasis_(&localBasis)
+    {}
+
+    template <class LocalContext, class QuadratureRule>
+    ShapeValuesContainer const& values(LocalContext const& context, QuadratureRule const& quad)
+    {
+      ValuesKey key{typename ValuesKey::type{context.type(),quad.order(),quad.size()}};
+      return ShapeValuesCache::get(key, [&](ShapeValuesContainer* data, ValuesKey const&)
+      {
+        data->resize(quad.size());
+        for (std::size_t iq = 0; iq < quad.size(); ++iq)
+          localBasis_->evaluateFunction(context.local(quad[iq].position()), (*data)[iq]);
+      });
+    }
+
+    template <class LocalContext, class QuadratureRule>
+    ShapeGradientsContainer const& gradients(LocalContext const& context, QuadratureRule const& quad)
+    {
+      GradientsKey key{typename GradientsKey::type{context.type(),quad.order(),quad.size()}};
+      return ShapeGradientsCache::get(key, [&](ShapeGradientsContainer* data, GradientsKey const&)
+      {
+        data->resize(quad.size());
+        for (std::size_t iq = 0; iq < quad.size(); ++iq)
+          localBasis_->evaluateJacobian(context.local(quad[iq].position()), (*data)[iq]);
+      });
+    }
+
+  private:
+    using ShapeValuesCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, Policy>;
+    using ShapeGradientsCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, Policy>;
+
+    LocalBasisType const* localBasis_;
+  };
+
+} // end namespace AMDiS

src/amdis/assembler/SecondOrderGradTestGradTrial.hpp

@@ -117,21 +117,22 @@ namespace AMDiS
 
       using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
       using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
-      using JacobianType = typename LocalBasisType::Traits::JacobianType;
 
-      for (auto const& qp : quad) {
+      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+      auto const& shapeGradientsCache = cache.gradients(context, quad);
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(qp.position());
+        decltype(auto) local = context.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = context.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
+        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
+                                                      * quad[iq].weight();
 
         // The gradients of the shape functions on the reference element
-        thread_local std::vector<JacobianType> shapeGradients;
-        localFE.localBasis().evaluateJacobian(local, shapeGradients);
+        auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
         thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients;

src/amdis/assembler/StokesOperator.hpp

@@ -53,25 +53,30 @@ namespace AMDiS
       using VelocityLocalBasisType = typename std::decay_t<decltype(velocityLocalFE)>::Traits::LocalBasisType;
       using PressureLocalBasisType = typename std::decay_t<decltype(pressureLocalFE)>::Traits::LocalBasisType;
 
+      LocalBasisCache<VelocityLocalBasisType> velocityCache(velocityLocalFE.localBasis());
+      LocalBasisCache<PressureLocalBasisType> pressureCache(pressureLocalFE.localBasis());
+
       using PressureRange = typename PressureLocalBasisType::Traits::RangeType;
       using RangeFieldType = typename VelocityLocalBasisType::Traits::RangeFieldType;
       using VelocityJacobian = typename VelocityLocalBasisType::Traits::JacobianType;
 
       auto const& quad = this->getQuadratureRule(context.type(), tree, tree);
-      for (auto const& qp : quad) {
+
+      auto const& shapeGradientsCache = velocityCache.gradients(context, quad);
+      auto const& pressureValuesCache = pressureCache.values(context, quad);
+
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(qp.position());
+        decltype(auto) local = context.local(quad[iq].position());
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = context.geometry().jacobianInverseTransposed(local);
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = context.integrationElement(qp.position()) * qp.weight();
+        const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
         const auto vel_factor = Super::coefficient(local) * factor;
 
-        // The gradients of the shape functions on the reference element
-        thread_local std::vector<VelocityJacobian> shapeGradients;
-        velocityLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
+        auto const& shapeGradients = shapeGradientsCache[iq];
 
         // Compute the shape function gradients on the real element
         thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients;
@@ -79,8 +84,7 @@ namespace AMDiS
         for (std::size_t i = 0; i < gradients.size(); ++i)
           jacobian.mv(shapeGradients[i][0], gradients[i]);
 
-        thread_local std::vector<PressureRange> pressureValues;
-        pressureLocalFE.localBasis().evaluateFunction(local, pressureValues);
+        auto const& pressureValues = pressureValuesCache[iq];
 
         // <viscosity * grad(u), grad(v)>
         for (std::size_t i = 0; i < velocityLocalFE.size(); ++i) {

src/amdis/assembler/ZeroOrderTest.hpp

@@ -3,6 +3,7 @@
 #include <type_traits>
 
 #include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalBasisCache.hpp>
 #include <amdis/common/ValueCategory.hpp>
 
 namespace AMDiS
@@ -44,19 +45,19 @@ namespace AMDiS
       auto const& localFE = node.finiteElement();
 
       using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
-      using RangeType = typename LocalBasisType::Traits::RangeType;
+      LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
 
       auto const& quad = this->getQuadratureRule(context.type(), node);
-      for (auto const& qp : quad) {
+      auto const& shapeValuesCache = cache.values(context,quad);
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
         // Position of the current quadrature point in the reference element
-        decltype(auto) local = context.local(qp.position());
+        decltype(auto) local = context.local(quad[iq].position());
 
         // The multiplicative factor in the integral transformation formula
-        const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight();
-
-        thread_local std::vector<RangeType> shapeValues;
-        localFE.localBasis().evaluateFunction(local, shapeValues);
+        const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position())
+                                                      * quad[iq].weight();
 
+        auto const& shapeValues = shapeValuesCache[iq];
         for (std::size_t i = 0; i < localFE.size(); ++i) {
           const auto local_i = node.localIndex(i);
           elementVector[local_i] += factor * shapeValues[i];