From f72e58740f3c8071b5acd26f31e8a98506345445 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 3 Nov 2017 22:01:16 +0100
Subject: [PATCH] [staggered] First attempt to introduce localFaceVars
---
dumux/assembly/staggeredlocalassembler.hh | 41 ++++----
.../discretization/staggered/facesolution.hh | 93 +++++++++++++++++++
.../staggered/freeflow/facevariables.hh | 86 ++++++++++++++++-
.../staggered/globalfacevariables.hh | 32 ++++++-
dumux/discretization/staggered/properties.hh | 4 +
dumux/freeflow/staggered/fluxvariables.hh | 81 +++++++++-------
dumux/freeflow/staggered/localresidual.hh | 8 +-
dumux/freeflow/staggered/velocityoutput.hh | 4 +-
dumux/implicit/staggered/newtoncontroller.hh | 2 +-
9 files changed, 283 insertions(+), 68 deletions(-)
create mode 100644 dumux/discretization/staggered/facesolution.hh
diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
index f7ab28422e..8351df9199 100644
--- a/dumux/assembly/staggeredlocalassembler.hh
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -31,6 +31,7 @@
#include <dumux/assembly/diffmethod.hh>
#include <dumux/implicit/staggered/primaryvariables.hh>
+#include <dumux/discretization/staggered/facesolution.hh>
namespace Dumux {
@@ -77,9 +78,6 @@ class StaggeredLocalAssembler<TypeTag,
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
- // using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
- // typename DofTypeIndices::CellCenterIdx cellCenterIdx;
- // typename DofTypeIndices::FaceIdx faceIdx;
static constexpr bool enableGlobalFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache);
@@ -144,7 +142,6 @@ public:
for(auto&& scvf : scvfs(fvGeometry))
{
- // res[faceIdx][scvf.dofIndex()] = 0.0;
faceResidualCache[scvf.localFaceIdx()] = localResidual.evalFace(problem,
element,
fvGeometry,
@@ -371,25 +368,30 @@ static void dCCdFace_(Assembler& assembler,
const auto& problem = assembler.problem();
auto& localResidual = assembler.localResidual();
- // auto& gridVariables = assembler.gridVariables();
- // build derivatives with for cell center dofs w.r.t. cell center dofs
- const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
-
- const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+ // build derivatives with for cell center dofs w.r.t. cell center dofs
+ const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
- for(const auto& globalJ : connectivityMap(cellCenterIdx, faceIdx, cellCenterGlobalI))
+ for(const auto& scvfJ : scvfs(fvGeometry))
{
+ const auto globalJ = scvfJ.dofIndex();
+
// get the faceVars of the face with respect to which we are going to build the derivative
- auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
- auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
+ auto origFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
+ auto& curFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
for(auto pvIdx : PriVarIndices(faceIdx))
{
PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(curSol[faceIdx][globalJ]));
+
+ auto faceSolution = StaggeredFaceSolution<TypeTag>(scvfJ, curSol[faceIdx], assembler.fvGridGeometry());
+
+
const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, faceIdx);
priVars[pvIdx] += eps;
- curFaceVars.update(priVars[faceIdx]);
+
+ faceSolution[globalJ][pvIdx] += eps;
+ curFaceVars.update(scvfJ, faceSolution);
const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry,
prevElemVolVars, curElemVolVars,
@@ -507,16 +509,17 @@ static void dFacedFace_(Assembler& assembler,
for(const auto& globalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
{
// get the faceVars of the face with respect to which we are going to build the derivative
- auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
- auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
+ auto origFaceVars = curGlobalFaceVars.faceVars(scvf.index());
+ auto& curFaceVars = curGlobalFaceVars.faceVars(scvf.index());
for(auto pvIdx : PriVarIndices(faceIdx))
{
- PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(curSol[faceIdx][globalJ]));
+ auto faceSolution = StaggeredFaceSolution<TypeTag>(scvf, curSol[faceIdx], assembler.fvGridGeometry());
- const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, faceIdx);
- priVars[pvIdx] += eps;
- curFaceVars.update(priVars[faceIdx]);
+ const Scalar eps = numericEpsilon(faceSolution[globalJ][pvIdx], faceIdx, faceIdx);
+
+ faceSolution[globalJ][pvIdx] += eps;
+ curFaceVars.update(scvf, faceSolution);
const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
prevElemVolVars, curElemVolVars,
diff --git a/dumux/discretization/staggered/facesolution.hh b/dumux/discretization/staggered/facesolution.hh
new file mode 100644
index 0000000000..386dbbece0
--- /dev/null
+++ b/dumux/discretization/staggered/facesolution.hh
@@ -0,0 +1,93 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The global volume variables class for cell centered models
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_FACE_SOLUTION_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_FACE_SOLUTION_HH
+
+#include <dumux/common/basicproperties.hh>
+#include <dumux/discretization/staggered/elementvolumevariables.hh>
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class StaggeredFaceSolution
+{
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+ using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+
+ StaggeredFaceSolution(const SubControlVolumeFace& scvf, const FaceSolutionVector& sol,
+ const FVGridGeometry& fvGridGeometry)
+ {
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ const auto& stencil = connectivityMap(faceIdx, faceIdx, scvf.index());
+
+ facePriVars_.clear();
+ map_.clear();
+
+ for(const auto dofJ : stencil)
+ {
+ map_.push_back(dofJ);
+ facePriVars_.push_back(sol[dofJ]);
+ }
+ }
+
+ //! bracket operator const access
+ template<typename IndexType>
+ const FacePrimaryVariables& operator [](IndexType globalFaceDofIdx) const
+ {
+ const auto pos = std::find(map_.begin(), map_.end(), globalFaceDofIdx);
+ assert (pos != map_.end());
+ return facePriVars_[pos - map_.begin()];
+ }
+
+ //! bracket operator
+ template<typename IndexType>
+ FacePrimaryVariables& operator [](IndexType globalFaceDofIdx)
+ {
+ const auto pos = std::find(map_.begin(), map_.end(), globalFaceDofIdx);
+ assert (pos != map_.end());
+ return facePriVars_[pos - map_.begin()];
+ }
+
+
+private:
+
+ std::vector<FacePrimaryVariables> facePriVars_;
+ std::vector<unsigned int> map_;
+
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index 7f8a7b0022..e8bdf8b5d3 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -27,25 +27,105 @@
namespace Dumux
{
+
+namespace Properties
+{
+ NEW_PROP_TAG(StaggeredFaceSolution);
+}
+
template<class TypeTag>
class StaggeredFaceVariables
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FaceSolution = typename GET_PROP_TYPE(TypeTag, StaggeredFaceSolution);
+
+ struct SubFaceData
+ {
+ Scalar velocityNormalInside;
+ Scalar velocityNormalOutside;
+ Scalar velocityParallelInside;
+ Scalar velocityParallelOutside;
+ };
+
public:
void update(const FacePrimaryVariables &facePrivars)
{
- velocity_ = facePrivars[0];
+ velocitySelf_ = facePrivars;
+ }
+
+ void update(const SubControlVolumeFace& scvf, const FaceSolutionVector& sol)
+ {
+ velocitySelf_ = sol[scvf.dofIndex()];
+ velocityOpposite_ = sol[scvf.dofIndexOpposingFace()];
+
+ subFaceVelocities_.resize(scvf.pairData().size());
+
+
+ for(int i = 0; i < scvf.pairData().size(); ++i)
+ {
+ subFaceVelocities_[i].velocityNormalInside = sol[scvf.pairData(i).normalPair.first];
+
+ if(scvf.pairData(i).normalPair.second >= 0)
+ subFaceVelocities_[i].velocityNormalOutside = sol[scvf.pairData(i).normalPair.second];
+
+ subFaceVelocities_[i].velocityParallelInside = sol[scvf.dofIndex()];
+
+ if(scvf.pairData(i).outerParallelFaceDofIdx >= 0)
+ subFaceVelocities_[i].velocityParallelOutside = sol[scvf.pairData(i).outerParallelFaceDofIdx];
+ }
+
+ }
+
+ void update(const SubControlVolumeFace& scvf, const FaceSolution& faceSol)
+ {
+ velocitySelf_ = faceSol[scvf.dofIndex()];
+ velocityOpposite_ = faceSol[scvf.dofIndexOpposingFace()];
+
+ subFaceVelocities_.resize(scvf.pairData().size());
+
+ for(int i = 0; i < scvf.pairData().size(); ++i)
+ {
+ subFaceVelocities_[i].velocityNormalInside = faceSol[scvf.pairData(i).normalPair.first];
+
+ if(scvf.pairData(i).normalPair.second >= 0)
+ subFaceVelocities_[i].velocityNormalOutside = faceSol[scvf.pairData(i).normalPair.second];
+
+ subFaceVelocities_[i].velocityParallelInside = faceSol[scvf.dofIndex()];
+
+ if(scvf.pairData(i).outerParallelFaceDofIdx >= 0)
+ subFaceVelocities_[i].velocityParallelOutside = faceSol[scvf.pairData(i).outerParallelFaceDofIdx];
+ }
+ }
+
+ Scalar velocitySelf() const
+ {
+ return velocitySelf_;
+ }
+
+ Scalar velocityOpposite() const
+ {
+ return velocityOpposite_;
+ }
+
+ auto& subFaceData() const
+ {
+ return subFaceVelocities_;
}
- Scalar velocity() const
+ auto& subFaceData(const int localSubFaceIdx) const
{
- return velocity_;
+ return subFaceVelocities_[localSubFaceIdx];
}
private:
Scalar velocity_;
+ Scalar velocitySelf_;
+ Scalar velocityOpposite_;
+ std::vector<SubFaceData> subFaceVelocities_;
};
} // end namespace
diff --git a/dumux/discretization/staggered/globalfacevariables.hh b/dumux/discretization/staggered/globalfacevariables.hh
index 7582a3bc26..0c0caa7234 100644
--- a/dumux/discretization/staggered/globalfacevariables.hh
+++ b/dumux/discretization/staggered/globalfacevariables.hh
@@ -24,6 +24,7 @@
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FACEVARIABLES_HH
#include <dumux/common/basicproperties.hh>
+#include <dumux/discretization/staggered/facesolution.hh>
namespace Dumux
{
@@ -51,15 +52,36 @@ public:
const auto& faceSol = sol[faceIdx];
- const auto numFaceDofs = fvGridGeometry.gridView().size(1);
- faceVariables_.resize(numFaceDofs);
- assert(faceVariables_.size() == faceSol.size());
- for(int i = 0; i < numFaceDofs; ++i)
+ // const auto numFaceDofs = fvGridGeometry.gridView().size(1);
+
+ // faceVariables_.resize(numFaceDofs);
+ faceVariables_.resize(fvGridGeometry.numScvf());
+
+ for(auto&& element : elements(fvGridGeometry.gridView()))
{
- faceVariables_[i].update(faceSol[i]);
+ auto fvGeometry = localView(fvGridGeometry);
+ fvGeometry.bindElement(element);
+
+ for(auto&& scvf : scvfs(fvGeometry))
+ {
+ faceVariables_[scvf.index()].update(scvf, faceSol);
+ auto test = StaggeredFaceSolution<TypeTag>(scvf, faceSol, fvGridGeometry);
+ // std::cout << "test" << test[scvf.dofIndex()] << std::endl;
+ }
}
+
+ // for(auto& faceVariable : faceVariables_)
+ // {
+ // faceVariable.update()
+ // }
+ // // assert(faceVariables_.size() == faceSol.size());
+ //
+ // for(int i = 0; i < numFaceDofs; ++i)
+ // {
+ // faceVariables_[i].update(faceSol[i]);
+ // }
}
const FaceVariables& faceVars(const IndexType facetIdx) const
diff --git a/dumux/discretization/staggered/properties.hh b/dumux/discretization/staggered/properties.hh
index 5ea414b538..67cb5b475c 100644
--- a/dumux/discretization/staggered/properties.hh
+++ b/dumux/discretization/staggered/properties.hh
@@ -44,6 +44,7 @@
#include <dumux/discretization/staggered/fvgridgeometry.hh>
#include <dumux/discretization/staggered/fvelementgeometry.hh>
#include <dumux/discretization/staggered/globalfacevariables.hh>
+#include <dumux/discretization/staggered/facesolution.hh>
#include <dumux/common/intersectionmapper.hh>
#include <dune/istl/multitypeblockvector.hh>
@@ -60,6 +61,7 @@ namespace Properties
NEW_PROP_TAG(CellCenterSolutionVector);
NEW_PROP_TAG(FaceSolutionVector);
+NEW_PROP_TAG(StaggeredFaceSolution);
//! Type tag for the box scheme.
NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(FiniteVolumeModel, NumericModel));
@@ -109,6 +111,8 @@ SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, StaggeredLocalResidual<TypeTag>
SET_TYPE_PROP(StaggeredModel, IntersectionMapper, Dumux::ConformingGridIntersectionMapper<TypeTag>);
+SET_TYPE_PROP(StaggeredModel, StaggeredFaceSolution, StaggeredFaceSolution<TypeTag>);
+
//! Definition of the indices for cell center and face dofs in the global solution vector
SET_PROP(StaggeredModel, DofTypeIndices)
{
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index b4d5878016..a3de0a7ccc 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -111,7 +111,8 @@ public:
const auto& outsideVolVars = scvf.boundary() ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
CellCenterPrimaryVariables flux(0.0);
- const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+ // const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+ const Scalar velocity = globalFaceVars.faceVars(scvf.index()).velocitySelf();
const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
@@ -201,8 +202,8 @@ public:
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideVolVars = elemVolVars[insideScvIdx];
- const Scalar velocitySelf = globalFaceVars.faceVars(scvf.dofIndex()).velocity() ;
- const Scalar velocityOpposite = globalFaceVars.faceVars(scvf.dofIndexOpposingFace()).velocity();
+ const Scalar velocitySelf = globalFaceVars.faceVars(scvf.index()).velocitySelf() ;
+ const Scalar velocityOpposite = globalFaceVars.faceVars(scvf.index()).velocityOpposite();
FacePrimaryVariables normalFlux(0.0);
if(navierStokes)
@@ -253,20 +254,24 @@ public:
{
FacePrimaryVariables tangentialFlux(0.0);
- // convenience function to get the velocity on a face
- auto velocity = [&globalFaceVars](const int dofIdx)
- {
- return globalFaceVars.faceVars(dofIdx).velocity();
- };
+ // // convenience function to get the velocity on a face
+ // auto velocity = [&globalFaceVars](const int dofIdx)
+ // {
+ // return globalFaceVars.faceVars(scvf.dofIdx()).velocity();
+ // };
+ auto& faceVars = globalFaceVars.faceVars(scvf.index());
+
+ const int numSubFaces = scvf.pairData().size();
// account for all sub-faces
- for(const auto& subFaceData : scvf.pairData())
+ // for(const auto& subFaceData : scvf.pairData())
+ for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
- const auto& normalFace = fvGeometry.scvf(eIdx, subFaceData.localNormalFaceIdx);
+ const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
// Check if we have a symmetry boundary condition. If yes, the tangental part of the momentum flux can be neglected.
- if(subFaceData.outerParallelFaceDofIdx < 0)
+ if(scvf.pairData()[localSubFaceIdx].outerParallelFaceDofIdx < 0)
{
// lambda to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
auto makeGhostFace = [eIdx] (const GlobalPosition& pos)
@@ -275,32 +280,33 @@ public:
};
// use the ghost face to check if there is a symmetry boundary condition and skip any further steps if yes
- const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos));
+ const auto bcTypes = problem.boundaryTypes(element, makeGhostFace(scvf.pairData()[localSubFaceIdx].virtualOuterParallelFaceDofPos));
if(bcTypes.isSymmetry())
continue;
}
// if there is no symmetry boundary condition, proceed to calculate the tangential momentum flux
if(navierStokes)
- tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, subFaceData, elemVolVars, velocity);
+ tangentialFlux += computeAdvectivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
- tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, subFaceData, elemVolVars, velocity);
+ tangentialFlux += computeDiffusivePartOfTangentialMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx);
}
return tangentialFlux;
}
private:
- template<class SubFaceData, class VelocityHelper>
+ template<class FaceVars>
FacePrimaryVariables computeAdvectivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
- const SubFaceData& subFaceData,
const ElementVolumeVariables& elemVolVars,
- VelocityHelper velocity)
+ const FaceVars& faceVars,
+ const int localSubFaceIdx)
{
- const Scalar transportingVelocity = velocity(subFaceData.normalPair.first);
+ // const Scalar transportingVelocity = velocity(subFaceData.normalPair.first);
+ const Scalar transportingVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
const auto insideScvIdx = normalFace.insideScvIdx();
const auto outsideScvIdx = normalFace.outsideScvIdx();
@@ -317,12 +323,14 @@ private:
Scalar transportedVelocity(0.0);
if(innerElementIsUpstream)
- transportedVelocity = velocity(scvf.dofIndex());
+ transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelInside;
+ // transportedVelocity = velocity(scvf.dofIndex());
else
{
- const int outerDofIdx = subFaceData.outerParallelFaceDofIdx;
+ const int outerDofIdx = scvf.pairData(localSubFaceIdx).outerParallelFaceDofIdx;
if(outerDofIdx >= 0)
- transportedVelocity = velocity(outerDofIdx);
+ // transportedVelocity = velocity(outerDofIdx);
+ transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelOutside;
else // this is the case when the outer parallal dof would lie outside the domain TODO: discuss which one is better
// transportedVelocity = problem.dirichlet(makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
transportedVelocity = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
@@ -334,14 +342,14 @@ private:
return transportingVelocity * momentum * sgn * normalFace.area() * 0.5;
}
- template<class SubFaceData, class VelocityHelper>
+ template<class FaceVars>
FacePrimaryVariables computeDiffusivePartOfTangentialMomentumFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
- const SubFaceData& subFaceData,
const ElementVolumeVariables& elemVolVars,
- VelocityHelper velocity)
+ const FaceVars& faceVars,
+ const int localSubFaceIdx)
{
FacePrimaryVariables tangentialDiffusiveFlux(0.0);
@@ -362,35 +370,38 @@ private:
const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
// the normal derivative
- const int innerNormalVelocityIdx = subFaceData.normalPair.first;
- const int outerNormalVelocityIdx = subFaceData.normalPair.second;
+ // const int innerNormalVelocityIdx = subFaceData.normalPair.first;
+ const int outerNormalVelocityIdx = scvf.pairData(localSubFaceIdx).normalPair.second;
- const Scalar innerNormalVelocity = velocity(innerNormalVelocityIdx);
+ // const Scalar innerNormalVelocity = velocity(innerNormalVelocityIdx);
+ const Scalar innerNormalVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
const Scalar outerNormalVelocity = outerNormalVelocityIdx >= 0 ?
- velocity(outerNormalVelocityIdx) :
- problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterNormalFaceDofPos))[faceIdx][normalDirIdx];
+ faceVars.subFaceData(localSubFaceIdx).velocityNormalOutside :
+ problem.dirichlet(element, makeGhostFace(scvf.pairData(localSubFaceIdx).virtualOuterNormalFaceDofPos))[faceIdx][normalDirIdx];
const Scalar normalDeltaV = scvf.normalInPosCoordDir() ?
(outerNormalVelocity - innerNormalVelocity) :
(innerNormalVelocity - outerNormalVelocity);
- const Scalar normalDerivative = normalDeltaV / subFaceData.normalDistance;
+ const Scalar normalDerivative = normalDeltaV / scvf.pairData(localSubFaceIdx).normalDistance;
tangentialDiffusiveFlux -= muAvg * normalDerivative;
// the parallel derivative
- const Scalar innerParallelVelocity = velocity(scvf.dofIndex());
+ // const Scalar innerParallelVelocity = velocity(scvf.dofIndex());
+ const Scalar innerParallelVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelInside;
- const int outerParallelFaceDofIdx = subFaceData.outerParallelFaceDofIdx;
+ const int outerParallelFaceDofIdx = scvf.pairData(localSubFaceIdx).outerParallelFaceDofIdx;
const Scalar outerParallelVelocity = outerParallelFaceDofIdx >= 0 ?
- velocity(outerParallelFaceDofIdx) :
- problem.dirichlet(element, makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
+ faceVars.subFaceData(localSubFaceIdx).velocityParallelOutside :
+ // velocity(outerParallelFaceDofIdx) :
+ problem.dirichlet(element, makeGhostFace(scvf.pairData(localSubFaceIdx).virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
const Scalar parallelDeltaV = normalFace.normalInPosCoordDir() ?
(outerParallelVelocity - innerParallelVelocity) :
(innerParallelVelocity - outerParallelVelocity);
- const Scalar parallelDerivative = parallelDeltaV / subFaceData.parallelDistance;
+ const Scalar parallelDerivative = parallelDeltaV / scvf.pairData(localSubFaceIdx).parallelDistance;
tangentialDiffusiveFlux -= muAvg * parallelDerivative;
const Scalar sgn = sign(normalFace.outerNormalScalar());
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index f64775a4a7..606f5fe287 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -182,7 +182,7 @@ public:
const GlobalFaceVars& globalFaceVars)
{
FacePrimaryVariables storage(0.0);
- const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+ const Scalar velocity = globalFaceVars.faceVars(scvf.index()).velocity();
storage[0] = volVars.density() * velocity;
return storage;
}
@@ -327,7 +327,8 @@ protected:
// set a fixed value for the velocity for Dirichlet boundary conditions
if(bcTypes.isDirichlet(momentumBalanceIdx))
{
- const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
+ // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
+ const Scalar velocity = faceVars.faceVars(scvf.index()).velocitySelf();
const Scalar dirichletValue = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
residual = velocity - dirichletValue;
}
@@ -336,7 +337,8 @@ protected:
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
if(bcTypes.isSymmetry())
{
- const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
+ // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
+ const Scalar velocity = faceVars.faceVars(scvf.index()).velocitySelf();
const Scalar fixedValue = 0.0;
residual = velocity - fixedValue;
}
diff --git a/dumux/freeflow/staggered/velocityoutput.hh b/dumux/freeflow/staggered/velocityoutput.hh
index d23742c857..e3f7e333f3 100644
--- a/dumux/freeflow/staggered/velocityoutput.hh
+++ b/dumux/freeflow/staggered/velocityoutput.hh
@@ -113,9 +113,9 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
{
- auto& origFaceVars = gridVariables_.curGridFaceVars().faceVars(scvf.dofIndex());
+ auto& origFaceVars = gridVariables_.curGridFaceVars().faceVars(scvf.index());
auto dirIdx = scvf.directionIndex();
- velocity[dofIdxGlobal][dirIdx] += 0.5*origFaceVars.velocity();
+ velocity[dofIdxGlobal][dirIdx] += 0.5*origFaceVars.velocitySelf();
}
}
}
diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh
index fc00686af3..2cf09627ea 100644
--- a/dumux/implicit/staggered/newtoncontroller.hh
+++ b/dumux/implicit/staggered/newtoncontroller.hh
@@ -137,7 +137,7 @@ public:
BlockVector y;
y.resize(numRows);
- // printmatrix(std::cout, M, "", "");
+ // printmatrix(std::cout, M, "old", "");
// solve
const bool converged = ls.solve(M, y, bTmp);
--
GitLab