From fb276ed0390a350e29c17459e5512a7c1ef6e808 Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 29 Sep 2022 17:46:44 +0200 Subject: [PATCH 1/2] [disc][cvfe] Gradient calculation for cvfe schemes --- dumux/discretization/evalgradients.hh | 96 +++++++++++++++++++++++---- 1 file changed, 84 insertions(+), 12 deletions(-) diff --git a/dumux/discretization/evalgradients.hh b/dumux/discretization/evalgradients.hh index e891f65957..c3e6e423fd 100644 --- a/dumux/discretization/evalgradients.hh +++ b/dumux/discretization/evalgradients.hh @@ -31,13 +31,17 @@ #include #include #include +#include +#include #include "evalsolution.hh" namespace Dumux { +// some implementation details +namespace Detail{ /*! - * \brief Evaluates the gradient of a given box element solution to a given global position. + * \brief Evaluates the gradient of a control-volume finite element solution to a given global position. * \ingroup Discretization * * \param element The element @@ -51,13 +55,13 @@ namespace Dumux { * the PrimaryVariables object (i.e. numEq). Each entry is * a GlobalCoordinate object holding the priVar gradient. */ -template -auto evalGradients(const Element& element, - const typename Element::Geometry& geometry, - const typename FVElementGeometry::GridGeometry& gridGeometry, - const BoxElementSolution& elemSol, - const typename Element::Geometry::GlobalCoordinate& globalPos, - bool ignoreState = false) +template +auto evalCVFEGradients(const Element& element, + const typename Element::Geometry& geometry, + const GridGeometry& gridGeometry, + const CVFEElemSol& elemSol, + const typename Element::Geometry::GlobalCoordinate& globalPos, + bool ignoreState = false) { // determine if all states are the same at all vertices using HasState = decltype(isValid(Detail::hasState())(elemSol[0])); @@ -80,15 +84,15 @@ auto evalGradients(const Element& element, const auto jacInvT = geometry.jacobianInverseTransposed(localPos); // interpolate the gradients - Dune::FieldVector result( GlobalPosition(0.0) ); - for (int i = 0; i < element.subEntities(Element::Geometry::mydimension); ++i) + Dune::FieldVector result( GlobalPosition(0.0) ); + for (int i = 0; i < shapeJacobian.size(); ++i) { // the global shape function gradient GlobalPosition gradN; jacInvT.mv(shapeJacobian[i][0], gradN); // add gradient to global privar gradients - for (unsigned int pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx) + for (unsigned int pvIdx = 0; pvIdx < CVFEElemSol::PrimaryVariables::dimension; ++pvIdx) { GlobalPosition tmp(gradN); tmp *= elemSol[i][pvIdx]; @@ -100,10 +104,78 @@ auto evalGradients(const Element& element, } else { - DUNE_THROW(Dune::NotImplemented, "Element vertices have different phase states. Enforce calculation by setting ignoreState to true."); + DUNE_THROW(Dune::NotImplemented, "Element dofs have different phase states. Enforce calculation by setting ignoreState to true."); } } +} // end namespace Detail + +/*! + * \brief Evaluates the gradient of a given box element solution to a given global position. + * \ingroup Discretization + * + * \param element The element + * \param geometry The element geometry + * \param gridGeometry The finite volume grid geometry + * \param elemSol The primary variables at the dofs of the element + * \param globalPos The global position + * \param ignoreState If true, the state of primary variables is ignored + */ +template +auto evalGradients(const Element& element, + const typename Element::Geometry& geometry, + const typename FVElementGeometry::GridGeometry& gridGeometry, + const BoxElementSolution& elemSol, + const typename Element::Geometry::GlobalCoordinate& globalPos, + bool ignoreState = false) +{ + return Detail::evalCVFEGradients(element, geometry, gridGeometry, elemSol, globalPos, ignoreState); +} + +/*! + * \brief Evaluates the gradient of a given diamond scheme element solution to a given global position. + * \ingroup Discretization + * + * \param element The element + * \param geometry The element geometry + * \param gridGeometry The finite volume grid geometry + * \param elemSol The primary variables at the dofs of the element + * \param globalPos The global position + * \param ignoreState If true, the state of primary variables is ignored + */ +template +auto evalGradients(const Element& element, + const typename Element::Geometry& geometry, + const typename FVElementGeometry::GridGeometry& gridGeometry, + const FaceCenteredDiamondElementSolution& elemSol, + const typename Element::Geometry::GlobalCoordinate& globalPos, + bool ignoreState = false) +{ + return Detail::evalCVFEGradients(element, geometry, gridGeometry, elemSol, globalPos, ignoreState); +} + +/*! + * \brief Evaluates the gradient of a given pq1bubble scheme element solution to a given global position. + * \ingroup Discretization + * + * \param element The element + * \param geometry The element geometry + * \param gridGeometry The finite volume grid geometry + * \param elemSol The primary variables at the dofs of the element + * \param globalPos The global position + * \param ignoreState If true, the state of primary variables is ignored + */ +template +auto evalGradients(const Element& element, + const typename Element::Geometry& geometry, + const typename FVElementGeometry::GridGeometry& gridGeometry, + const PQ1BubbleElementSolution& elemSol, + const typename Element::Geometry::GlobalCoordinate& globalPos, + bool ignoreState = false) +{ + return Detail::evalCVFEGradients(element, geometry, gridGeometry, elemSol, globalPos, ignoreState); +} + /*! * \ingroup Discretization * \brief Evaluates the gradient of a given CCElementSolution to a given global position. -- GitLab From 3cccb97b2b003d2c41b456d34f4dc0ca7642a3fa Mon Sep 17 00:00:00 2001 From: Martin Schneider Date: Thu, 29 Sep 2022 18:00:50 +0200 Subject: [PATCH 2/2] [test][ff][unstructured] Use evalGradients interface to calculate stress --- .../navierstokes/unstructured/problem.hh | 64 +++++-------------- 1 file changed, 16 insertions(+), 48 deletions(-) diff --git a/test/freeflow/navierstokes/unstructured/problem.hh b/test/freeflow/navierstokes/unstructured/problem.hh index d769b9e55a..fa18189fd5 100644 --- a/test/freeflow/navierstokes/unstructured/problem.hh +++ b/test/freeflow/navierstokes/unstructured/problem.hh @@ -32,6 +32,7 @@ #include #include #include +#include #include namespace Dumux { @@ -210,7 +211,6 @@ public: { auto fvGeometry = localView(this->gridGeometry()); auto elemVolVars = localView(gridVariables.curGridVolVars()); - auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache()); BoundaryFluxes forces(0.0); Scalar cylinderSurface = 0.0; @@ -223,24 +223,17 @@ public: { if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal())) { - elemVolVars.bind(element, fvGeometry, v); - elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); - const auto& fluxVarCache = elemFluxVarsCache[scvf]; - - Dune::FieldMatrix stress(0.0); - for (const auto& scv : scvs(fvGeometry)) - { - const auto velocity = elemVolVars[scv].velocity(); - for (int dir = 0; dir < dim; ++dir) - stress[dir].axpy(velocity[dir], fluxVarCache.gradN(scv.indexInElement())); - } - + const auto elemSol = elementSolution(element, elemVolVars, fvGeometry); + auto stress = evalGradients(element, element.geometry(), fvGeometry.gridGeometry(), elemSol, scvf.ipGlobal()); stress *= viscosity_; + BoundaryFluxes normalStress(0.0); for (int dir = 0; dir < dim; ++dir) + { stress[dir][dir] -= this->pressure(element, fvGeometry, scvf); + normalStress[dir] = stress[dir]*scvf.unitOuterNormal(); + } - BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal()); forces.axpy(scvf.area(), normalStress); cylinderSurface += scvf.area(); @@ -373,54 +366,29 @@ public: if (scvf.boundary() && onCylinderBoundary_(scvf.ipGlobal())) { elemVolVars.bind(element, fvGeometry, v); - Dune::FieldMatrix stress(0.0); - const auto& localBasis = fvGeometry.feLocalBasis(); - Scalar area = 0.0; + Dune::FieldVector stress(GlobalPosition(0.0)); const auto geometry = element.geometry(); const auto isgeometry = fvGeometry.geometry(scvf); const auto& quad = Dune::QuadratureRules::mydimension>::rule(isgeometry.type(), 5); for (auto&& qp : quad) { - Dune::FieldMatrix stressLocal(0.0); - using FeLocalBasis = typename GridGeometry::FeCache::FiniteElementType::Traits::LocalBasisType; - using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType; - using ShapeValue = typename Dune::FieldVector; - using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed; - std::vector gradN; - std::vector shapeJacobian; - std::vector shapeValues; const auto ipGlobal = isgeometry.global(qp.position()); - const auto ipLocal = geometry.local(ipGlobal); - JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(ipLocal); - localBasis.evaluateJacobian(ipLocal, shapeJacobian); - localBasis.evaluateFunction(ipLocal, shapeValues); - - // compute the gradN at for every scv/dof - gradN.resize(fvGeometry.numScv()); - for (const auto& scv: scvs(fvGeometry)) - jacInvT.mv(shapeJacobian[scv.localDofIndex()][0], gradN[scv.indexInElement()]); - - for (const auto& scv : scvs(fvGeometry)) - { - const auto velocity = elemVolVars[scv].velocity(); - for (int dir = 0; dir < dim; ++dir) - stressLocal[dir].axpy(velocity[dir], gradN[scv.indexInElement()]); - } - area += qp.weight()*isgeometry.integrationElement(qp.position()); - stressLocal += getTransposed(stressLocal); + const auto elemSol = elementSolution(element, elemVolVars, fvGeometry); + auto stressLocal = evalGradients(element, geometry, fvGeometry.gridGeometry(), elemSol, ipGlobal); stressLocal *= qp.weight()*isgeometry.integrationElement(qp.position()); - - stress += stressLocal; + stress = stress + stressLocal; } stress *= viscosity_; stress /= scvf.area(); - // Add pressure contribuition + BoundaryFluxes normalStress(0.0); for (int dir = 0; dir < dim; ++dir) + { + // Add pressure contribuition stress[dir][dir] -= this->pressure(element, fvGeometry, scvf); - - BoundaryFluxes normalStress = mv(stress, scvf.unitOuterNormal()); + normalStress[dir] = stress[dir]*scvf.unitOuterNormal(); + } forces.axpy(scvf.area(), normalStress); cylinderSurface += scvf.area(); -- GitLab