Commit 4db9083a authored by Praetorius, Simon's avatar Praetorius, Simon

Corrected errors in ConvectionDiffusionOperator and added example

parent 36d6a365
add_custom_target(examples) add_custom_target(examples)
set(projects2d "ellipt" "heat" "vecellipt" "stokes0" "stokes1" "stokes3" "navier_stokes") set(projects2d "ellipt" "heat" "vecellipt" "stokes0" "stokes1" "stokes3" "navier_stokes" "convection_diffusion")
foreach(project ${projects2d}) foreach(project ${projects2d})
add_executable(${project}.2d ${project}.cc) add_executable(${project}.2d ${project}.cc)
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <iostream>
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/assembler/ConvectionDiffusionOperator.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
// 1 component with polynomial degree 1
//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using ElliptParam = YaspGridBasis<AMDIS_DIM, 1>;
using ElliptProblem = ProblemStat<ElliptParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using namespace Dune::Indices;
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
// -div(A*grad(u)) + div(b*u) + c*u = f
auto opCD = convectionDiffusion(/*A=*/1.0, /*b=*/0.0, /*c=*/1.0, /*f=*/1.0);
prob.addMatrixOperator(opCD, _0, _0);
prob.addVectorOperator(opCD, _0);
// set boundary condition
auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
auto dbcValues = [](auto const& x){ return 0.0; }; // set value
prob.addDirichletBC(predicate, _0, _0, dbcValues);
AdaptInfo adaptInfo("adapt");
prob.buildAfterCoarsen(adaptInfo, Flag(0));
prob.solve(adaptInfo);
prob.writeFiles(adaptInfo, true);
AMDiS::finalize();
return 0;
}
...@@ -14,7 +14,7 @@ namespace AMDiS ...@@ -14,7 +14,7 @@ namespace AMDiS
* @{ * @{
**/ **/
/// convection-diffusion operator /// convection-diffusion operator, see \ref convectionDiffusion
/// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or /// <A*grad(u),grad(v)> - <b*u, grad(v)> + <c*u, v> = <f, v> (conserving) or
/// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving) /// <A*grad(u),grad(v)> + <b*grad(u), v> + <c*u, v> = <f, v> (non conserving)
template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true> template <class LocalContext, class GridFctA, class GridFctB, class GridFctC, class GridFctF, bool conserving = true>
...@@ -71,7 +71,7 @@ namespace AMDiS ...@@ -71,7 +71,7 @@ namespace AMDiS
this->getDegree(1,coeffOrder(localFctB),rowNode,colNode), this->getDegree(1,coeffOrder(localFctB),rowNode,colNode),
this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)}); this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)});
using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::dimension>; 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) { for (auto const& qp : quad) {
...@@ -93,30 +93,33 @@ namespace AMDiS ...@@ -93,30 +93,33 @@ namespace AMDiS
localFE.localBasis().evaluateJacobian(local, shapeGradients); localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients(shapeGradients.size()); using WorldVector = FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i) for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]); jacobian.mv(shapeGradients[i][0], gradients[i]);
const auto A = localFctA(local); const auto A = localFctA(local);
const auto b = localFctB(local); WorldVector b; b = localFctB(local);
const auto c = localFctC(local); const auto c = localFctC(local);
IF_CONSTEXPR(conserving) { IF_CONSTEXPR(conserving) {
WorldVector gradAi, gradBi;
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowNode.localIndex(i); const int local_i = rowNode.localIndex(i);
const auto gradAi = A * gradients[i]; gradAi = A * gradients[i];
const auto 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 < localFE.size(); ++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;
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowNode.localIndex(i); const int local_i = rowNode.localIndex(i);
const auto 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 < localFE.size(); ++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;
...@@ -136,10 +139,10 @@ namespace AMDiS ...@@ -136,10 +139,10 @@ namespace AMDiS
Node const& node, Node const& node,
ElementVector& elementVector) ElementVector& elementVector)
{ {
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<RowNode>::Traits::LocalBasisType::Traits::RangeType; using RangeType = typename FiniteElementType_t<Node>::Traits::LocalBasisType::Traits::RangeType;
auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element()); auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
...@@ -147,7 +150,7 @@ namespace AMDiS ...@@ -147,7 +150,7 @@ namespace AMDiS
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::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) { for (auto const& qp : quad) {
...@@ -164,8 +167,8 @@ namespace AMDiS ...@@ -164,8 +167,8 @@ namespace AMDiS
const auto f = localFctF(local); const auto f = localFctF(local);
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowNode.localIndex(i); const int local_i = node.localIndex(i);
elementVector[local_i][local_j] += f * shapeValues[i] * factor; elementVector[local_i] += f * shapeValues[i] * factor;
} }
} }
...@@ -180,22 +183,23 @@ namespace AMDiS ...@@ -180,22 +183,23 @@ namespace AMDiS
template <class LocalFct> template <class LocalFct>
int coeffOrder(LocalFct const& localFct) int coeffOrder(LocalFct const& localFct)
{ {
using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>; using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFct>;
return Dune::Hybrid::ifElse(Concept{}, return Dune::Hybrid::ifElse(Concept{},
[&](auto id) { return order(id(localFct)); }, [&](auto id) { return order(id(localFct)); },
[] (auto) { return 0; }); [] (auto) { return 0; });
} }
private: private:
GridFctA gridFctA; GridFctA gridFctA_;
GridFctB gridFctB; GridFctB gridFctB_;
GridFctC gridFctC; GridFctC gridFctC_;
GridFctF gridFctF; GridFctF gridFctF_;
}; };
template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, bool conserving> template <class PreGridFctA, class PreGridFctB, class PreGridFctC, class PreGridFctF, class c>
struct PreConvectionDiffusionOperator struct PreConvectionDiffusionOperator
{ {
static constexpr bool conserving = c::value;
PreGridFctA gridFctA; PreGridFctA gridFctA;
PreGridFctB gridFctB; PreGridFctB gridFctB;
PreGridFctC gridFctC; PreGridFctC gridFctC;
...@@ -207,20 +211,21 @@ namespace AMDiS ...@@ -207,20 +211,21 @@ namespace AMDiS
PreGridFctC const& gridFctC, PreGridFctF const& gridFctF, PreGridFctC const& gridFctC, PreGridFctF const& gridFctF,
bool_t<conserving> = {}) bool_t<conserving> = {})
{ {
using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, conserving>; using Pre = PreConvectionDiffusionOperator<PreGridFctA, PreGridFctB, PreGridFctC, PreGridFctF, bool_t<conserving>>;
return Pre{gridFctA, gridFctB, gridFctC, gridFctF}; return Pre{gridFctA, gridFctB, gridFctC, gridFctF};
} }
template <class LocalContext, class... GrdFcts, bool conserving, class GridView> template <class LocalContext, class... T, class GridView>
auto makeLocalOperator(PreConvectionDiffusionOperator<GridFcts..., conserving> const& pre, GridView const& gridView) auto makeLocalOperator(PreConvectionDiffusionOperator<T...> const& pre, GridView const& gridView)
{ {
using Pre = PreConvectionDiffusionOperator<T...>;
auto gridFctA = makeGridFunction(pre.gridFctA, gridView); auto gridFctA = makeGridFunction(pre.gridFctA, gridView);
auto gridFctB = makeGridFunction(pre.gridFctB, gridView); auto gridFctB = makeGridFunction(pre.gridFctB, gridView);
auto gridFctC = makeGridFunction(pre.gridFctC, gridView); auto gridFctC = makeGridFunction(pre.gridFctC, gridView);
auto gridFctF = makeGridFunction(pre.gridFctF, gridView); auto gridFctF = makeGridFunction(pre.gridFctF, gridView);
using GridFctOp = ConvectionDiffusionOperator<LocalContext, using GridFctOp = ConvectionDiffusionOperator<LocalContext,
decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), conserving>; decltype(gridFctA), decltype(gridFctB), decltype(gridFctC), decltype(gridFctF), Pre::conserving>;
GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF}; GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF};
return localOperator; return localOperator;
......
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