Commit 18ecdb30 authored by Praetorius, Simon's avatar Praetorius, Simon

Intel compiler compatibility

parent ce0fece2
Pipeline #953 failed with stage
in 2 minutes and 6 seconds
...@@ -60,13 +60,22 @@ namespace AMDiS ...@@ -60,13 +60,22 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const; bool asmMatrix, bool asmVector) const;
template <class ElementContainer, class Container, class Operators, class Geometry, class... Bases> template <class ElementContainer, class Container, class Operators, class Geometry, class Basis>
void assembleElementOperators( void assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
Geometry const& geometry, Geometry const& geometry,
Bases const&... subBases) const; Basis const& subBasis) const;
template <class ElementContainer, class Container, class Operators, class Geometry, class RowBasis, class ColBasis>
void assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
RowBasis const& rowBasis,
ColBasis const& colBasis) const;
/// Finish insertion into the matrix and assembles boundary conditions /// Finish insertion into the matrix and assembles boundary conditions
......
...@@ -94,23 +94,23 @@ void Assembler<Traits>::assemble( ...@@ -94,23 +94,23 @@ void Assembler<Traits>::assemble(
template <class Traits> template <class Traits>
template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews> template <class ElementContainer, class Container, class Operators, class Geometry, class LocalView>
void Assembler<Traits>::assembleElementOperators( void Assembler<Traits>::assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
Geometry const& geometry, Geometry const& geometry,
LocalViews const&... localViews) const LocalView const& localView) const
{ {
auto const& element = getElement(localViews...); auto const& element = getElement(localView);
auto const& gridView = getGridView(localViews...); auto const& gridView = getGridView(localView);
bool add = false; bool add = false;
auto assemble_operators = [&](auto const& context, auto& operator_list) { auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) { for (auto scaled : operator_list) {
scaled.op->bind(element, geometry); scaled.op->bind(element, geometry);
bool add_op = scaled.op->assemble(context, elementContainer, localViews.tree()...); bool add_op = scaled.op->assemble(context, elementContainer, localView.tree());
scaled.op->unbind(); scaled.op->unbind();
add = add || add_op; add = add || add_op;
} }
...@@ -133,6 +133,45 @@ void Assembler<Traits>::assembleElementOperators( ...@@ -133,6 +133,45 @@ void Assembler<Traits>::assembleElementOperators(
} }
template <class Traits>
template <class ElementContainer, class Container, class Operators, class Geometry, class RowLocalViews, class ColLocalView>
void Assembler<Traits>::assembleElementOperators(
ElementContainer& elementContainer,
Container& container,
Operators& operators,
Geometry const& geometry,
RowLocalViews const& rowLocalView, ColLocalView const& colLocalView) const
{
auto const& element = getElement(rowLocalView, colLocalView);
auto const& gridView = getGridView(rowLocalView, colLocalView);
bool add = false;
auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) {
scaled.op->bind(element, geometry);
bool add_op = scaled.op->assemble(context, elementContainer, rowLocalView.tree(), colLocalView.tree());
scaled.op->unbind();
add = add || add_op;
}
};
// assemble element operators
assemble_operators(element, operators.element);
// assemble intersection operators
if (!operators.intersection.empty()
|| (!operators.boundary.empty() && element.hasBoundaryIntersections()))
{
for (auto const& intersection : intersections(gridView, element)) {
if (intersection.boundary())
assemble_operators(intersection, operators.boundary);
else
assemble_operators(intersection, operators.intersection);
}
}
}
template <class Traits> template <class Traits>
template <class SystemMatrixType, class SystemVectorType> template <class SystemMatrixType, class SystemVectorType>
void Assembler<Traits>::initMatrixVector( void Assembler<Traits>::initMatrixVector(
......
...@@ -265,7 +265,7 @@ namespace AMDiS ...@@ -265,7 +265,7 @@ namespace AMDiS
{ {
using RawExpr = Underlying_t<Expr>; using RawExpr = Underlying_t<Expr>;
static_assert(Concepts::HasGridFunctionOrder<RawExpr>, static_assert(Concepts::HasGridFunctionOrder<RawExpr>,
"Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in `makeOperator()`."); "Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in 'makeOperator()'.");
return ExpressionPreOperator<Tag, Expr, tag::deduce>{t, expr}; return ExpressionPreOperator<Tag, Expr, tag::deduce>{t, expr};
} }
......
...@@ -137,6 +137,11 @@ void ProblemStat<Traits>::addMatrixOperator( ...@@ -137,6 +137,11 @@ void ProblemStat<Traits>::addMatrixOperator(
RowTreePath row, ColTreePath col, RowTreePath row, ColTreePath col,
double* factor, double* estFactor) double* factor, double* estFactor)
{ {
static_assert( Concept::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
static_assert( Concept::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row)); auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col)); auto j = child(globalBasis->localView().tree(), makeTreePath(col));
...@@ -156,6 +161,11 @@ void ProblemStat<Traits>::addMatrixOperator( ...@@ -156,6 +161,11 @@ void ProblemStat<Traits>::addMatrixOperator(
RowTreePath row, ColTreePath col, RowTreePath row, ColTreePath col,
double* factor, double* estFactor) double* factor, double* estFactor)
{ {
static_assert( Concept::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
static_assert( Concept::PreTreePath<ColTreePath>,
"col must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(row)); auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col)); auto j = child(globalBasis->localView().tree(), makeTreePath(col));
using Intersection = typename GridView::Intersection; using Intersection = typename GridView::Intersection;
...@@ -175,6 +185,9 @@ void ProblemStat<Traits>::addVectorOperator( ...@@ -175,6 +185,9 @@ void ProblemStat<Traits>::addVectorOperator(
TreePath path, TreePath path,
double* factor, double* estFactor) double* factor, double* estFactor)
{ {
static_assert( Concept::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path)); auto i = child(globalBasis->localView().tree(), makeTreePath(path));
auto op = makeGridOperator(preOp, globalBasis->gridView()); auto op = makeGridOperator(preOp, globalBasis->gridView());
...@@ -193,6 +206,9 @@ void ProblemStat<Traits>::addVectorOperator( ...@@ -193,6 +206,9 @@ void ProblemStat<Traits>::addVectorOperator(
TreePath path, TreePath path,
double* factor, double* estFactor) double* factor, double* estFactor)
{ {
static_assert( Concept::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
auto i = child(globalBasis->localView().tree(), makeTreePath(path)); auto i = child(globalBasis->localView().tree(), makeTreePath(path));
using Intersection = typename GridView::Intersection; using Intersection = typename GridView::Intersection;
......
...@@ -116,7 +116,7 @@ namespace AMDiS ...@@ -116,7 +116,7 @@ namespace AMDiS
} }
/// \brief Return the LocalFunction of the AnalyticGridFunction. /// \brief Return the LocalFunction of the AnalyticGridFunction.
friend auto localFunction(AnalyticGridFunction const& gf) friend LocalFunction localFunction(AnalyticGridFunction const& gf)
{ {
return LocalFunction{gf.fct_}; return LocalFunction{gf.fct_};
} }
......
...@@ -47,19 +47,23 @@ namespace AMDiS ...@@ -47,19 +47,23 @@ namespace AMDiS
return 0; return 0;
} }
/// \relates ConstantLocalFunction template <class R_, class D_, class L_, class T_>
friend auto derivative(ConstantLocalFunction const& lf) friend auto derivative(ConstantLocalFunction<R_(D_), L_, T_> const& lf);
{
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
DerivativeRange diff(0);
return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
}
private: private:
T value_; T value_;
}; };
/// \relates ConstantLocalFunction
template <class R, class D, class LocalContext, class T>
auto derivative(ConstantLocalFunction<R(D), LocalContext, T> const& lf)
{
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
DerivativeRange diff(0);
return ConstantLocalFunction<DerivativeRange(D),LocalContext,DerivativeRange>{diff};
}
/** /**
* \addtogroup GridFunctions * \addtogroup GridFunctions
...@@ -101,7 +105,7 @@ namespace AMDiS ...@@ -101,7 +105,7 @@ namespace AMDiS
} }
/// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction /// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction
friend auto localFunction(ConstantGridFunction const& gf) friend LocalFunction localFunction(ConstantGridFunction const& gf)
{ {
return LocalFunction{gf.value_}; return LocalFunction{gf.value_};
} }
......
...@@ -8,6 +8,68 @@ ...@@ -8,6 +8,68 @@
namespace AMDiS namespace AMDiS
{ {
namespace Concepts
{
/** \addtogroup Concepts
* @{
**/
namespace Definition
{
template <class TP>
struct IsPreTreePath
: any_of_t<std::is_same<TP, int>::value, std::is_same<TP,std::size_t>::value>
{};
template <int I>
struct IsPreTreePath<int_t<I>>
: std::true_type
{};
template <std::size_t I>
struct IsPreTreePath<index_t<I>>
: std::true_type
{};
template <std::size_t... I>
struct IsPreTreePath<Indices<I...>>
: std::true_type
{};
} // end namespace Definition
template <class TP>
constexpr bool PreTreePath = Dune::TypeTree::IsTreePath<TP>::value || Definition::IsPreTreePath<TP>::value;
/** @} **/
} // end namespace Concepts
#ifdef DOXYGEN
/// \brief Converts a (pre)TreePath into a HybridTreePath
/**
* Converts an integer, an integralconstant or a Dune TreePath into an
* \ref Dune::TypeTree::HybridTreePath that is used in GlobalBasis traversal.
*
* **Requirements:**
* - `PreTreePath` one of
* + integer type (`int, std::size_t`),
* + integral constant (`std::integral_constant<[int|std::size_t], i>`)
* + any Dune TreePath (`TreePath<std::size_t...>, DynamicTreePath, HybridTreePath<class... T>`)
*
* **Example:**
* ```
* makeTreePath(1),
* makeTreePath(int_<2>),
* makeTreePath(treepath(1, int_<2>))
* ```
**/
template <class PreTreePath>
auto makeTreePath(PreTreePath tp);
#else // DOXYGEN
auto makeTreePath(int i) { return Dune::TypeTree::hybridTreePath(std::size_t(i)); } auto makeTreePath(int i) { return Dune::TypeTree::hybridTreePath(std::size_t(i)); }
auto makeTreePath(std::size_t i) { return Dune::TypeTree::hybridTreePath(i); } auto makeTreePath(std::size_t i) { return Dune::TypeTree::hybridTreePath(i); }
...@@ -28,61 +90,46 @@ namespace AMDiS ...@@ -28,61 +90,46 @@ namespace AMDiS
template <class TP> template <class TP>
auto makeTreePath(TP const&) { auto makeTreePath(TP const&) {
static_assert( Dune::TypeTree::IsTreePath<TP>::value, static_assert( Concepts::PreTreePath<TP>,
"Argument must be a valid treepath, or an integer/index-constant"); "Argument must be a valid treepath, or an integer/index-constant");
return Dune::TypeTree::hybridTreePath(); return Dune::TypeTree::hybridTreePath();
} }
#endif // DOXYGEN
namespace Impl namespace Impl
{ {
template <std::size_t I, std::size_t N> template <class TreePath, std::size_t... I>
struct PrintTreePath void printTreePath(std::ostream& os, TreePath const& tp, std::index_sequence<I...>)
{
template <class... T>
static void apply(std::ostream& os, Dune::TypeTree::HybridTreePath<T...> const& tp)
{
os << "," << treePathIndex(tp, index_<I>);
PrintTreePath<I+1,N>::apply(os, tp);
}
};
template <std::size_t N>
struct PrintTreePath<N,N>
{
template <class... T>
static void apply(std::ostream&, Dune::TypeTree::HybridTreePath<T...> const&) {}
};
template <std::size_t N>
struct PrintTreePath<0,N>
{
template <class... T>
static void apply(std::ostream& os, Dune::TypeTree::HybridTreePath<T...> const& tp)
{
os << treePathIndex(tp, _0);
PrintTreePath<1,N>::apply(os, tp);
}
};
template <>
struct PrintTreePath<0,0>
{ {
template <class... T> (void)std::initializer_list<int>{
static void apply(std::ostream&, Dune::TypeTree::HybridTreePath<T...> const&) {} ((void)(os << treePathIndex(tp, index_<I>) << ","), 0)...
}; };
}
} }
template <class... T> template <class T0, class... T>
std::string to_string(Dune::TypeTree::HybridTreePath<T...> const& tp) std::string to_string(Dune::TypeTree::HybridTreePath<T0,T...> const& tp)
{ {
std::stringstream ss; std::stringstream ss;
Impl::PrintTreePath<0,sizeof...(T)>::apply(ss, tp); Impl::printTreePath(ss, tp, std::make_index_sequence<sizeof...(T)>{});
ss << back(tp);
return ss.str(); return ss.str();
} }
// alias for hybridTreePath(...) std::string to_string(Dune::TypeTree::HybridTreePath<> const& tp)
{
return "";
}
/// \brief Generate a TreePath from a sequence of integers and integral-constants
/**
* Converts a sequence of arguments to a \ref Dune::TypeTree::HybridTreePath.
* The arguments can be one of
* - integer type (`int, std::size_t`)
* - integral constant (`std::integral_constant<std::size_t,i>, index_t<i>`)
**/
template <class... T> template <class... T>
constexpr Dune::TypeTree::HybridTreePath<T...> treepath(T const&... t) constexpr Dune::TypeTree::HybridTreePath<T...> treepath(T const&... t)
{ {
......
...@@ -50,53 +50,42 @@ int main(int argc, char** argv) ...@@ -50,53 +50,42 @@ int main(int argc, char** argv)
auto _p = 1_c; auto _p = 1_c;
// <1/tau * u, v> // <1/tau * u, v>
Op opTime; auto opTime = makeOperator(tag::testvec_trialvec{}, density);
opTime.addZOT(density);
prob.addMatrixOperator(opTime, _v, _v); prob.addMatrixOperator(opTime, _v, _v);
// <1/tau * u^old, v> // <1/tau * u^old, v>
Op opTimeOld; auto opTimeOld = makeOperator(tag::testvec{}, density * prob.getSolution(_v));
opTimeOld.addZOT(valueOf(prob.getSolution(_v), density));
prob.addVectorOperator(opTimeOld, _v); prob.addVectorOperator(opTimeOld, _v);
// sum_i <grad(u_i),grad(v_i)> // sum_i <grad(u_i),grad(v_i)>
Op opL; for (std::size_t i = 0; i < DOW; ++i) {
opL.addSOT(viscosity); auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
prob.addMatrixOperator(opL, _v, _v); prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
}
// <p, div(v)> // <p, div(v)>
Op opP; auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
opP.addFOT(1.0, GRD_PSI);
prob.addMatrixOperator(opP, _v, _p); prob.addMatrixOperator(opP, _v, _p);
// <div(u), q> // <div(u), q>
Op opDiv; auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
opDiv.addFOT(1.0, GRD_PHI);
prob.addMatrixOperator(opDiv, _p, _v); prob.addMatrixOperator(opDiv, _p, _v);
for (std::size_t i = 0; i < DOW; ++i) { // <(u * nabla)u_i^old, v_i>
// <(u * nabla)u_i^old, v_i> auto opNonlin1 = makeOperator(tag::testvec_trialvec{}, density * trans(gradientAtQP(prob.getSolution(_v))));
for (std::size_t j = 0; j < DOW; ++j) { prob.addMatrixOperator(opNonlin1, _v, _v);
Op opNonlin1;
opNonlin1.addZOT(derivativeOf(prob.getSolution(treepath(_v,i)), j, density));
prob.addMatrixOperator(opNonlin1, treepath(_v,i), treepath(_v,j));
}
for (std::size_t i = 0; i < DOW; ++i) {
// <(u^old * nabla)u_i, v_i> // <(u^old * nabla)u_i, v_i>
Op opNonlin2; auto opNonlin2 = makeOperator(tag::test_gradtrial{}, density * prob.getSolution(_v));
opNonlin2.addFOT(valueOf(prob.getSolution(_v)), GRD_PHI);
prob.addMatrixOperator(opNonlin2, treepath(_v,i), treepath(_v,i)); prob.addMatrixOperator(opNonlin2, treepath(_v,i), treepath(_v,i));
// <(u^old * grad(u_i^old)), v_i>
Op opNonlin3;
opNonlin3.addZOT( func([density](Dune::FieldVector<double,DOW> const& u, Dune::FieldVector<double,DOW> const& grad_ui)
{
return density * (u * grad_ui);
}, valueOf(prob.getSolution(_v)), gradientOf(prob.getSolution(treepath(_v,i)))) );
prob.addVectorOperator(opNonlin, treepath(_v,i));
} }
// <(u^old * grad(u_i^old)), v_i>
auto opNonlin3 = makeOperator(tag::testvec{}, trans(gradientAtQP(prob.getSolution(_v))) * prob.getSolution(_v));
prob.addVectorOperator(opNonlin3, _v);
// define boundary regions // define boundary regions
auto left = [](auto const& x) { return x[0] < 1.e-8; }; auto left = [](auto const& x) { return x[0] < 1.e-8; };
auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; }; auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
...@@ -119,6 +108,9 @@ int main(int argc, char** argv) ...@@ -119,6 +108,9 @@ int main(int argc, char** argv)
// set point constraint for pressure // set point constraint for pressure
prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0); prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);
// set initial conditions
prob.getSolution(_v).interpolate(parabolic_y);
prob.getSolution(_p).interpolate(0.0);
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo); AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt(); adapt.adapt();
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment