Commit b01c33f7 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

caching in local operators and dofvectors

parent 7d87fdfe
...@@ -20,4 +20,4 @@ foreach(project ${projects3d}) ...@@ -20,4 +20,4 @@ foreach(project ${projects3d})
target_link_libraries(${project}.3d amdis) target_link_libraries(${project}.3d amdis)
target_compile_definitions(${project}.3d PRIVATE AMDIS_DIM=3 AMDIS_DOW=3) target_compile_definitions(${project}.3d PRIVATE AMDIS_DIM=3 AMDIS_DOW=3)
add_dependencies(examples ${project}.3d) add_dependencies(examples ${project}.3d)
endforeach() endforeach()
\ No newline at end of file
...@@ -32,7 +32,7 @@ int main(int argc, char** argv) ...@@ -32,7 +32,7 @@ int main(int argc, char** argv)
auto _p = Dune::Indices::_1; auto _p = Dune::Indices::_1;
auto opStokes = makeOperator(tag::stokes{}, viscosity); auto opStokes = makeOperator(tag::stokes{}, viscosity);
prob.addMatrixOperator(opStokes, treepath(), treepath()); prob.addMatrixOperator(opStokes);
auto opZero = makeOperator(tag::test_trial{}, 0.0); auto opZero = makeOperator(tag::test_trial{}, 0.0);
prob.addMatrixOperator(opZero, _p, _p); prob.addMatrixOperator(opZero, _p, _p);
......
...@@ -2,117 +2,163 @@ ...@@ -2,117 +2,163 @@
#include <vector> #include <vector>
#include <dune/common/hash.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh> #include <dune/geometry/type.hh>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/utility/ConcurrentCache.hpp> #include <amdis/utility/ConcurrentCache.hpp>
namespace AMDiS namespace AMDiS
{ {
namespace Impl template <class LocalBasisType>
class LocalBasisCache
{ {
struct HashableGeometryType public:
{ using Traits = typename LocalBasisType::Traits;
using type = std::tuple<unsigned int,unsigned int,bool>;
type value;
HashableGeometryType(Dune::GeometryType const& type) using DomainType = typename Traits::DomainType;
: value{type.id(), type.dim(), type.isNone()} using RangeType = typename Traits::RangeType;
{} using RangeFieldType = typename Traits::RangeFieldType;
using JacobianType = typename Traits::JacobianType;
friend std::size_t hash_value(HashableGeometryType const& t) using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
{ using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
std::size_t seed = 0;
Impl::HashTupleImpl<HashableGeometryType::type>::apply(seed, t.value);
return seed;
}
friend bool operator==(HashableGeometryType const& lhs, HashableGeometryType const& rhs) private:
{
return lhs.value == rhs.value;
}
};
/// Pair of GeometryType and quadrature order and size /// Pair of GeometryType and quadrature order and size
template <class LocalBasisType, class Tag> template <class Tag>
struct LocalBasisCacheKey struct QuadKey
{ {
using type = std::tuple<HashableGeometryType,int,int>; unsigned int id; // topologyId
type value; int p; // quadrature order
std::size_t size; // number of quadrature points
friend std::size_t hash_value(LocalBasisCacheKey const& t) struct hasher
{ {
std::size_t seed = 0; std::size_t operator()(QuadKey const& t) const
Impl::HashTupleImpl<LocalBasisCacheKey::type>::apply(seed, t.value); {
return seed; std::size_t seed = 0;
Dune::hash_combine(seed, t.id);
Dune::hash_combine(seed, t.p);
Dune::hash_combine(seed, t.size);
return seed;
}
};
friend bool operator==(QuadKey const& lhs, QuadKey const& rhs)
{
return std::tie(lhs.id,lhs.p,lhs.size) == std::tie(rhs.id,rhs.p,rhs.size);
} }
};
/// Pair of GeometryType and local coordinates
template <class Tag>
struct CoordsKey
{
unsigned int id; // topologyId
DomainType local; // local coordinate
friend bool operator==(LocalBasisCacheKey const& lhs, LocalBasisCacheKey const& rhs) struct hasher
{
std::size_t operator()(CoordsKey const& t) const
{
std::size_t seed = 0;
Dune::hash_combine(seed, t.id);
Dune::hash_range(seed, t.local.begin(), t.local.end());
return seed;
}
};
friend bool operator==(CoordsKey const& lhs, CoordsKey const& rhs)
{ {
return lhs.value == rhs.value; return std::tie(lhs.id,lhs.local) == std::tie(rhs.id,rhs.local);
} }
}; };
}
} // end namespace AMDiS
DUNE_DEFINE_STD_HASH( , AMDiS::Impl::HashableGeometryType) using LocalValuesKey = CoordsKey<struct ValuesTag>;
DUNE_DEFINE_HASH(DUNE_HASH_TEMPLATE_ARGS(class LocalBasisType, class Tag),DUNE_HASH_TYPE(AMDiS::Impl::LocalBasisCacheKey<LocalBasisType,Tag>)) using LocalGradientsKey = CoordsKey<struct GradientsTag>;
namespace AMDiS using ValuesKey = QuadKey<struct ValuesTag>;
{ using GradientsKey = QuadKey<struct GradientsTag>;
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>; public:
using ShapeValues = std::vector<RangeType>;
using ShapeGradients = std::vector<JacobianType>;
private: using ShapeValuesContainer = std::vector<ShapeValues>;
using ValuesKey = Impl::LocalBasisCacheKey<LocalBasisType, struct ValuesTag>; using ShapeGradientsContainer = std::vector<ShapeGradients>;
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) LocalBasisCache(LocalBasisType const& localBasis)
: localBasis_(&localBasis) : localBasis_(&localBasis)
{} {}
template <class LocalContext, class QuadratureRule> // evaluate local basis-function at all quadrature points
ShapeValuesContainer const& values(LocalContext const& context, QuadratureRule const& quad) template <class LocalContext>
ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, QuadratureRule const& quad) const
{ {
ValuesKey key{typename ValuesKey::type{context.type(),quad.order(),quad.size()}}; ValuesKey key{context.type().id(),quad.order(),quad.size()};
return ShapeValuesCache::get(key, [&](ShapeValuesContainer* data, ValuesKey const&) return ShapeValuesContainerCache::get(key, [&](ValuesKey const&)
{ {
data->resize(quad.size()); ShapeValuesContainer data(quad.size());
for (std::size_t iq = 0; iq < quad.size(); ++iq) for (std::size_t iq = 0; iq < quad.size(); ++iq)
localBasis_->evaluateFunction(context.local(quad[iq].position()), (*data)[iq]); localBasis_->evaluateFunction(context.local(quad[iq].position()), data[iq]);
return data;
}); });
} }
template <class LocalContext, class QuadratureRule> // evaluate local basis-gradients at all quadrature points
ShapeGradientsContainer const& gradients(LocalContext const& context, QuadratureRule const& quad) template <class LocalContext>
ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, QuadratureRule const& quad) const
{ {
GradientsKey key{typename GradientsKey::type{context.type(),quad.order(),quad.size()}}; GradientsKey key{context.type().id(),quad.order(),quad.size()};
return ShapeGradientsCache::get(key, [&](ShapeGradientsContainer* data, GradientsKey const&) return ShapeGradientsContainerCache::get(key, [&](GradientsKey const&)
{ {
data->resize(quad.size()); ShapeGradientsContainer data(quad.size());
for (std::size_t iq = 0; iq < quad.size(); ++iq) for (std::size_t iq = 0; iq < quad.size(); ++iq)
localBasis_->evaluateJacobian(context.local(quad[iq].position()), (*data)[iq]); localBasis_->evaluateJacobian(context.local(quad[iq].position()), data[iq]);
return data;
});
}
// evaluate local basis-function at point
template <class Element>
ShapeValues const& evaluateFunction(Element const& element, DomainType const& local) const
{
LocalValuesKey key{element.type().id(),local};
return ShapeValuesCache::get(key, [&](LocalValuesKey const&)
{
ShapeValues data;
localBasis_->evaluateFunction(local, data);
return data;
});
}
// evaluate local basis-gradients at point
template <class Element>
ShapeGradients const& evaluateJacobian(Element const& element, DomainType const& local) const
{
LocalGradientsKey key{element.type().id(),local};
return ShapeGradientsCache::get(key, [&](LocalGradientsKey const&)
{
ShapeGradients data;
localBasis_->evaluateJacobian(local, data);
return data;
}); });
} }
private: private:
using ShapeValuesCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, Policy>; using ShapeValuesCache = ConcurrentCache<LocalValuesKey, ShapeValues, ThreadLocalPolicy,
using ShapeGradientsCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, Policy>; std::unordered_map<LocalValuesKey, ShapeValues, typename LocalValuesKey::hasher>>;
using ShapeGradientsCache = ConcurrentCache<LocalGradientsKey, ShapeGradients, ThreadLocalPolicy,
std::unordered_map<LocalGradientsKey, ShapeGradients, typename LocalGradientsKey::hasher>>;
using ShapeValuesContainerCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, ThreadLocalPolicy,
std::unordered_map<ValuesKey, ShapeValuesContainer, typename ValuesKey::hasher>>;
using ShapeGradientsContainerCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, ThreadLocalPolicy,
std::unordered_map<GradientsKey, ShapeGradientsContainer, typename GradientsKey::hasher>>;
LocalBasisType const* localBasis_; LocalBasisType const* localBasis_;
}; };
......
#pragma once
#include <vector>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/utility/ConcurrentCache.hpp>
namespace AMDiS
{
namespace Impl
{
/// Pair of GeometryType and quadrature order and size
template <class LocalBasisType, class Tag>
struct EvaluatorCacheKey
{
using type = std::tuple<HashableGeometryType,typename LocalBasisType::Traits::DomainType>;
type value;
friend std::size_t hash_value(EvaluatorCacheKey const& t)
{
std::size_t seed = hash_value(std::get<0>(t.value));
Dune::hash_range(seed, std::get<1>(t.value).begin(), std::get<1>(t.value).end());
return seed;
}
friend bool operator==(EvaluatorCacheKey const& lhs, EvaluatorCacheKey const& rhs)
{
return lhs.value == rhs.value;
}
};
}
} // end namespace AMDiS
DUNE_DEFINE_HASH(DUNE_HASH_TEMPLATE_ARGS(class LocalBasisType, class Tag),DUNE_HASH_TYPE(AMDiS::Impl::EvaluatorCacheKey<LocalBasisType,Tag>))
namespace AMDiS
{
template <class LocalBasisType>
class LocalBasisEvaluatorCache
{
public:
using Traits = typename LocalBasisType::Traits;
using DomainType = typename Traits::DomainType;
using RangeType = typename Traits::RangeType;
using RangeFieldType = typename Traits::RangeFieldType;
using JacobianType = typename Traits::JacobianType;
using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
private:
using ValuesKey = Impl::EvaluatorCacheKey<LocalBasisType, struct ValuesTag>;
using GradientsKey = Impl::EvaluatorCacheKey<LocalBasisType, struct GradientsTag>;
using Policy = tag::thread_local_policy;
public:
using ShapeValues = std::vector<RangeType>;
using ShapeGradients = std::vector<JacobianType>;
LocalBasisEvaluatorCache(LocalBasisType const& localBasis)
: localBasis_(&localBasis)
{}
template <class LocalContext>
ShapeValues const& evaluateFunction(LocalContext const& context, DomainType const& local)
{
ValuesKey key{typename ValuesKey::type{context.type(),local}};
return ShapeValuesCache::get(key, [&](ShapeValues* data, ValuesKey const&)
{
localBasis_->evaluateFunction(local, *data);
});
}
template <class LocalContext>
ShapeGradients const& evaluateJacobian(LocalContext const& context, DomainType const& local)
{
GradientsKey key{typename GradientsKey::type{context.type(),local}};
return ShapeGradientsCache::get(key, [&](ShapeGradients* data, GradientsKey const&)
{
localBasis_->evaluateJacobian(local, *data);
});
}
private:
using ShapeValuesCache = ConcurrentCache<ValuesKey, ShapeValues, Policy>;
using ShapeGradientsCache = ConcurrentCache<GradientsKey, ShapeGradients, Policy>;
LocalBasisType const* localBasis_;
};
} // end namespace AMDiS
...@@ -56,16 +56,14 @@ namespace AMDiS ...@@ -56,16 +56,14 @@ namespace AMDiS
"Galerkin operator requires equal finite elements for test and trial space." ); "Galerkin operator requires equal finite elements for test and trial space." );
using LocalBasisType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType; using LocalBasisType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType;
using RangeType = typename LocalBasisType::Traits::RangeType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
using JacobianType = typename LocalBasisType::Traits::JacobianType;
auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element()); auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element());
auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element()); auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element());
auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element()); auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element());
auto const& localFE = rowNode.finiteElement(); auto const& localFE = rowNode.finiteElement();
std::size_t numLocalFe = localFE.size();
int quadDeg = std::max({this->getDegree(2,coeffOrder(localFctA),rowNode,colNode), int quadDeg = std::max({this->getDegree(2,coeffOrder(localFctA),rowNode,colNode),
this->getDegree(1,coeffOrder(localFctB),rowNode,colNode), this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
...@@ -74,23 +72,25 @@ namespace AMDiS ...@@ -74,23 +72,25 @@ namespace AMDiS
using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>; using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>;
auto const& quad = QuadratureRules::rule(context.type(), quadDeg); auto const& quad = QuadratureRules::rule(context.type(), quadDeg);
for (auto const& qp : quad) { LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // 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 // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // 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();
// the values of the shape functions on the reference element at the quadrature point // the values of the shape functions on the reference element at the quadrature point
thread_local std::vector<RangeType> shapeValues; auto const& shapeValues = shapeValuesCache[iq];
localFE.localBasis().evaluateFunction(local, shapeValues);
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
using WorldVector = FieldVector<RangeFieldType,Context::dow>; using WorldVector = FieldVector<RangeFieldType,Context::dow>;
...@@ -106,22 +106,22 @@ namespace AMDiS ...@@ -106,22 +106,22 @@ namespace AMDiS
IF_CONSTEXPR(conserving) { IF_CONSTEXPR(conserving) {
WorldVector gradAi, gradBi; WorldVector gradAi, gradBi;
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < numLocalFe; ++i) {
const int local_i = rowNode.localIndex(i); const int local_i = rowNode.localIndex(i);
gradAi = A * gradients[i]; gradAi = A * gradients[i];
gradBi = b * gradients[i]; gradBi = b * gradients[i];
for (std::size_t j = 0; j < localFE.size(); ++j) { for (std::size_t j = 0; j < numLocalFe; ++j) {
const int local_j = colNode.localIndex(j); const int local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor; elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
} }
} }
} else { } else {
WorldVector grad_i; WorldVector grad_i;
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < numLocalFe; ++i) {
const int local_i = rowNode.localIndex(i); const int local_i = rowNode.localIndex(i);
grad_i = A * gradients[i]; grad_i = A * gradients[i];
grad_i+= b * shapeValues[i]; grad_i+= b * shapeValues[i];
for (std::size_t j = 0; j < localFE.size(); ++j) { for (std::size_t j = 0; j < numLocalFe; ++j) {
const int local_j = colNode.localIndex(j); const int local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor; elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor;
} }
...@@ -143,33 +143,34 @@ namespace AMDiS ...@@ -143,33 +143,34 @@ namespace AMDiS
static_assert(Node::isLeaf, static_assert(Node::isLeaf,
"Operator can be applied to Leaf-Nodes only." ); "Operator can be applied to Leaf-Nodes only." );
using RangeType = typename FiniteElementType_t<Node>::Traits::LocalBasisType::Traits::RangeType; using LocalBasisType = typename FiniteElementType_t<Node>::Traits::LocalBasisType;
auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element()); auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
auto const& localFE = node.finiteElement(); auto const& localFE = node.finiteElement();
std::size_t numLocalFe = localFE.size();
int quad_order = this->getDegree(0,coeffOrder(localFctF),node); int quad_order = this->getDegree(0,coeffOrder(localFctF),node);
using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>; using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>;
auto const& quad = QuadratureRules::rule(context.type(), quad_order); auto const& quad = QuadratureRules::rule(context.type(), quad_order);
for (auto const& qp : quad) { LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
// Position of the current quadrature point in the reference element auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad);
decltype(auto) local = context.local(qp.position());
// The multiplicative factor in the integral transformation formula for (std::size_t iq = 0; iq < quad.size(); ++iq) {
const auto factor = context.integrationElement(qp.position()) * qp.weight(); // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(quad[iq].position());
// the values of the shape functions on the reference element at the quadrature point // the values of the shape functions on the reference element at the quadrature point
std::vector<RangeType> shapeValues; auto const& shapeValues = shapeValuesCache[iq];
localFE.localBasis().evaluateFunction(local, shapeValues);
const auto f = localFctF(local); // The multiplicative factor in the integral transformation formula
const auto factor = localFctF(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < numLocalFe; ++i) {
const int local_i = node.localIndex(i); const int local_i = node.localIndex(i);
elementVector[local_i] += f * shapeValues[i] * factor; elementVector[local_i] += shapeValues[i] * factor;
} }
} }
......
...@@ -119,7 +119,7 @@ namespace AMDiS ...@@ -119,7 +119,7 @@ namespace AMDiS
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.gradients(context, quad); auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(quad[iq].position()); decltype(auto) local = context.local(quad[iq].position());
......
...@@ -48,7 +48,10 @@ namespace AMDiS ...@@ -48,7 +48,10 @@ namespace AMDiS
using namespace Dune::Indices; using namespace Dune::Indices;
auto const& velocityLocalFE = tree.child(_0,0).finiteElement(); auto const& velocityLocalFE = tree.child(_0,0).finiteElement();